DEGRADATION OF POLYMERIC ADHESIVES IN ADVERSE ENVIRONMENTAL
CONDITIONS

By

Yang Chen

A DISSERTATION

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

Civil Engineering—Doctor of Philosophy

2025

ABSTRACT

Cross-linked elastomers, known for their exceptional flexibility, toughness, formability, and

versatility, play a vital role in a wide array of engineering applications across aerospace, construc-

tion, transportation, marine, aeronautics, and automotive sectors. These materials are expected to

maintain high performance throughout their service life, even when exposed to aggressive environ-

mental conditions such as water infusion, temperature fluctuations, and ultraviolet (UV) radiation.

These environmental factors pose significant challenges, as they can gradually deteriorate the ma-

terial’s properties and reduce its overall durability.

Among the most critical forms of environmental degradation are thermal aging under oxygen-

deficient conditions and oxidative aging caused by elevated oxygen concentrations. The former

typically results in uniform thermal degradation of the material, while the latter, known as diffusion-

limited oxidation (DLO), induces spatially heterogeneous damage, primarily at the surfaces of the

polymer where oxygen diffusion is more pronounced. Together, these two forms of aging represent

fundamentally inverse degradation conditions—one being inert and volumetric, and the other being

oxidative and surface-driven. Accurate prediction of the long-term behavior of elastomers under

both these conditions is essential for designing reliable rubber components that resist early failure

in service.

To address this need, high-fidelity constitutive models are essential for simulating the effects of

aging on the mechanical, thermal, and failure characteristics of cross-linked polymers. Historically,

most aging models have employed hyperelastic constitutive laws coupled with single-kinetic degra-

dation equations to model the evolution of material properties over time. While these approaches

have been useful in capturing basic degradation trends, they often fall short in representing the

complex interactions between microstructural evolution and macroscale mechanical behavior un-

der realistic service conditions.

This dissertation presents a comprehensive multi-physics modeling framework to capture the

distinct and coupled degradation behaviors of cross-linked polymers under both diffusion-limited

oxidation and inert thermal aging. These models incorporate the effects of oxygen diffusion, re-

action kinetics, and thermally activated chain scission and cross-linking processes to simulate the

evolution of polymer microstructure over time. By resolving the spatial and temporal development

of aging parameters, the framework successfully reproduces both uniform and spatially heteroge-

neous degradation phenomena, providing critical insight into how these opposing environmental

factors influence long-term material performance.

The modeling approach is grounded in continuum mechanics and integrates finite strain the-

ory with micro-mechanically motivated degradation mechanisms. Rubber elasticity is described

across three scales: statistical mechanics at the microscale to account for molecular chain behav-

ior, network-based phenomenological modeling at the mesoscale, and continuum theory at the

macroscale. This multi-scale approach enables the representation of polymer network reconfig-

uration during aging, which is essential for predicting stiffness loss, permanent deformation, and

eventual failure.

While the core of this dissertation focuses on physics-based modeling of aging mechanisms,

recent advancements in machine learning (ML), particularly physics-informed neural networks

(PINNs), offer promising opportunities to enhance traditional modeling techniques. Although early

generations of data-driven black-box models were limited by their need for large datasets and lack

of physical constraints, hybrid approaches—where experimental macroscopic data are used to infer

underlying microstructural behavior within physics-informed frameworks—have begun to bridge

this gap. In this work, such techniques are explored in a supporting capacity to augment the pre-

dictive capabilities of the proposed aging models, without compromising the underlying physical

consistency.

By coupling the two degradation models and incorporating both mechanistic understanding and

data-driven insights, this dissertation delivers a unified computational framework capable of simu-

lating the long-term behavior of elastomeric materials under varied environmental exposures. The

results not only enhance our understanding of how oxidative and inert aging conditions uniquely

and jointly affect polymer durability, but also provide a practical foundation for the design of next-

generation elastomers with improved resistance to environmental degradation.

Copyright by
YANG CHEN
2025

This dissertation is dedicated to my advisor, Dr. Dargazany, for his invaluable guidance.
To my respectful parents, Wuyun Chen and Weiwei Zhang, for always believing in me.
And to my caring and beautiful wife Yiwen Hu, for her unwavering support.

v

ACKNOWLEDGMENTS

I would like to express my heartfelt gratitude to my advisor, Dr. Roozbeh Dargazany, for his

unwavering support. His insightful guidance, invaluable advice, and endless assistance have been

crucial in making this journey possible.

Next, I extend my deepest appreciation to my committee members, Dr. Nizar Lajnef, Dr. Sara

Roccabianca, and Dr. Venkatesh K. R. Kodur. It has been a privilege to benefit from their vast

experiences over the years. The honor of attending their classes and drawing from their expertise

has profoundly influenced my intellectual development as a graduate student.

I am also grateful for the assistance and camaraderie of my dear peers and friends in the High-

Performance Materials (HPM) group. Special thanks to Vahid Morovati, Hamid Mohammadi, Amir

Bahrololoumi, Aref Ghaderi, Sharif Alazhary, and Mamoon Shaafaey for sharing ideas and provid-

ing invaluable help throughout the experimental procedures.

vi

TABLE OF CONTENTS

CHAPTER 1: INTRODUCTION .

.

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

1

CHAPTER 2: CONTINUUM MECHANICS . . . . . . . . . . . . . . . . . . . . . . . .

. 11

CHAPTER 3: CONSTITUTIVE RELATIONS . . . . . . . . . . . . . . . . . . . . . . . . 24

CHAPTER 4: MACHINE-LEARNED CONSTITUTIVE MODE FOR SINGE . . . . . . . 29

CHAPTER 5: FINITE ELEMENT ANALYSIS . . . . . . . . . . . . . . . . . . . . . . . . 46

CHAPTER 6: INERT THERMAL AGING . . . . . . . . . . . . . . . . . . . . . . . . . . 54

CHAPTER 7: CONSTITUTIVE MODELING OF DIFFUSION-LIMITED OXIDATION

(DLO) COUPLED WITH A LARGE DEFORMATION THEORY FOR
POLYMER DEGRADATION . . . . . . . . . . . . . . . . . . . . . . . . . . 79

CHAPTER 8: GRAPHIC USER INTERFACE . . . . . . . . . . . . . . . . . . . . . . .

. 104

BIBLIOGRAPHY .

APPENDIX .

.

.

.

.

.

.

.

.

.

.

. .

. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

. .

. .

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

. 131

vii

CHAPTER 1: INTRODUCTION

1.1: Motivation, Challenges and Objectives

Challenges: Despite the widespread use of polymer adhesives, accurately predicting their long-

term performance remains a significant challenge, especially under harsh environmental condi-

tions. One critical factor is thermal aging, which can drastically alter the mechanical and chemical

integrity of the polymer network. This aging is typically accompanied by competing processes

such as chain scission and additional cross-linking, both of which influence the stiffness, strength,

and durability of the material.

Characterizing and modeling these phenomena is complicated by the nonlinear, time-dependent

behavior of polymers and the interplay between physical and chemical degradation mechanisms.

Moreover, the multiscale nature of the material response — from molecular changes at the chain

level to macroscopic deformation behavior — necessitates a modeling approach that captures both

microstructural evolution and continuum-scale performance.

Traditional constitutive models often fail to account for such micro-mechanical transformations,

limiting their predictive capabilities. In addition, implementing these models within a finite ele-

ment framework poses numerical challenges, particularly when large deformations and nonlinear

elasticity are involved.

Objectives: The primary objective of this research is to develop a mechanistically informed, micro-

mechanical model that captures the constitutive behavior of cross-linked polymer adhesives under

inert thermal aging. Specifically, the model aims to:

• Capture the competitive and time-dependent effects of chain scission and cross-linking dur-

ing thermal exposure;

• Relate molecular-scale damage mechanisms to macroscopic mechanical response;

• Implement the model within a finite element framework under plane stress conditions, incor-

porating consistent tangent stiffness for large deformation analysis;

• Validate the model using experimental data and numerical simulations, including compar-

1

isons with custom UMAT subroutines and continuum solvers;

• Provide insights to guide the design and life prediction of adhesive joints exposed to elevated

temperatures in oxygen-limited environments.

By addressing these objectives, this work seeks to bridge the gap between chemical degradation

processes and mechanical modeling, enabling a more accurate prediction of adhesive behavior over

the service life of engineered systems.

1.2: Fundamentals of Cross-linked Polymers

1.2.1: Fundamentals of Polymer Physics

The mechanical behavior of polymers is highly sensitive to temperature and the time scale of de-

formation. Polymers are inherently viscoelastic, exhibiting characteristics of both viscous fluids

and elastic solids. At low temperatures or high loading rates, polymers behave in a glassy manner,

often fracturing or deforming irreversibly under modest strains. At elevated temperatures or low

strain rates, however, they can adopt a rubber-like response, sustaining large deformations (stretch

ratios exceeding 2) without permanent set. Further increases in temperature result in behavior akin

to viscous liquids, where deformation becomes irreversible and time-dependent.

Phenomenological models typically describe material response at a macroscopic level without

incorporating explicit microstructural details. These models are often derived by: (i) proposing

plausible strain energy functions, or (ii) extending small-strain elasticity theory to accommodate

large strain regimes.

For hyperelastic rubber-like materials, two key assumptions are commonly made: (i) The mate-

rial is isotropic in its undeformed state, exhibiting identical behavior in all directions, and (ii) The

material is nearly incompressible, with negligible volume change upon deformation. Accordingly,

modeling efforts frequently focus on homogeneous pure strains under these assumptions.

2

Figure 1.1: The Mullins effect with permanent set: (a) Two uniaxial tensile loading cycles; (b)
Corresponding stress-stretch curves.

Cross-linked polymers generally exhibit pronounced nonlinear elasticity coupled with inelas-

tic effects such as stress-softening (Mullins effect). Amorphous elastomers in particular demon-

strate rate-dependent elastic-plastic responses due to their complex microstructure, as illustrated

in Fig. 1.2. These behaviors arise from molecular mechanisms such as chain rearrangement, kink

unfolding, and segment reorientation during loading.

Figure 1.2: Schematic representation of the microstructure in filled elastomers.

Micro-stretch λd:

In micromechanical network models, polymer chains oriented in direction d experience direc-

tional stretching governed by the deformation gradient F. The chain-level (micro) stretch is defined

as:

λd = 1 +

∆Ld
Ld

,

where ∆Ld is the elongation, and Ld is the initial chain length. The maximum sustainable micro-

3

123341stloading2nd loadingunloadingunloading𝜆𝑟𝑒(1) 𝜆𝑟𝑒(2) StressStretcha)b)𝜆𝑟𝑒(1) 𝜆𝑟𝑒(2) 𝜆(2) 𝜆(1) HCrosslinkerCH=CH2Polymer ChainChemical crosslinks CH2=CH2FillerCyclic siloxanesSiloxanes (low molecular weight)Polysiloxanes (medium molecular weight)stretch for a chain is limited by its geometry:

λd
limit =

Lc
νR0

.

Here, Lc is the contour length, R0 is the end-to-end distance, and ν > 1 is the threshold for me-

chanical rupture [31]. Shorter chains are more susceptible to breakage due to higher stretch ratios.

Mechanisms of Breakage:

Chain failure can occur in various forms, including polymer-polymer detachment, cross-link

rupture, desorption from filler surfaces, and chain slippage (see Fig.1.3). These microstructural

mechanisms play a central role in damage evolution and are essential for understanding the Mullins

effect and stress-softening behavior.

History Dependence via Maximum Micro-stretch λd

max:

During cyclic loading, polymer chains progressively debond from aggregates, starting with the

shortest chains. Once debonded, chains do not reattach, making the history of deformation critical.

The maximum micro-stretch in a given direction d over the material’s loading history is defined as:

λd
max = max

τ ∈(−∞,t]

λd(τ ).

This parameter serves as a key internal variable for modeling direction-dependent damage.

4

Figure 1.3: Mechanisms of chain damage in cross-linked polymers: (a) Cross-link breakage; (b)
Chain scission; (c) Rupture of rubber–filler bonds; (d) Desorption from filler surface; (e) Chain
slippage at entanglement sites; (f) Slippage over particle surfaces; (g) Re-adsorption onto particle
surfaces; (h) Filler rupture.

1.3: Non-linear Features in Cross-Linked Polymers

Cross-linked polymers are commonly modeled as hyperelastic materials, characterized by their

ability to undergo large, reversible deformations. However, their mechanical behavior also includes

5

BreakageBreakageLonger chainBreakageBreakageLonger chaindesorptionEntangelmentLonger chainShorter chainShorter chainLonger chainSlippageabsorptionRuptura)b)c)d)e)f)g)h)Figure 1.4: Schematic illustration of deformation-induced damage mechanisms in cross-linked
polymers.

prominent inelastic features, particularly under cyclic or large-strain loading. One notable inelastic

phenomenon is the Mullins effect, which refers to the stress-softening observed after the initial

loading cycle [10, 85]. This effect is observed in both filled and unfilled cross-linked polymers.

To better understand the microstructural origins of such inelastic effects, Fig.1.4 illustrates sev-

eral deformation-induced damage mechanisms in amorphous polymer networks. These include:

• Chain scission [13]

• Chain disentanglement [53]

• Molecular slippage at junctions [62]

• Filler–matrix debonding and rupture of filler clusters [72]

These mechanisms collectively lead to permanent changes in the microstructure of the polymer

network. As a result, materials often exhibit a residual strain upon unloading—referred to as the

permanent set. While the permanent set is minimal in unfilled elastomers, it becomes significant

in filled systems due to interfacial debonding and filler rupture, as depicted in Fig.1.4.

Traditional modeling of such nonlinear and inelastic behaviors has relied on two main ap-

6

deformationdeformationdeformationdeformationBond ruptureMolecules slippingFiller ruptureDisentanglement(a)(b)(c)(d)StretchReloadingLoadingStress(e)Mullins EffectPermanent SetMullins effect + Permanent set proaches:

• Phenomenological models, which offer computational efficiency but lack physical inter-

pretability and generalizability outside calibration conditions.

• Physics-based models, which provide greater fidelity but are often too complex for real-time

or large-scale simulations.

To bridge this gap, we propose a hybrid framework that couples physically informed micro-

mechanical theory with a machine-learned agent framework. Specifically, we develop a knowledge-

driven ML model where internal variables evolve based on known physical laws and history-

dependent damage. The learning agents (L-agents) are designed with memory structures to reflect

material behavior, particularly relevant in rubbery polymers where damage evolution depends on

the maximum historical strain rather than the full deformation path.

A microsphere-based formulation is adopted to account for directional dependence in the poly-

mer matrix, enabling accurate modeling of inelastic effects such as stress-softening, hysteresis, and

permanent set. The architecture of the neural network—including activation functions, number of

hidden layers, and neurons—is tailored to balance complexity and accuracy. These modeling strate-

gies are detailed in the following section.

1.4: Frame Independency

Frame objectivity, during rigid body motion, requires strain energy of the material remains un-

changed. Thus, the material response should not depend on the choice of the reference frame. The

strain energy frame independency can be written as

Ψm(QF) = Ψm(F),

(1.1)

where Q is the rotation tensor. So, a constitutive law is frame independent if energy is left

rotationally invariant. The mentioned condition is satisfied when the strain energy is a function of

the right Cauchy-Green deformation tensor C, due to

7

C+ = (F+)T F+ = FT QT QF = FT F = C,

(1.2)

which F+ = FQ. The proposed model is a function of right Cauchy-Green deformation tensor. So,

the frame independency condition is satisfied automatically.

1.5: Thermodynamic Consistency

1.5.1: Polyconvexity

Polyconvexity is one of the known conditions which ensure the thermodynamic consistency. In

this section, we briefly describe sufficient but not necessary free energy function conditions which

guarantee the existence of minimizers of some variational principles. In order to understand poly-

convexity, we start with some properties of convexity. Consider that Ψm (F) is the strain energy

function on set of K. We can say Ψm (F) is convex on set of K if hessian matrix of Ψm (F) be

positive in that set.

D2Ψm (F). (H, H) ⩾ 0,

(1.3)

and for proof of polyconvexity we can mention that F → Ψm (F) is polyconvex if and only if

there exist a function G such that

Ψm(F) = G(F, adjF, det F),

(1.4)

and the function G is convex. Besides, F = F−1

detF and the implication chain shows relations from

convexity to ellipticity.

convexity → polyconvexity → quasiconvexity → ellipticity.

The Hessian matrix of the strain energy is positive if

8

2

∂Ψm
di
∂2λj

NdX

=

wi

i=1

∂ψj
∂2λj

2

NdX

di

=

di

wi

i=1

∂AN Nj(Wj, λj
di
∂λj

⃗di, λj max

⃗di)

> 0,

f or j = 1, 2, ..., Ns.

(1.5)

If weights which connect the input of λj to other neurons be positive, the proposed model holds

the condition of polyconvexity.

1.5.2: Second Law of Thermodynamic

Because all of the constitutive models should satisfy the second law of thermodynamic, the sat-

isfaction of this law should be checked for the proposed model. On the other hand, checking

Clausius-Duhem inequality would be enough for this. Because λj max is internal variables in the

strain energy function of Cross-linked polymers, we can reduce the second law of thermodynamics

to Clausius-Duhem inequality that shows thermodynamic consistency of the model in direction di.

This inequality can be written as

∂Ψm
∂λjmax

di

≤ 0

∀ d

f or j = 1, 2, ..., Ns.

(1.6)

If we consider the energy of matrix as

which

NdX

NsX

(cid:0)

(cid:1)

di

ψj

wi,

Ψm =

i=1

j=1

d
ψj = AN Nj(Wj, λj

⃗di, λj max

⃗di),

(1.7)

(1.8)

thus, Clausius-Duhem can be written as

9

∂Ψm
∂λj max

di

=

NdX

i=1

wi

di

∂ψj
∂λj max

=

di

NdX

i=1

wi

∂AN Nj(Wj, λj
∂λj max

di

⃗di, λj max

⃗di)

≤ 0,

f or j = 1, 2, ..., Ns.

(1.9)

If weights that connect the input of λj max to other neurons be negative, the proposed model holds

the condition of thermodynamic consistency.

10

CHAPTER 2: CONTINUUM MECHANICS

2.1: Introduction

Continuum mechanics serves as a foundational framework for nonlinear finite element analysis,

focusing on modeling solids and fluids with properties and responses characterized by smooth

functions of spatial variables. This approach overlooks inhomogeneities such as molecular, grain,

or crystal structures.

Consider a body in its initial state at time t = 0; the domain of the body in this initial state

is denoted as Ω0 and referred to as the initial configuration. When describing the motion and

deformation of the body, a configuration is essential as a reference for various equations; this is

termed the reference configuration.

Figure 2.1: A schematic figure of the hierarchical micro-macro multi-scale approach.

In various instances, it becomes necessary to define a configuration known as the undeformed

configuration, occupying Ω0. This undeformed configuration is synonymous with the initial con-

dition. The domain of the body’s current configuration is denoted as Ω and is often referred to as

the deformed configuration.

11

ΓuΓtInitial configurationtuFDeformed configurationC, SxxxxxxxxxxxxxxxxGaussIntegrationpointMulti-scale model2.2: Eulerian and Lagrangian Coordinates

The position vector of a material point in the initial configuration is denoted as X, where

nSDX

X = Xiei =

Xiei,

i=1

(2.1)

where nSD stands for the number of space dimensions, Xi represents the components of the posi-

tion vector in the initial configuration, and ei are the unit base vectors of a rectangular Cartesian

coordinate system. The vector variable X for a specific material point remains constant over time;

these variables X are termed material coordinates or Lagrangian coordinates. The position of a

point in the current configuration is expressed as

nSDX

x = xiei =

xiei,

i=1

(2.2)

where xi are the components of the position vector in the current configuration.

2.2.1: Eulerian and Lagrangian Description

Two methods are employed to depict the deformation and response of a continuum. In the first

method, the material coordinates Xi and time t serve as the independent variables, represented

as x = ￿(X, t); this method is known as the material description or Lagrangian description. In

the second method, the independent variables are the spatial coordinates x and time t. This is

referred to as the spatial or Eulerian description. Distinct symbols are used to differentiate between

independent variables.

2.2.2: Deformation

Consider a continuum body B existing in three-dimensional Euclidean space at a given time t0, as

shown in Fig.2.2. Using the notation X ∈ E3, any point P0 on body B can be uniquely specified

relative to any basis. As the body B moves through three-dimensional Euclidean space from time

t0 to time t, it occupies an alternate geometric configuration referred to as β’s configuration. Since

12

the mapping from body B to β is one-to-one, point P0 on body B corresponds to point P .

To simplify the understanding of the mapping process, we label the regions of body B and β

as the reference and current configurations, respectively. The positions of points P and P0 can be

expressed using coordinates in a system denoted by θi (where i = 1, 2, 3).

(cid:0)

x = ˆx

θ1, θ2, θ3, t

(cid:1)

, X = ˆX

(cid:0)

θ1, θ2, θ3

(cid:1)

(cid:0)

= ˆx

θ1, θ2, θ3, t0

(cid:1)

,

i = 1, 2, 3,

(2.3)

respectively. As a result, the displacement of point P may be calculated by

(cid:0)

u = ˆu

θ1, θ2, θ3, t

(cid:1)

= x − X.

(2.4)

Here, u represents the displacement vector. When we introduce Euclidean space with a set of

pairwise orthogonal unit vectors, denoted as ei, we have

X = X iei, X j = X · ej,

j = 1, 2, 3,

u = uiei,

uj = u · ej,

j = 1, 2, 3,

(2.5)

x = xiei,

xj = x · ej = X j + uj,

j = 1, 2, 3,

where the Einstein summation convention is used to repeated indices.

Figure 2.2: Description of body motion.

Body deformations can be easily described by vectors tangent to the coordinate lines in both the

13

reference and current configurations, particularly in the context of significant deformations. These

tangent vectors are defined as follows

Gi =

∂X
∂θi ,

gi =

∂x
∂θi ,

i = 1, 2, 3.

Furthermore, Gi and gi may be defined as the bases dual to Gi and gi, implying that

Gi · Gj = δj
i ,

and

gi

· gj = δj
i .

Also, the deformation gradient tensor F can be represented by,

F = Gradx,

where F can be written in a matrix format in two dimensions as

2

6
4

∂x1
∂X1

∂x1
∂X1

3

7
5 =

2

6
4

∂x
∂X

∂y
∂X

∂x1
∂X1

∂x1
∂X1

3

7
5 .

∂x
∂Y

∂y
∂Y

F =

(2.6)

(2.7)

(2.8)

(2.9)

In deformation kinematics, the deformation gradient tensor is a fundamental second-order tensor

that describes the changes in material constituents during deformation. To this goal, consider dX

and dx in reference and current configurations as infinitesimal material components, respectively.

Then, one has

dx = FdX,

dX = F−1dx.

Additionally, the changes in volume J may be derived by

J =

dV
dV0

(cid:12)
(cid:12)F i
·j

(cid:12)
(cid:12) = det F > 0.

=

In the current setup, the lengths of the material components are provided by

(2.10)

(2.11)

14

∥dx∥2 = dx · dx = dX
(cid:0)

dX = dXCdX,
(cid:1)

∥dX∥2 = dX · dX = dx

F−TF−1

dx = dxb−1dx,

(cid:0)

(cid:1)
FTF

where

C = FTF = gijGi ⊗ Gj,

b = FT = Gijgi

⊗ gj,

(2.12)

(2.13)

are referred to as the right and left Cauchy-Green tensor, respectively. The principal invariant of C

mentioned as

I1 (C) = tr (C) ,

I2 (C) =

1
2

((I1(C))2 − tr(C2)),

I3 (C) = det (C) .

(2.14)

To ascertain the stretch of a material element, divide the deformed length by the referenced length

of the material element. To achieve this, the stretch in the direction of N is obtained by summing

the unit vectors N and N along the material element dX and its counterpart dx.

s

λ(N) =

dx
dX

=

∥dx∥2
∥dX∥2 =

r

dXNCNdX
dX 2

= (NCN)

1
2 .

(2.15)

2.3: Strain Measures

Unlike in linear elasticity, nonlinear continuum mechanics employs various measures of strain and

strain rate. However, only two of these measures are taken into consideration here.

1. The Green (Green-Lagrange) strain, E.

2. The rate-of-deformation tensor, D.

A strain measure should exhibit zero values for any rigid body motion, particularly for rigid body

rotation. If a strain measure does not satisfy this condition, it would result in predicting non-zero

strains and consequently non-zero stresses during rigid body rotation. This fundamental reason

underlies the abandonment of the conventional linear strain displacement equations in nonlinear

theory.

15

2.3.1: Green Strain Tensor

The Green strain tensor E is defined by

ds2 − dS2 = 2dX · E · dX.

(2.16)

The Green strain is defined as the disparity between the square of an infinitesimal segment in

the current (deformed) configuration and its counterpart in the initial (undeformed) configuration.

This concept becomes more evident when expressed in indicial notation.

dx · dx = dxidxi = FijdXjFikdXk = dX · (FT · F) · dX.

(2.17)

As dX · FT · F · dX − dX · I · dX − dX · 2E · dX = 0, it follows that

E =

(cid:1)

FT · F − I

(cid:0)

1
2

or Eij =

ikFkj − δij
F T

(cid:1)

.

(cid:0)

1
2

(2.18)

2.3.2: Rate-of-Deformation

The second kinematic measure to be discussed is the rate-of-deformation tensor D. Unlike the

Green strain tensor, it serves as a rate measure of deformation. The velocity gradient is defined as:

L =

∂v
x

= (gradv)T or Lij =

∂vi
∂xj

.

(2.19)

The velocity gradient tensor can be decomposed into symmetric and skew-symmetric components

through the following

L =

(cid:0)

1
2

L + LT

(cid:1)

+

(cid:0)

1
2

(cid:1)

L − LT

or Lij =

1
2

(Lij + Lji) +

1
2

(Lij − Lji) .

(2.20)

16

The rate-of-deformation tensor D is defined as the symmetric part of L, while the spin tensor W is

the skew-symmetric part of L. With these definitions, we can express

L = D + W or Lij = vi,j = Dij + Wij

D =

(cid:1)

L + LT

(cid:0)

1
2

or Dij =

(cid:18)

1
2

(cid:19)

+

∂vj
∂xi

∂vi
∂xj
(cid:1)

W =

(cid:0)

1
2

L − LT

or Wij =

(cid:19)

(cid:18)

1
2

∂vi
∂xj

− ∂vj
∂xi

.

(2.21)

For deriving a unified expression connecting these two measures of strain rate, we have

D =

(cid:0)

1
2

L + LT

(cid:1)

=

(cid:0)

1
2

˙F · F−1 + F−T · ˙FT

(cid:1)

.

Taking time derivative of the Green strain, we got

˙E =

(cid:0)

1
2

D
Dt

FT · F − I

(cid:1)

=

1
2

(cid:0)

(cid:1)
tensF · ˙F + ˙FT · F

.

The above can easily yield

(2.22)

(2.23)

D = F−T · ˙E · F−1 or Dij = F −T

ij

˙EklF −1

lj

.

(2.24)

The time integral of the Green strain rate is path-independent, whereas the time integral of the

rate-of-deformation is not path-independent.

2.4: Stress Measures

In nonlinear problems, different stress measures can be defined. Three stress measures are under

consideration.

1. The Cauchy stress, σ

2. The nominal stress tensor, P, closely related to the first Piola-Kirchhoff stress

3. The second-Piola Kirchhoff stress tensor, S.

17

The stress are defined by Cauchy’s law

n · σdΓ = tdΓ,

where t is the traction. In the initial configuation the counterpart is

n0 · PdΓ0 = t0dΓ.

(2.25)

(2.26)

Please note that the normal is consistently oriented to the left. The definition of the second-Piola

Kirchhoff stress is given by

n0 · SdΓ0 = F−1 · t0dΓ0.

(2.27)

Cauchy stress σ

Nominal stress P 2nd Piola-Kirchhoff stress S

σ
P = JF−1 · σ
S = JF−1 · σ · F−T
τ = Jσ

J−1F · P
-
P · F−T
F · P

J−1F · S · F
S · FT
-
F · S · FT

2.5: Conservation Laws

After formulating the kinematics to represent potential deformed states of a continuous medium,

it’s crucial to note that these fields alone are insufficient for predicting the specific configuration

a body will assume under a given applied load. Achieving such predictions requires an extension

of the principles of mechanics to encompass a continuum medium, coupled with the application

of thermodynamic laws. In this context, we are considering four conservation laws relevant to

thermomechanical systems.

Gauss’s Theorem: In the derivation of conservation equations, Gauss’s theorem is a frequently

employed tool. This theorem establishes a connection between integrals over a domain and an

integral over the boundary of that domain. Specifically, Gauss’s theorem states that when f (x) is

18

piecewise continuously differentiable—a C 0 function then

Z

∂f (x)
∂xi

Ω

Z

dΩ =

nif (x) dΓ.

Γ

(2.28)

The theorem remains applicable to any domain, including the reference domain, where for a C 0

function f (X), we have

Z

∂f (X)
∂xi

Ω0

dΩ0 =

Z

Γ0

n0
i f (X) dΓ0.

(2.29)

Reynolds’ Transport Theorem: The material time derivative of an integral signifies how fast

the integral changes over a material domain. This domain moves with the material, ensuring that

points on its boundary remain fixed on the boundary, with no mass flux across the boundaries. The

various formulations for material time derivatives of integrals are referred to as Reynolds’ transport

theorem. The definition of the material time derivative of an integral is

Z

D
Dt

f dΩ = lim
∆t→0

Ω

1
∆t

 Z

Z

!

f (x, τ + ∆t)dΩ −

f (x, τ )dΩ

,

(2.30)

Ωτ +∆t

Ωτ

here, Ωτ represents the spatial domain at time τ , and Ωτ +∆t denotes the spatial domain occupied by

the same material points at time τ + ∆t. The derivation of Reynolds’ transport theorem is achieved

through the definition of the material time derivative, yielding

Z

D
Dt

Z

(cid:18)

f dΩ =

Ω

Ω

∂f
∂t

+ vi

∂f
∂xi

+

∂vi
∂xif

(cid:19)

Z

(cid:18)

dΩ =

Ω

∂f
∂t

+

∂(vif )
∂xi

(cid:19)

,

(2.31)

which can be written in tensor form as

Z

D
Dt

Z

(cid:18)

f dΩ =

Ω

Ω

∂f
∂t

(cid:19)

+ div(vf )

dΩ.

(2.32)

Mass conservation: A fundamental principle in classical mechanics asserts that mass is an invari-

ant quantity, immune to creation or annihilation; instead, it can only experience deformation in

19

response to applied loads. The mass m(Ω) of a material domain Ω is determined by

Z

m (Ω) =

ρ (x, t) dΩ,

Ω

(2.33)

here, ρ(x, t) represents the density. Mass conservation dictates that the mass of any material domain

remains constant, as no material flows through the boundaries of a material domain, and there is

no consideration for mass-to-energy conversion. Following the principle of mass conservation, the

material time derivative of m(Ω) is zero, indicating

Dm
Dt

=

D
Dt

Z

Ω

ρdΩ = 0.

(2.34)

Applying Reynolds’s theorem, Eqn.2.34 yields

Z

(cid:18)

=

Ω

Dρ
Dt

(cid:19)

+ ρdiv(v)

dΩ = 0.

(2.35)

Since this holds for any subdomain Ω, it follows that

Dρ
Dt

+ ρdiv(v) = 0; or

Dρ
Dt

+ ρvi,i = 0 or ˙ρ + ρvi,i = 0.

(2.36)

The aforementioned equation is the expression of mass conservation, commonly referred to as the

continuity equation. It stands as a first-order partial differential equation.

Conservation of Linear Momentum: The equation derived from the principle of linear momen-

tum conservation holds significance in nonlinear finite element procedures. Linear momentum

conservation aligns with Newton’s second law of motion, linking the forces acting on a body to its

acceleration. We consider an arbitrary domain Ω with a boundary Γ, subjected to body forces ρb

and surface reactions t (force per unit area). The total force is determined by

Z

Z

f (t) =

ρb(x, t)dΩ =

t(x, t)dΓ.

Ω

Γ

(2.37)

20

The linear momentum is given by

Z

p(t) =

ρv(x, t)dΩ,

Ω

(2.38)

where ρv is the linear momentum per unit volume. Using the material time derivative of the linear

momentum, this gives

Dp
Dt

= f ⇒ D
Dt

Z

Ω

Z

Z

ρvdΩ =

ρbdΩ +

tdΓ.

Ω

Γ

(2.39)

Reynolds’ transport theorem applied to LHS integral gives

Z

Ω

D
Dt

ρvdΩ =

Z

(cid:18)

Ω

D
Dt

(cid:19)

Z

(cid:20)

(ρv) + div(v)ρv

=

Ω

ρ

Dv]
Dt

+ v

(cid:18)

Dρ
Dt

(cid:19)(cid:21)

+ ρdiv(v)

dΩ,

(2.40)

where the second equality is obtained by using the product rule of derivatives for the first term

of the integrand and rearranging terms. To convert it into a domain integral, we invoke Cauchy’s

relation and Gauss’s theorem in sequence, giving

Z

Z

Z

tdΓ =

n · σdΓ =

∇ · σdΩ or

Γ

Γ

Ω

Z

Γ

Z

tjdΓ =

niσijdΓ =

Γ

Z

∂σij
∂xi

Ω

dΩ.

(2.41)

In nonlinear continuum mechanics, unlike in linear continuum mechanics, it is crucial to develop

the habit of placing the divergence operator appropriately. This is because certain stress tensors,

like the nominal stress, lack symmetry. In cases where the stress is asymmetric, the left and right

divergence operators yield distinct results. This distinction is illustrated by

Z

(cid:18)

Ω

ρ

Dv
Dt

(cid:19)

− ρb − ∇ · σ

dΩ.

(2.42)

Therefore, if the integrand C −1, since upper equation holds for an arbitrary domain, it yields

ρ

Dv
Dt

= ∇ · σ + ρb ≡ divσ + ρb or ρ

Dvi
Dt

=

∂σij
∂xj

+ ρbi.

(2.43)

21

This is referred to as the momentum equation, also known as the balance of linear momentum.

Equilibrium Equation: In numerous scenarios, loads are applied gradually, and interfacial forces

are negligible. Under such circumstances, the acceleration in the momentum equation diminishes,

leading to

∇ · σ + ρb = 0; or

∂σij
∂xj

+ ρbj = 0.

(2.44)

This equation is termed the equilibrium equation. Situations to which the equilibrium equation

applies are commonly referred to as static problems.

Conservation of Energy: The principle of conservation of energy, also known as the energy bal-

ance principle, asserts that the rate of change of total energy is equal to the work done by the body

forces and surface reactions, in addition to the heat energy supplied to the body by the heat flux

and other sources of heat.

The rate of change of total energy in the body is expressed as

Ptot = Pint + Pkin, Pint =

Z

D
Dt

Ωρ ˙wintdΩ, Pkin =

Z

1
2

Ω

D
Dt

ρv · vdΩ,

(2.45)

here, Pint represents the rate of change of internal energy, and Pkin represents the rate of change of

kinetic energy. The rate of work performed by the body forces in the domain and the reactions on

the surface is given by

Z

Z

Pext =

v · ρbdΩ +

v · tdΓ =

Ω

Γ

Z

Ω

Z

viρbidΩ +

vitidΓ.

Γ

(2.46)

The power supplied by heat sources s and the heat flux q is

Z

Z

Z

Pheat =

ρsdΩ −

n · qdΓ =

ρsdΩ −

Ω

Ω

Z

Γ

niqidΓ,

(2.47)

In the heat flux term, the negative sign signifies that positive heat flow is directed outward from

22

Eulerian description
Mass conservation
Linear momentum conservation

Energy conservation
Lagrangian description
Mass conservation
Linear momentum conservation

Energy conservation

ρ Dwint

Tensor form
Dρ
Dt + ρdiv(v) = 0
Dt = ∇ · σ + ρb
ρ Dv
Dt = D : σ − ∇ · q + ρs
Tensor form
ρ(X, t)J(X, t) = ρ0(X)
∂t = ∇0 · P + ρ0b
ρ0
ρ0 ˙wint = ρ0

∂wint(X,t)
∂t

∂v(X,t)

= ˙FT : P − ∇0 · ˆq + ρ0s

ρ0

Indicial Form
˙ρ + ρvi,i
Dt = ∂σij

∂xj

+ ρbi

ρ Dvi

Indicial Form
ρJ = ρ0
∂t = ∂Pij

∂Xj

∂vi(X,t)

+ ρ0bi

Table 2.1: Conservation equations.

the body. The conservation of energy is expressed as

Ptot = Pext + Pheat.

(2.48)

In other words, the first law of thermodynamics states that the rate of change of the total energy

in the body, comprising internal and kinetic energies, is equivalent to the rate of work performed

by external forces and the rate of work supplied by heat flux and energy sources.

The conservation equations are succinctly presented in the following Table 2.1, providing both

tensor and indicial forms.

23

CHAPTER 3: CONSTITUTIVE

RELATIONS

3.1: Linear elasticity and Hooke’s Law

Consider the stress strain curve σ = f (ϵ) of a linear elastic material subjected to uni-axial stress

loading conditions. For a given value of the strain ϵ, the strain energy density per unit volume ψ =

ˆψ(ϵ), is defined as the area under the curve (see Fig.3.1), which can be expressed mathematically

as

Figure 3.1: Stress-strain curve for a linear elastic material subject to uni-axial stress σ.

According to this definition,

σ =

∂ ˆψ
∂ϵ

= Eϵ.

In general, for possibly nonlinear elastic materials

σij = σij(ϵ) =

∂ ˆψ
∂ϵij

.

(3.1)

(3.2)

It helps define Hooke’s Law, the most general linear relation among all the components of the stress

and strain tensor

σij = Cijklϵkl.

(3.3)

24

E1σεIn this expression, Cijkl are the components of the fourth-order stiffness tensor of material prop-

erties of elastic moduli. The fourth-order tensor has 81 and 16 components for three-dimensional

and two-dimensional problems, respectively. The strain energy density in this case is a quadratic

function of the strain

ˆψ(ϵ) =

1
2

Cijklϵijϵkl.

The stiffness tensor can be written in two different orthonormal basis as

C = Cijklei ⊗ ej ⊗ ek ⊗ el.

(3.4)

(3.5)

The stiffness tensor has the following minor symmetries which result from the symmetry of the

stress and strain tensors

σij = σji ⇒ Cijkl = Cjikl.

(3.6)

This symmetric characterization reduces the number of material constants from 81 = 3×3×3×3 →

54 = 6 × 3 × 3. In a similar fashion we can make use of the symmetry of the strain tensor

ϵij = ϵji ⇒ Cijkl =⇒ Cijlk.

(3.7)

This further reduces the number of material constants to 36 = 6 × 6.

3.2: Isotropic linear elastic materials

The constitutive law for an isotropic, linear elastic and homogeneous material is given by the fol-

lowing equation

Cijkl = λδijδkl + µ (δikδjl + δilδjk) ,

(3.8)

where δij is the Kronecker-delta, and λ, µ are the Lame constants. Replacing in σij = Cijjklϵkl,

gives

σij = λδijϵkk + µ (ϵij + ϵji) .

(3.9)

25

The compliance form of Hooke’s law can be written in indicial notation as

ϵij =

1
E

[(1 + ν) σij − νσkkδkk] .

(3.10)

It can be written in matrix form

2

6
6
6
6
6
6
6
6
6
6
6
6
6
6
4

=

1
E

3

7
7
7
7
7
7
7
7
7
7
7
7
7
7
5

2

6
6
6
6
6
6
6
6
6
6
6
6
6
6
4

ϵ11

ϵ22

ϵ33

2ϵ23

2ϵ13

2ϵ12

1 −ν

1

−ν

−ν

1

0

0

0

2(1 + ν)

0

0

0

0

2(1 + ν)

0

0

0

0

0

symm

2(1 + ν)

3

7
7
7
7
7
7
7
7
7
7
7
7
7
7
5

3

2

7
7
7
7
7
7
7
7
7
7
7
7
7
7
5

σ11

σ22

σ33

σ23

σ13

6
6
6
6
6
6
6
6
6
6
6
6
6
6
4

σ12

.

(3.11)

Invert and compare with

2

6
6
6
6
6
6
6
6
6
6
6
6
6
6
4

=

1
E

3

7
7
7
7
7
7
7
7
7
7
7
7
7
7
5

2

6
6
6
6
6
6
6
6
6
6
6
6
6
6
4

σ11

σ22

σ33

σ23

σ13

σ12

λ + 2µ

λ

λ + 2ν

λ

λ

0 0 0

0 0 0

λ + 2µ 0 0 0

ν 0 0

ν 0

ν

symm

3

2

7
7
7
7
7
7
7
7
7
7
7
7
7
7
5

6
6
6
6
6
6
6
6
6
6
6
6
6
6
4

3

7
7
7
7
7
7
7
7
7
7
7
7
7
7
5

ϵ11

ϵ22

ϵ33

2ϵ23

2ϵ13

2ϵ12

,

(3.12)

and concluded that

λ =

Eν
(1 + ν)(1 − 2ν)

; µ = G =

E
2(1 + ν)

.

(3.13)

3.3: 2D Formulations for Solid Mechanics

It is assumed that 2D means the xy plane, and that z is the out-of-plane direciton. Plane-stress and

Plane-strain are two terms commonly come across in finite element analysis. These are essentially

two variations of 2D analysis, with difference in their assumptions.

26

3.3.1: Hooke’s Law for Plane Stress

For the simplification of plane stress, where the stresses in the z direction are considered to be

negligible. σ33 = σ23 = σ13 = 0, the stress-strain compliance relationship for an isotropic material

reduces to a simple 3 × 3 matrix, the compliance matrix is written as

3

7
7
7
7
5

2

6
6
6
6
4

ϵ11

ϵ22

2ϵ12

=

1
E

2

6
6
6
6
4

1 −ν

−ν

0

1

0

3

2

7
7
7
7
5

6
6
6
6
4

3

7
7
7
7
5

.

σ11

σ22

σ12

0

0

2(1 + ν)

(3.14)

The stiffness matrix for plane stress is found by inverting the plane stress compliance matrix, and

is given by

3

7
7
7
7
5

2

6
6
6
6
4

σ11

σ22

σ12

=

E
1 − ν2

2

6
6
6
6
4

1 ν

ν 1

3

2

7
7
7
7
5

6
6
6
6
4

0

0

0 0 1−ν
2

3

7
7
7
7
5

.

ϵ11

ϵ22

2ϵ12

(3.15)

3.3.2: Hooke’s Law for Plane Strain

In the case of plane strain, where the strains in the z direction are considered to be negligible,

ϵ33 = ϵ23 = ϵ13 = 0, the stiffness matrix reduces to a simple 3 matrix

2

6
6
6
6
4

=

1
E

3

7
7
7
7
5

2

6
6
6
6
4

ϵ11

ϵ22

2ϵ11

1 − ν2

−ν(1 + ν)

−ν(1 + ν)

1 − ν2

0

0

0

0

2(1 + ν)

3

2

7
7
7
7
5

6
6
6
6
4

3

7
7
7
7
5

.

σ11

σ22

σ12

Inverting stiffness matrix gives compliance matrix, which is written as

=

E
(1 + ν)(1 − 2ν)

2

6
6
6
6
4

3

7
7
7
7
5

2

6
6
6
6
4

σ11

σ22

σ12

ν

0

0

0

3

2

7
7
7
7
5

6
6
6
6
4

3

7
7
7
7
5

.

ϵ11

ϵ22

2ϵ12

1 − ν

0

1−2ν
2

1 − ν

ν

(3.16)

(3.17)

27

3.3.3: Relations between Plane-stress and Plane-strain

Differences

• In plane-stress analysis, there can be strain in the thickness of the element. If the element

is stretched, it wll get thinner, and if compressed, it will get thicker. Plane-stress is gener-

ally suitable for the analysis of an element with relatively limited depth in the out-of-plane

direction. This analysis is typically offered only in structural or generic FE analysis software.

• In plane-strain analysis, the deformation in the out-of-plane (thickness) direction is fully fixed

so no deformation could take place. Plane-strain is generally suitable for cross-sectional

analysis of elements of almost infinite depth in the out-of-plane direction.

28

CHAPTER 4: MACHINE-LEARNED

CONSTITUTIVE MODE FOR SINGE

4.1: Introduction

This chapter introduces a data-driven modeling framework that extends the environmental damage

modeling presented in Chapter 2. The objective is to describe the coupled effects of mechanical

and environmental damage on elastomeric matrices using a hybrid micro-mechanical and machine

learning approach.

The elastomer matrix is modeled as a cooperative multi-agent system, where each agent is

implemented as a conditionally activated neural network. These agents are rigorously constrained

by the principles of thermodynamics and continuum mechanics. Compared to traditional physics-

based methods, the integration of neural network agents within the micro-mechanical framework

enables a new generation of hybrid models. These models retain the interpretability and simplicity

of classical micro-mechanics while leveraging the adaptability and computational efficiency of

machine learning.

It is important to note that the predictive accuracy of the proposed hybrid model depends heavily

on the quality and diversity of the training data.

The formulation of this model adheres to several key design constraints:

• The model is defined in terms of strain energy to ensure physical consistency;

• A micro-sphere approach is used to reduce the dimensionality from 3D to 1D representations;

• A network decomposition strategy is adopted to isolate different inelastic mechanisms;

• Each one-dimensional subnetwork is represented by an independent learning agent.

These elements together constitute a cooperative multi-agent framework for simulating the com-

plex response of elastomeric materials under coupled damage conditions.

29

4.2: Cooperative Multi-Agent System

Continuum Mechanics

In the finite deformation regime, strain energy-based constitutive mod-

els offer significant advantages over classical stress-strain relations. These models inherently sat-

isfy several key physical and mathematical conditions—such as normalization, growth behavior,

isotropy, objectivity, and polyconvexity—ensuring uniqueness and stability of the solution. More-

over, strain energy-based formulations can be directly verified against the second law of thermo-

dynamics, particularly through the Clausius–Duhem inequality. They are commonly constructed

with respect to various work-conjugate stress-strain pairs, including:

Two-point tensors:

(F : P) → F: deformation gradient, P: first Piola-Kirchhoff stress;

Material tensors:

(E : S) → E: Green-Lagrange strain, S: second Piola-Kirchhoff stress;

Spatial tensors:

(L : τ ) → L: Hencky strain, τ : Kirchhoff stress.

To address the challenge of incomplete data in modeling the finite deformation behavior of

elastomers, we adopt a physics-informed order-reduction strategy. This includes the microsphere

concept, network decomposition, continuum mechanics constraints, and polymer physics insights.

The result is a micro-mechanical model where deep-learned agents, embedded as neural networks,

represent 1D constitutive responses. These agents collectively capture the macroscopic behavior

of the matrix based on learned local rules.

Building on our prior work [44], we reuse the knowledge integration framework developed

for mechanically induced deformations, as illustrated in Fig. 4.1. The major innovation here is

the replacement of the original network structure with a new class of machine learning models—

conditional neural networks (CondNNs)—to define the learning agents (L-agents).

30

Figure 4.1: Schematic of the proposed model from order-reduction to model fusion.

3D-to-1D Transition via Microsphere Concept

In an undeformed (virgin) matrix, polymer

chains are assumed to be uniformly distributed in all spatial directions, reflecting initial isotropy.

31

Polymer MatrixSub-Network 2Sub-Network 1Sub-Network 2Sub-Network 1MicrosphereSub-Element NsSub-Element 1ElementSub-Element NsSub-Element 1Sub-Element 1Sub-Element NsSub-Element jIW11IW2111L-Agent𝝍j𝒅𝒊 𝝀j𝒅𝒊 𝝀j𝒎𝒂𝒙𝒅𝒊 Polymer ChainsL-AgentAgent TypeL-AgentAgent Team1NdMulti-Agent System213Order-ReductionModel FusionL-AgentNsL-Agent1L-AgentNdL-Agent111L-AgentNdL-Agent1NsNsSub-ElementsCooperative Learning   Multi-Agent SystemKnowledgeModelingUsing the microsphere approach, we represent the 3D matrix as a superposition of many 1D ele-

ments distributed over a unit sphere. The total strain energy is computed by integrating the direc-

tional contributions over the surface of the sphere:

Ψm =

Z

S

1
4π

m dSd ∼
Ψd
=

NdX

i=1

wiΨdi

m, with Ψdi

m = Bdi,

(4.1)

where Nd denotes the number of integration directions, wi are the associated weights, and Bdi

represents a team of L-agents modeling the energy in direction di. Each team comprises multiple

sub-elements:

Bdi =

NsX

j=1

Ai
j,

(4.2)

where Ns is the number of sub-elements per direction and Ai

j is the L-agent representing a sub-

element’s contribution. In the virgin state, all directions are identical: Bdi = Bdj . As deformation

progresses, directional differences emerge, enabling the model to naturally describe anisotropy,

damage initiation, and failure propagation.

Figure 4.2: Microsphere-based 3D-to-1D transition model.

Network Decomposition To capture the complex behavior of each 1D element, we apply net-

work decomposition, representing the total element response as the superposition of Ns simpler

sub-element responses:

32

ee3e2e1diSub-element 1Sub-element jPolymer MatrixMicrosphereElementΨdi

m =

NsX

j=1

Ψdi

j , where Ai

j := Ψdi
j .

Combining with Eq. 4.1, the total matrix strain energy becomes:

NdX

NsX

Ψm

∼
=

wiΨdi
j

:=

NdX

NsX

wiAi
j.

i=1

j=1

i=1

j=1

(4.3)

(4.4)

Each sub-element behaves as a scalar-to-scalar function and is implemented via a simplified

feed-forward neural network. The resulting first Piola–Kirchhoff stress tensor is derived using the

chain rule with respect to the deformation gradient:

P =

∂Ψm
∂F

− pF−T =

NdX

NsX

i=1

j=1

wi

∂Ai
j
∂F

− pF−T ,

(4.5)

where p is the Lagrange multiplier enforcing incompressibility.

4.3: Mechanical and Environmental Damages in Polymers

In general, elastomer degradation is caused by two types of stresses: mechanical, such as deforma-

tion, vibration, and relaxation, and environmental, such as thermal aging, hydrolysis, and photo-

oxidation.

Deformation-induced damage: One good example of mechanical damages is the Mullins ef-

fect, where the matrix softens after the first stretch [34, 109, 120]. This phenomenon happens in

filled rubber polymers, not all polymers. There are multiple sub-structural changes in the matrix

that are formed due to mechanical damages, such as chain breakage in the filler interface, chain

disentanglement, molecules slipping, and rupture in the cluster of fillers [13, 53], which are sum-

marized in Fig. 4.3.

Environmental damage occurs mostly due to changes in the chemical micro-structure of the

matrix, which occurs due to competition or collaboration of multiple sub-mechanisms. For exam-

ple, thermal aging is widely considered as the competition between cross-link formation mechanism

and chain scission mechanisms, where the former makes the matrix stiffer by inducing new bonds.

33

Figure 4.3: Schematic of polymer matrix change alongside with deformation-induced damage,
which include breakages of chain crosslinks, chain scission, bond rupture, desorption of chains,
slippage of chains, absorption of chains, and filler rupture.

In contrast, the latter makes it softer by debonding polymer chains from the system (Fig. 4.4). The

rate of polymer chain scission and cross-links formation defines that the material should become

brittle or ductile, although in both cases, the matrix experiences damage [86].

34

ABCBond ruptureMolecules slippingB'A'C'FillerruptureFigure 4.4: Schematic of polymer matrix change alongside with environmental-induced damage.
Environmental-induced damages include chain scission, cross-links formation and reduction.

Thermo-oxidative kinetics: During thermo-oxidative aging, the rate of chemical oxidation

can be characterized as

− d[P ]
dt

= k[P ]q,

(4.6)

[P ] is the chemical compound concentration of P , k is the reaction rate coefficient, and q is the

reaction order. The chemical reaction governing a degradation process is commonly described

using first-order kinetic equations (q = 1). Furthermore, k is just a function of temperature in the

presence of homogenous conditions, low stretches during aging, and no diffusion-limited oxidation.

The major cause of alterations in the polymer matrix is chemical interactions between the polymer

backbone and oxygen during thermal-induced aging [43, 69, 88].

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

scribes a different energy-dissipating damage mechanism. Using the concept of micro-sphere, each

35

Aged matrixChain scissionCross-links formationChain scissionCross-links reductionVirgin matrixAged matrixnetwork is considered as a 3D composition of infinite 1D subnetworks that are distributed in all

spatial directions. Subnetworks can only be subjected to uniaxial deformation and thus will expe-

rience different deformations based on their directions. To develop a model for the subnetwork

in direction, only a simplified form of entropic energy is needed with respect to uniaxial deforma-

tion. Integrating a subnetwork in all directions, the consequent network is a representation of that

concept in a 3D configuration. The similar concept will be used for training of the neural network.

As a result, network decomposition helps us describing the complex behaviour of a 1D element

as the super-position of the simple behaviour of Ns sub-elements, where Ψm

P

di =

Ns
j=1 Ψj

di. Here,

each sub-element is responsible for one or no inelastic behaviour. Assigning L-agent to each sub-

element, Ai
P

Ns
j=1

Ai

j := Ψj

di, an element behaviour is described by a team of cooperative L-agents Bdi =

j. Consequently, by substituting Eq. 4.1, we can extract the energy of the matrix directly

with regard to sub-elements and the L-agents as

Ψm =

Z

S

1
4π

ddS d ∼
=

Ψm

NdX

NsX

wiψj

di :=

NdX

NsX

wiAi
j.

i=1

j=1

i=1

j=1

(4.7)

Considering each sub-element is represented by an L-agent, the super-simplified scalar-to-scalar

mapping behavior of a sub-element can be modeled by a simplified feed-forward neural network.

Summarizing the implemented constraints, the first Piola—Kirchhoff stress tensor P can be derived

on the basis of Eqn. 4.6 and Eqn.4.7, as

P =

∂Ψm
∂F

− pF−T :=

NdX

NsX

i=1

j=1

wi

∂Ai
j
∂F

− pF−T ,

(4.8)

where p signifies the Lagrange multiplier to ensure incompressibility of the material.

4.3.1: Conditional Neural Network (CondNN) L-Agent

There is a new family of hybrid machine learning algorithms known as conditional neural networks

(CondNNs) for problems in which the outputs are not only dependent on past occurrences, e.g., de-

formation effects on the matrix, but also on external actions, e.g., temperature and time of aging

36

effects on the polymer matrix. This new algorithm has been derived by the combination of general

neural network and decision tree concepts (see [63,110]). A routed behavior is a feature of decision

trees; the data is sent to one or more children based on some learned routing function. This condi-

tional computation means that we can infuse the knowledge in the model. Meanwhile, throughout

many tasks, NNs achieve industry-leading precision for learning, but decision trees have the ability

to infuse the knowledge into NNs. In other words, CondNNs are decision trees with the exception

that instead of moving the data as is, each node applies a non-linear transformation to it. We can

also forward the data to one or more of its children using routers.

In order to design an L-agent which can model mechanical and environmental damage in the

polymer matrix, we use CondNNs as a new hybrid framework. The feed-forward L-agents consist

of NNs followed by a fully connected layers, which the first branch represents the mechanical

damage in the polymer matrix. The resulting features are then combined multiplicatively with the

second branch of NN followed by fully connected layer that represents environmental damage in

the matrix (see Fig. 4.6).

- Neural Network Architecture The architecture of the neural network (Fig. 4.6) is determined

through the summary of the hyperparameters: (i) number of hidden layers, nl (network depth), (ii)

number of neurons per hidden layer, nn (network width), (iii) activation function.

Activation function : The activation function of a node in artificial neural networks determines

the output of that node given input or a set of inputs. This is similar to the linear perceptron’s

behavior in neural networks. However, only nonlinear activation functions allow such networks,

using only a small number of nodes, to measure non-trivial problems, and such activation functions

are called non-linearities. Choosing the activation function has a significant impact on the neural

network’s capacity and efficiency, and various activation functions can be used in different sections

of the model. Based on the type of neural network architecture, the activation function used in

hidden layers is usually chosen. If we are unsure which activation function to use for our network,

we try and compare the results to find the best ones. We have done the same procedure in this study.

Note that the choice of activation function in the hidden layer will control how well the network

37

model learns the training dataset. Also, the effect of choosing different activation functions has not

a significant effect on training time [91].

Figure 4.5: a) ReLU(Rectified Linear Unit), ELU(Exponential Linear Unit) and tanh activation
functions, b) first derivatives of activation functions, and c) second derivatives of activation
functions [61].

For training of the described problem, it is crucial to use activation functions tailored to the

underlying physical output quantities. Using physics-based rationale, an activation function with

a positive output parameter for tension and a negative output parameter for compression and non-

zero derivative (which is important in physics-based neural networks) is an obvious choice. For

all hidden layers, the hyperbolic tangent function is used, which is a preferable activation function

due to its smoothness and non-zero derivative.

Note that Rectified Linear Unit activation function or ReLU for short is a piecewise linear

function that will output the input directly if it is positive, otherwise, it will output zero. Also, The

Exponential Linear Unit (ELU) is an activation function for neural networks. In contrast to ReLUs,

ELUs have negative values which allows them to push mean unit activation closer to zero.

The L-agent response is computed using a feed -forward algorithm for a given set of hyper-

parameters (nl, nn). Each L-agent can be represented by a CondNNs as

Ai

j := Ddi(Ei)ψdi

j (Mj

i, Sj

i)

Ddi = CN Ne(We, Ei),

ψdi
j = CN N m(Wj

m, Mj

i, Sj

i), (4.9)

where ψdi
j

represents the energy of deformation-induced damage part and Ddi is related to envi-

38

ronmental damage of one sub-element. Similarly, two weight matrices Wm

(cid:2)

(cid:3)

(cid:2)

j =

(cid:3)

Wm

j
1...Wm

j
nl+1

and We =

We1...Wenl+1

are related to weight matrices of mechanical and environmental damage

CondNNs, respectively. Here, ψi

j(Mj

i, Sj

i) is trained on the basis of a non-kinematic input sets Mj

i

and internal parameters Sj

i for the mechanical damage CondNN. Selection of internal parameters

depend on (full or recent) material memory. Also, the input vector Ei should represent the setting

of the problem related to environmental damage such as temperature and time of degradation.

Normalization, conditions of growth, isotropy, objectivity, and poly-convexity are already sat-

isfied in the proposed equation (Appendix).

Figure 4.6: Schematic of a neural network with physics-informed engineered features to learn the
effect of environmental aging on the mechanical Behaviour. This is referred to as CNN in this
study.

- Model parameter identification The parameters of the CondNN model are identified using

a Gradient descent algorithm. A data collection consisting of n data points is the basis for model

recognition. Loss function was defined with respect to the difference between target values, AKA

experimental data, and approximated value, AKA CondNN output. Here, we define the loss func-

tion L Mean Squared Error (MSE) for a total of ntot data points as

L(W1

m, ...., Wj

m, We) =

ntotX

0

@

2
4g1

NdX

NsX

n=1

i=1

j=1

1
2

∂Ai
j
∂F

wi

− pF−T

1
A g1

3

2

5

,

− P 11
n

(4.10)

where P 11
n

:= g1Pg1 is the first component of the experimental macro-scale stress tensor Pn in

loading direction g1 for point n.

39

IIIII111   𝝀max𝝀𝝍 RouterW1m1nl = 1nn = 2nl = 1nn = 4I1𝑡 𝜃W1m2We1We24.4: Validation: Environmental + Mechanical Damages in Rub-

ber

Here, the proposed hybrid model is developed and validated for two different loading scenarios

1. Thermal-induced aging + Mechanical deformation (Mullins effect)

2. Hydrolysis (desalinated water) + Mechanical deformation (Mullins effect)

In both cases, the model was benchmarked against experimental data to predict the inelastic behav-

ior of cross-linked elastomers for different states of deformation at different stages of aging. We

also validated the model predictions for different coupling scenarios for mechanical and environ-

mental damages, namely

- Intermittent test where the sample has been aged first with no stress, and then the mechanical

behaviour is characterized

-Relaxation test where the samples were loaded and aged at constant deformation and the drop

in the stress has been recorded.

In both cases, the predictions were validated against four main variables the effect of (i) defor-

mation, (ii) deformation history, (iii) aging time (t), and (iv) aging temperature (θ).

4.4.1: Network Architecture

Wwe used identical engines; a relatively simple engine built by Nd = 21 teams, where each team

had Ns = 2 agents [38]. Note that the number of teams and their related agents can be chosen

based on the trade-off between accuracy and computational cost. In other words, we consider 21

teams because it is a small number for estimation of the integral with a summation. As a result, it

guarantees that the proposed model error in prediction is excellent even in a small number of teams.

To capture all deformation states, we used only two agents, each of which is representative of the

first and second invariants of the Green-Cauchy deformation tensor. For the CondNNs structure of

L-agents, we considered one input layer, one hidden layer with four neurons and three activation

functions soft plus ψ(•) = ln(1 + e•), sinusoid ψ(•) = sin(•) and hyperbolic tangent ψ(•) =

40

tanh(•).

The internal parameters of L-agents were built via λj−max parameters to capture the deformation

of the rubbers with full memory. In order to allow teams to predict various deformation states, the

first and second deformation invariants were supplied to each team [30, 74]. The condition was

satisfied by providing input sets into the first and second L-agents as

Mdi

1 = [λdi],

while

Sdi
1 = [λdi

max], Edi = [t, θ], Mdi

2 = [νdi],

Sdi
2 = [νdi

max], Edi = [t, θ], (4.11)

p

diCdi,

λdi =

p

νdi =

diC−1di,

C = FT F,

(4.12)

where [di]i=1...Nd

is integration directions in micro-sphere, λdi and νdi are related to I1 and I2, as the

first and second invariants of C, respectively.

In summary, the rubber matrix was represented by a cooperative game of 21 teams of 2 agents

through Ai

j, i ∈ {1, 21} , j ∈ {1, 2}. After agent fusion, the final cost function is given by

L(W1

m, W2

m, We) =

X

21X

2X

[g1(

wi

1
2

∂Ai
j
di

∂λj

di

∂λj
∂F

− pF−T )g1

− P 11

n ]2,

(4.13)

j=1
subjected to weights related to λmax and νmax ⩽ 0, and weights related to λ and ν ⩾ 0 to satisfy

i=1

n=1

thermodynamic consistency and poly-convexity, respectively. Accordingly, the energy of each sub-

element can be written with respect to the deformation gradient F as follows

21X

i=1

wi

∂Ai
1
∂λdi

∂λdi
∂F

=

21X

i=1

wi

∂Ai
1
∂λdi

1
λdi

F (di ⊗ di),

(4.14)

21X

i=1

wi

∂Ai
2
∂νdi

∂νdi
∂F

= −

21X

i=1

wi

∂Ai
2
∂νdi

1
νdi

F−1F−T F−1 (di ⊗ di).

(4.15)

4.4.2: Step 1: CondNN Training in the absence of aging

To investigate the performance of the proposed model in the material in the absence of aging, we

41

benchmarked the inelastic features in the behavior of natural rubber, namely Mullins effect and per-

manent set. Figure 4.7 shows stress-stretch curves for this cross-linked polymer with experimental

data of [80]. We used one set of bi-axial loading-unloading until λ = 2.7 for training and predicting

inelastic effects in different states of deformation, e.g., uniaxial and pure shear at increasing stretch

amplitudes which constitutes deformation histories.

Figure 4.7: Model training and prediction of bi-axial, pure shear, and uni-axial(Urayama’s
dataset [80]). Dash lines stand for fitting and solid lines stand for prediction. The biaxial dataset
which have been showed with the blue line used for training and of the model. The black lines in
bi-axial, shear, and uni-axial datasets show the performance of the model in prediction after
training.

4.4.3: Case Study : Thermo-Oxidation + Mechanical Damage

Polyurethane A bone form punch was used to punch dumbbell samples in compliance with

ASTM D412 standards. The reason that we are using dog-bone samples is that in this shape, when

we put the extensometer in a specific location of that sample, we are sure that it gives us one com-

ponent of the deformation matrix, for example, x-x direction. But in complex shapes, several com-

ponents of stress tensor participate; so, it is not possible to see the effect of each component with

an external extensometer. Three different temperatures, 60◦C, 80◦C, and 95◦C, in zero humidity

level, were used to age the samples. The performance of the invented model was determined using

the results of the uni-axial tensile failure tests. Fig. 4.8 depicts the evaluation’s related findings.

42

StretchNominal Stress [Mpa]5432106711.41.82.22.6311.522.53.53Experimental data - UniaxialPredictionNominal Stress [Mpa]Stretch543210543210611.41.82.22.6Experimental data - ShearPredictionNominal Stress [Mpa]StretchExperimental data - BiaxialTrainingPredictionFigure 4.8: Training and model predictions for polyurethane against intermittent test; (a) at 60◦C,
(b) at 80◦C, (c) at 95◦C. Dash lines stand for fitting and solid lines stand for prediction. The black
(unaged) and red points in the Fig. a), the black (unaged) and blue points in Fig. b), and the black
(unaged) and green points in the Fig. c) have been used for training. The solid lines show the
performance of model in prediction after the training.

Silicon Adhesive To assess the model’s capabilities, we compared its predictions to our exper-

imental results, which were specifically developed to demonstrate the effect of chemical aging

on the constitutive response of silicone adhesive. In this case, samples were placed in ovens at

temperatures 60◦C, 80◦C, and 100◦C with a relative humidity of zero (i.e. RH= 0%). All of the

specimens were aged under constant pressure, and after a set amount of time, they were removed

from the containers and dried with tissue paper. The predictions of the proposed model against the

experimental data for different types of aging and different amount of aging times, temperatures,

and deformations are plotted in figure 4.9.

43

30 days10 days1 daysUnagedPredictionTraining11.21.41.61.82StretchTemp. = 600246Stress [MPa]30 days10 days1 daysUnagedPredictionTraining11.21.41.61.82Stretch0246Stress [MPa]30 days10 days1 daysUnagedPredictionTrainingTemp. = 8011.21.41.61.82Stretch0246Stress [MPa]Temp. = 95a)b)c)Figure 4.9: Training and model predictions for silicone adhesive against intermittent test; (a) at
60◦C, (b) at 80◦C, (c) at 100◦C. Dash lines stand for fitting and solid lines stand for prediction.
The black(unaged) and green points in the figure a), the black(unaged) points in figure b), and the
black(unaged) and green points in the figure c) have been used for training. The solid lines show
the performance of model in prediction after the training.

Note that the model is able to predict the behavior of elastomers even in the stretch greater than

4. Thermo-oxidative aging makes the samples brittle. As a result, the failure point for samples

happens on the less stretch. So, because the focus of the proposed model is on the behavior of aged

materials we showed the performance of the model till stretch less than 4. In addition, the kernel

of model is the same even for higher stretches. So, the model works for higher stretch as well.

4.5: Conclusion

The aim of this chapter was to propose a novel physics-informed hybrid framework to capture the

relation between elastomeric network mechanics and environmental damage. The proposed model

is used to describe the effect of single-mechanism aging, such as thermal-induced or hydrolytic

44

21 days12 days72 hoursUnaged11.11.21.31.41.5012345678StretchStress [MPa]PredictionTrainingTemp. = 6011.11.21.31.41.5Stretch012345678Stress [MPa]Temp. = 808 days4 days36 hoursUnagedPredictionTraining11.11.21.31.41.5Stretch012345678Stress [MPa]Temp. = 1006 days72 hours48 hoursUnagedPredictionTraininga)b)c)aging, on the behavior of the material in this hybrid system. Huge chain scission, cross-link re-

duction, chain formations, and changes in polymer morphology are all examples of environmental

single-mechanism damages that alter the polymer matrix over time. In other words, this model is

a predictive model and not a descriptive model. As a result, it can provide prediction in a frac-

tion of time and cost to the end-users for the reliable design of components. Although, The model

cannot be used to understand or explore kinetics. By extracting micro-structural activity from

macroscopic experimental data, the data-driven approach aims to overcome the shortcomings of

both phenomenological and micro-mechanical models.

The model arises from polymer physics and an order reduction strategy ending to the con-

strained L-agents training. The polymer matrix was described by a cooperative multi-agents system,

in which each agent is represented by a simple deep-learned CondNN that is super-constrained by

laws derived from physics, thermodynamics, and continuum mechanics. While satisfying the con-

tinuum mechanics and thermodynamics rules, the proposed hybrid model is quite simple for the

most rubbery media in the extreme environment.

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 single mechanism agings. The efficiency of the model was found sat-

isfactory, and in some cases, excellent when compared with the experimental data. ”In general,

the error of our model in training and prediction is less than 10%.” Accuracy and simplicity of the

model make it a proper choice for commercial and industrial application because we do not need

to know the exact behavior and interaction of micro-structures; however, in the future, the model

can be extended to consider viscoelasticity and non-isotropic formation for better precision due to

platform of the model. Besides, note that the proposed model focus is on the mechanical behavior

of the material during degradation, not lifetime prediction. Also, the performance of the proposed

model is completely dependent on the dataset that we choose for training. It is axiomatic that if we

increase the available dataset for training, the predictionability of the model is more accurate.

45

CHAPTER 5: FINITE ELEMENT

ANALYSIS

5.1: Introduction

The study of polymer mechanics focuses on the relationship between the mechanical response of

polymers and externally applied forces. This field plays a crucial role in the design, optimization,

and prediction of the mechanical performance of polymer-based components used in a broad range

of engineering applications.

Historically, polymers were considered inexpensive materials with limited structural applica-

tions. However, advances in theoretical modeling and experimental techniques have repositioned

polymers as high-performance, reliable, and tunable materials. These developments, coupled with

a better understanding of polymer behavior under different loading and environmental conditions,

have expanded their role in critical engineering systems.

Traditional mechanical analysis techniques such as closed-form solutions were largely limited

to simple geometries and loading scenarios. With the advent of modern computational tools, the

Finite Element Method (FEM) has become the predominant technique for simulating and analyzing

the mechanical behavior of polymers. The widespread availability of commercial FEM software

packages has made it possible to solve complex geometries with nonlinear material behavior. Nev-

ertheless, the use of such tools demands careful modeling assumptions, validation, and physical

insight to avoid erroneous results.

This chapter provides a comprehensive overview of finite element modeling strategies for poly-

meric materials, emphasizing nonlinear formulation, hyperelastic constitutive models, and consis-

tent tangent stiffness derivation for robust computational implementation.

46

Analysis Method
Linear Elasticity

FE with Basic Material
Model

FE with Advanced Mate-
rial Model

analytical

efficient;
insight

Strengths
Computationally
provides
into basic deformation
nonlinear
Captures
ge-
ometries;
straightforward
implementation
Captures nonlinear geometry
and material responses with
higher accuracy

Limitations
Inadequate for complex ge-
ometries, nonlinearities, or
time-dependent behavior
Poor predictive accuracy for
large strain or rate-sensitive
effects
Requires extensive calibra-
tion and high computational
cost

Table 5.1: Comparison of common analysis methods for polymer mechanics.

5.2: Nonlinear Equilibrium and Directional Derivatives

The governing equations in finite deformation mechanics typically take the form of a nonlinear

boundary value problem. To solve these equations numerically, iterative approaches such as the

Newton-Raphson method are used. Central to this approach is the concept of the directional deriva-

tive, which enables the linearization of nonlinear equations.

5.2.1: Concept of Directional Derivative

Consider a two-spring system with total potential energy defined as

Π(x) =

1
2

kx2

1 +

1
2

k(x2 − x1)2 − F x2,

(5.1)

where x = [x1 x2]T and F is the applied load. For a small displacement increment u = [u1 u2]T ,

the change in potential energy can be approximated by the directional derivative:

DΠ(x)[u] ≈ Π(x + u) − Π(x).

(5.2)

Using Taylor expansion, the directional derivative becomes:

DΠ(x)[u] =

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

d
dϵ

ϵ=0

Π(x + ϵu) = uT (Kx − F),

(5.3)

47

where the stiffness matrix is given by:

2

6
4

K =

3

7
5 .

2k −k

−k

k

Setting DΠ = 0 for all admissible u leads to the equilibrium condition:

Kx = F.

(5.4)

(5.5)

5.3: Newton-Raphson Linearization

To solve a general nonlinear system R(x) = 0, the Newton-Raphson method approximates the

residual at each iteration using a first-order Taylor expansion:

R(xk+1) ≈ R(xk) + K(xk)u = 0,

⇒ K(xk)u = −R(xk),

xk+1 = xk + u.

Here, K is the tangent stiffness matrix:

Kij =

∂Ri
∂xj

.

(5.6)

(5.7)

(5.8)

5.4: Tangent Stiffness Matrix in Finite Element Analysis

For a typical finite element, the tangent stiffness matrix is assembled from element-level contribu-

tions:

K(e) = DT(e)(x)[u] =

2

6
4

Kaa Kab

Kba Kbb

3

7
5 ,

where T(e) is the internal force vector, and u is the displacement increment.

(5.9)

48

5.5: Hyperelastic Constitutive Models

In hyperelastic materials, stress is derived from a strain energy density function Ψ(F) that depends

on the deformation gradient F:

P =

∂Ψ
∂F

.

(5.10)

For incompressible materials, the constraint detF = 1 introduces a Lagrange multiplier p, lead-

ing to:

P = −pF−T +

∂Ψ
∂F

.

The second Piola-Kirchhoff stress S then becomes:

S = −pC−1 + 2

∂Ψ
∂C

.

(5.11)

(5.12)

5.6: Elasticity Tensors

To linearize constitutive models in Newton-based solvers, one defines the elasticity tensor C:

dS = C : dE, C = 4

∂2Ψ
∂C∂C

.

5.7: Isotropic Hyperelasticity

For isotropic materials, the strain energy depends only on the invariants of C:

Ψ = Ψ(IC, IIC, IIIC),

S = 2ΨII + 4ΨIIC + 2J 2ΨIIIC−1,

where ΨI = ∂Ψ/∂IC, and so on.

(5.13)

(5.14)

(5.15)

49

5.8: Spatial Form of Stress Tensor

The Cauchy stress tensor is recovered via:

σ =

1
J

FSFT .

Alternatively, using spatial derivatives of Ψ(b), one may compute:

σ = 2J −1 ∂Ψ
∂b

b.

5.9: Incompressible Hyperelasticity

In incompressible cases, the strain energy formulation becomes:

(5.16)

(5.17)

Ψ = Ψ(F) − p(J − 1),

(5.18)

and the constitutive relation for S must satisfy the incompressibility condition. Using the consis-

tency condition on allowable strain rates, the general form becomes:

S = 2

∂Ψ
∂C

+ γJC−1,

where γ is related to hydrostatic pressure p by:

p = γ +

2
3

J −1 ∂Ψ
∂C

: C.

(5.19)

(5.20)

5.10: Finite Element Implementation Using UMAT

The Finite Element Method (FEM) enables the integration of custom material models using user-

defined subroutines. In Abaqus, the UMAT interface (User MATerial) allows users to implement

constitutive laws written in Fortran, facilitating advanced material behavior beyond the built-in

50

options.

5.10.1: Role of UMAT in Nonlinear Finite Element Analysis

For finite deformation problems, accurate stress integration and consistent tangent stiffness are

essential for the convergence of Newton-Raphson iterations. The UMAT must return two key outputs

at each integration point:

• The second Piola-Kirchhoff stress tensor S.

• The algorithmic tangent modulus tensor C =

∂S
∂E

.

The input and output quantities are typically defined in the material (Lagrangian) configuration,

consistent with the total Lagrangian finite element formulation.

5.10.2: Inputs and Outputs of UMAT

The inputs passed to UMAT by the FEM solver include:

• Current deformation gradient F.

• Strain increment ∆E.

• State variables (internal variables).

• Temperature, time increment ∆t, and material constants.

The subroutine must compute:

• Updated stress tensor S.

• Updated internal state variables.

• Consistent tangent modulus C for global system assembly.

5.10.3: Constitutive Consistency in UMAT

The stress update in hyperelasticity is typically defined as:

S = −pC−1 + 2

∂Ψ
∂C

,

(5.21)

51

where C = FT F is the right Cauchy-Green deformation tensor, and Ψ is the strain energy density

function.

The consistent tangent modulus is computed as:

C = 4

∂2Ψ
∂C∂C

,

(5.22)

ensuring that the linearization of the residual vector R(x) matches the expected directional deriva-

tive:

DR(x)[u] = C : δE.

(5.23)

This guarantees quadratic convergence in Newton-Raphson iterations.

5.10.4: Benefits of Using UMAT

Feature
Custom Constitutive Modeling

Direct FEM Coupling

Research Flexibility

Robustness

Advantages
Implementation of rate-dependence, damage,
aging, or experimental formulations
Allows stress updates and tangent moduli to be
synchronized with global solvers
Enables prototyping of novel theoretical models
for polymers and composites
Improves convergence for complex nonlinear
problems with large deformations

Table 5.2: Advantages of implementing material models via UMAT.

5.10.5: Practical Considerations for UMAT Development

Developing a reliable UMAT routine requires:

• Rigorous derivation of the stress and tangent modulus consistent with the strain energy func-

tion.

• Verification against standard benchmarks: uniaxial, biaxial, and shear tests.

• Careful treatment of incompressibility, internal variables, and numerical stability.

52

Neglecting the consistent tangent modulus or implementing an inconsistent formulation can

lead to poor convergence behavior, especially for nearly incompressible or large-strain materials.

5.11: Summary

This chapter introduced the nonlinear finite element formulation for polymer mechanics, covering

the derivation of the directional derivative, Newton-Raphson linearization, tangent stiffness matri-

ces, and hyperelastic constitutive laws. The mathematical foundations established here are essential

for implementing robust user-defined material models (UMATs) and for the accurate simulation of

polymeric behavior under large deformations. The integration of user-defined constitutive models

via the UMAT subroutine extends the versatility of finite element solvers, allowing for more accu-

rate and physically motivated simulations, especially when dealing with complex phenomena such

as environmental degradation, aging, and damage evolution in polymers. When combined with

rigorous continuum theory and consistent linearization strategies, the UMAT interface serves as a

powerful bridge between advanced material model and commercial software.

53

CHAPTER 6: INERT THERMAL AGING

6.1: Abstract

This study presents a micro-mechanical model developed to investigate the effects of thermal

aging—excluding oxidative influences—on the constitutive behavior of cross-linked polymeric

systems, with a focus on elastomers. The model accounts for the simultaneous and competing

processes of chain scission and cross-linking that occur during aging. Although both mechanisms

are active, one typically dominates, leading to distinct changes in material properties. The cross-

link density may increase or decrease depending on the reaction kinetics, which are governed by

temperature and exposure duration.

The goal of this modeling effort is to develop a mechanistic understanding of how thermal ag-

ing influences the mechanical response of cross-linked polymers and to integrate this understand-

ing into a computational framework. By explicitly capturing the interplay between chain scission

and cross-linking, the model enables accurate prediction of the evolution of material properties

over time. To facilitate practical application, the constitutive model is linearized for implementa-

tion within a finite element method (FEM) framework. This includes derivation of the consistent

tangent modulus to ensure numerical stability and convergence during nonlinear simulations. The

resulting model has been implemented in a FEM environment, allowing for simulation of thermally

aged polymer behavior in complex geometries and loading conditions relevant to engineering ap-

plications.

By isolating inert thermal aging from oxidative effects, this work sheds light on fundamental

aging mechanisms specific to elastomers. The insights gained have practical implications for ma-

terial development and can inform improved design practices in applications where elastomeric

durability is critical.

54

6.2: Introduction

The space environment poses a complex set of challenges for material systems, characterized by

extreme vacuum, temperature fluctuations, radiation, plasma, and microgravity. Polymers exposed

to such conditions inevitably undergo aging and degradation, leading to progressive deterioration of

their mechanical, physical, and chemical properties. Over time, these changes compromise material

functionality, making durability and environmental resistance essential for polymeric components

in space applications.

This chapter introduces the harsh conditions encountered in space and their implications for

the long-term performance of polymeric materials, with a focus on adhesives and elastomers. Key

degradation mechanisms—such as thermal cycling, mechanical stress, electromagnetic radiation,

and ionizing particles—are reviewed in the context of their influence on the structural and func-

tional integrity of polymers. This background provides a foundation for material selection, perfor-

mance prediction, and reliability assessment in space systems.

Materials used in space must withstand unique stressors including thermal cycling, high vac-

uum, atomic oxygen, radiation, micrometeoroids, and debris impacts ( [18]). Elastomeric materials,

in particular, are susceptible to aging from various environmental influences such as heat, humidity,

radiation, and chemical exposure ( [114]; [23]). Among these, thermal aging—caused by prolonged

exposure to elevated temperatures—is recognized as one of the most significant contributors to

long-term degradation.

Extensive prior research has examined thermal aging in rubbers and cross-linked polymers.

The resulting changes in network structure are influenced by factors such as monomer composi-

tion, initial curing systems ( [22]), and aging conditions ( [93]). Moreover, the mechanical loading

state during aging has been shown to significantly affect material response ( [84]; [97]). Addi-

tional investigations have addressed aging-induced cross-linking, chain scission, and accelerated

degradation mechanisms ( [105]; [76]; [12]).

A series of works ( [45]; [46]; [66]; [65]; [64]) have modeled the thermo-oxidative behavior

55

of polymers by incorporating diffusion-reaction processes and mechanical coupling into thermo-

dynamically consistent frameworks. These models account for the evolution of internal variables,

such as the formation or degradation of cross-links and chain scission events, highlighting the com-

plex interplay between chemical and mechanical factors in aging ( [68]; [87]).

An alternative approach to understanding aging mechanisms is through monitoring molecu-

lar weight changes and distributions throughout degradation ( [73]). Under inert thermal condi-

tions, degradation can involve cross-reactions between macromolecules and low-molecular-weight

species, potentially enhancing or diminishing thermal stability depending on the reactivity of inter-

mediates.

Thermal aging induces microstructural transformations in polymers. For instance, high tem-

peratures can cause chain scission, leading to the coexistence of soft and hard phases within the

polymer matrix ( [1]; [104]). These changes manifest as modifications in chemical composition,

chain conformation, molecular weight, cross-link density, and the behavior of functional groups.

The following aspects are particularly affected:

• Chemical Composition: New species may form, altering the polymer’s chemistry.

• Chain Conformation: Spatial arrangements of polymer chains shift under heat.

• Molecular Weight and Distribution: Chain scission leads to shorter chains and broader

distributions.

• Cross-Linking and Branching: Additional cross-links or branched structures may form.

• Functional Groups: Their reactivity and distribution may change with temperature.

Cross-linked polymers generally exhibit improved thermal stability due to the need for multiple

bond breakages for degradation to occur ( [98]).

The goal of this research is to develop a predictive micro-mechanical model capable of es-

timating the inelastic behavior of cross-linked polymers subjected to thermal aging. The model

assumes homogeneous aging—an appropriate condition for ultra-thin samples aged in oxygen-free

environments—allowing a focused investigation into the aging phenomena specific to inert thermal

degradation.

56

6.3: Mechanical and Experimental Damages in Polymers

Cross-linked elastomers result from the extensive interconnection of polymer macromolecules

within a limited volume. Typically, the deterioration of elastomers can be linked to two primary

categories of forces: mechanical stresses like deformation, vibration, and relaxation, and deteriora-

tion stemming from diverse environmental factors, including elevated temperatures, immersion in

water, and exposure to ultraviolet (UV) radiation.

6.3.1: Deformation-Indeced Damage:

Material deformation frequently arises as an outcome of various mechanical stresses. These stresses

can induce several sub-structural alterations within the matrix, encompassing phenomena like chain

breakage, molecular displacement, and filler ruptures. These effects can be consolidated as follows:

Fig. 6.1.

Figure 6.1: Schematic of polymer matrix change alongside with deformation-induced damage,
which include breakages of chain crosslinks, chain scission, bond rupture, desorption of chains,
slippage of chains, absorption of chains, and filler rupture.

6.3.2: Environmental Damage:

This type of damage primarily stems from alterations in the chemical microstructure of the matrix.

It can manifest as a result of the interplay between multiple sub-mechanisms, such as chain scission

57

ABCBond ruptureMolecules slippingB'A'C'Fillerruptureand cross-linking, either in competition or collaboration. (see Fig.6.2). The relationship between

Figure 6.2: Schematic of chain scission and cross-linking processes.

tensile strength and crosslink density, before and after aging, is presented in Fig.6.3. The data reveal

a distinct peak in the tensile strength of NBR samples as crosslink density increases. This trend

can be attributed to two competing effects. Initially, increasing the crosslink density enhances the

number of junction points within the rubber network, which improves the uniform distribution of

applied stress and leads to higher tensile strength. However, beyond a certain threshold, further

crosslinking reduces the average molecular weight between crosslink points, resulting in stiffer

chain segments and the development of localized stress concentrations. This structural rigidity

ultimately diminishes the material’s ability to deform, causing a decline in tensile strength.

Figure 6.3: Tensile strength vs. Crosslink density of NBR vulcanizates before and after thermal
aging [54].

58

Cross-linkActive siteFillerCross-linkingScissionUnaged NetworkAged NetworkInert thermal aging20040060080010001200121314151617181920 Before Aging After AgingTensile Strength (MPa)Crosslink Density (mol.m-3)Figure 6.4: Mass loss curves as a function of temperature via TGA under inert and oxidative
environment [70].

6.4: Thermal Aging Mechanisms and Experimental Observations

6.4.1: Competing Mechanisms: Chain Scission vs. Crosslinking

A marginal increase in material hardness during aging presents an interpretive challenge: it is often

unclear whether this change stems from a subtle increase in crosslink density or from the simul-

taneous occurrence of both chain scission and crosslinking. In this study, we assume that both

mechanisms act concurrently, though one may dominate depending on the aging conditions.

As highlighted by Zaghdoudi et al. [119], chain scission reactions tend to reduce crosslink den-

sity, leading to material softening and degradation of elastic properties. In contrast, crosslinking

during aging typically results in increased stiffness and can lead to brittleness. Our model aims to

capture the combined effects of these opposing processes to accurately characterize material evolu-

tion under thermal aging. By identifying the dominant mechanism, we seek to improve predictions

of long-term mechanical performance in cross-linked polymer systems.

The theoretical basis for quantifying these effects can be traced to the work of P.J. Flory, who

used statistical mechanics to relate the elastic modulus of a polymer in the rubber phase to the num-

59

02004006008001000020406080100Weight (%)Temperature (oC) Nitrogen AirMass loss curves as a function of temperatureFigure 6.5: Evolution of the weight percentage loss as a function of the aging time (in days) for
CFRP samples [70].

ber of elastically effective chains. This relationship establishes a direct link between macroscopic

mechanical properties and underlying crosslink density.

6.4.2: Case Study: Thermal Vacuum Aging of Adhesives

Experimental data from Chen et al. [17] provide insight into how thermal aging under vacuum

conditions affects polymeric adhesives. In their study (see Fig.6.6), various commonly used curing

adhesives were exposed to a high vacuum (< 10−3 Pa) at a temperature of 90◦C for 48 hours. This

procedure simulated the combined effects of thermal and vacuum environments encountered in

aerospace applications.

Post-aging analyses revealed a marked reduction in tensile strength across all adhesives tested.

This degradation is attributed to the outward diffusion of small molecules, which led to microvoid

formation near the surface and weakened the adhesive’s resistance to tensile failure. The degree of

mass loss was found to correlate strongly with reductions in mechanical performance, indicating a

direct link between molecular migration and structural integrity.

Interestingly, the elastic modulus increased following aging, particularly in condensation-type

60

051015202530-8-6-4-20 Inert OxidativeAverage Mass Change (%)Aging Time (Days)silicone rubber and epoxy-based adhesives. This increase is attributed to additional crosslinking

reactions involving residual hydroxyl groups in the polymer network. Under thermal vacuum condi-

tions, these groups reacted with trace water molecules, promoting further intermolecular crosslink-

ing, increasing molecular weight, and resulting in a stiffer material.

Figure 6.6: Tensile stress and strain before and after thermal vacuum aging; (a) condensation-type
silicone rubber, (b) addition-type silicone rubber, (c) epoxy resin adhesive, (d) acrylic resin
adhesive [17].

The polymer utilized in the research conducted by Tcharkhtchi and colleagues [104] was a

polyurethane supplied by the RAIGI Company, located in Rouvray-Saint-Denis, France. To initiate

the aging process, the samples were placed within autoclaves that had been purged with argon

for a few minutes to establish an inert atmosphere. Subsequently, they were subjected to constant

temperatures of both 85◦C and 120◦C. To ensure temperature stability, the autoclaves’ temperatures

were closely monitored and recorded throughout the study.

Tensile tests were performed using an Instron machine equipped with 10 KN load cells, and the

61

0.00.51.01.52.02.53.00.00.10.20.30.40.5Tensile Stress (MPa)Tensile Strain (%) Before Aging After Aging0.00.20.40.60.81.0-0.20.00.20.40.60.81.01.21.41.6Stress (MPa)Strain (%) Before Aging After Aging0.000.010.020.030.040.05051015202530Stress (MPa)Strain (%) Before Aging After Aging0.000.020.040.060.080.10-50510152025303540Stress (MPa)Strain (%) Before Aging After Agingstrain rate was set at 50mm/min. The dimensions of the test samples were as follows: 115 mm

in length, 10 mm in width, and 4 mm in thickness. Figure 6.7 illustrates the impact of aging on

stress-strain curves obtained from these tensile tests. The variation in tensile stress levels can be

observed at different stages during the polymer aging process:

• During the initial phase, there was a notable increase in Young’s modulus and stress at the

breaking point after the post-cross-linking process. This observation indicates that the poly-

mer was not fully cross-linked prior to the aging process.

• In the subsequent phase, both Young’s modulus and stress at the breaking point exhibited a

decrease, ultimately returning to their initial values and stabilizing. This phenomenon can

be attributed to the rise in critical molecular weight and the reduction in cross-link density,

primarily as a consequence of chain scission within the polymer network.

Figure 6.7: Stress-strain curves of the polyurethane at different stages of aging; (a) at 85oC, (b) at
120oC [104].

The test procedures outlined by Ozmen and detailed in their study ( [90]) were carried out on PMR-

15 neat resin, with no pre-existing thermal aging. The experimental investigations encompassed

strain-controlled tests focused on tension-to-failure behavior, involving monotonic loading with

varying strain rates. Specifically, the influence of strain rate was examined in tension-to-failure

tests conducted at consistent strain rates of 106, s−1, 105, s−1, 104, s−1, and 103, s−1.

62

Stress (MPa)Strain (%)Aging Temperature = 85 oC Aging Temperature = 120 oC 00204060801001205101520253035400020406080100120510152025303540140Stress (MPa)Strain (%)Fitted CurveUnaged Experimental DataExperimental DataUnaged Fitted CurveAged for 62 weeksAged for 52 weeksAged for 9 weeksAged for 3 weeks13 weeks9 weeks4 weeks1 week2 daysFigure 6.8: Stress-strain curves for PMR-15 Specimens aged at 316oC in Argon Obtained in
Tensile Tests to Failure for aging time: (a) 50 hours, (b) 100 hours, (c) 250 hours, (d) 500 hours,
(e) 1000 hours [90].

63

012345051015202530 Loading at 10-3 s-1 Loading at 10-4 s-1 Loading at 10-5 s-1 Loading at 10-6 s-1Stress (MPa)Strain (%)Aging Temperature = 316oCAging Time = 50 h01234051015202530 Loading at 10-3 s-1 Loading at 10-4 s-1 Loading at 10-5 s-1 Loading at 10-6 s-1Stress (MPa)Strain (%)Aging Time = 100 hAging Temperature = 316oC0.00.51.01.52.02.50510152025 Loading at 10-3 s-1 Loading at 10-4 s-1 Loading at 10-5 s-1 Loading at 10-6 s-1Stress (MPa)Strain (%)Aging Temperature = 316oCAging Time = 250 h0.00.20.40.60.81.01.21.41.6051015202530 Loading at 10-3 s-1 Loading at 10-4 s-1 Loading at 10-5 s-1 Loading at 10-6 s-1Stress (MPa)Strain (%)Aging Temperature = 316oCAging Time = 500 h0.00.20.40.60.81.00246810121416 Loading at 10-3 s-1 Loading at 10-4 s-1 Loading at 10-5 s-1 Loading at 10-6 s-1Stress (MPa)Strain (%)Aging Temperature = 316oCAging Time = 1000 hFigure 6.9: Evolution of tensile strength as a function of thermal aging time (in days) of CFRP
samples [70].

6.5: Reaction Kinetics and Modeling Assumptions

To construct a predictive micromechanical model for the thermal aging of cross-linked polymers,

it is essential to characterize the underlying reaction kinetics and molecular architecture evolution

under inert (oxygen-free) conditions. Thermal degradation in such systems typically involves two

concurrent and competing mechanisms: chain scission and crosslinking. Both processes exhibit

strong temperature dependence and influence the material’s mechanical performance.

6.5.1: Temperature-Dependent Reaction Rate Modeling

In a finite element framework, temperature-dependent kinetics can be incorporated by defining

the reaction rate k(T ) as a function of temperature. Many thermal aging processes exhibit a non-

monotonic rate, increasing to a peak and then decreasing due to kinetic limitations or depletion of

reactive species. This behavior can be captured by a Gaussian-exponential product:

k(T ) = A · exp

(cid:19)

(cid:18)

− (T − µ)2
2σ2

· exp (−λ(T − Tpeak)) ,

(6.1)

64

010203040505510515520525530535540545550Tensile Strength (MPa)Aging Time (Days) Inert Oxidativewhere:

• A is the amplitude (maximum rate),

• µ is the mean temperature (center of the Gaussian),

• σ controls the spread (width of the Gaussian),

• λ is the decay constant post-peak,

• Tpeak is the temperature at which the rate is maximized.

This formulation enables a flexible representation of reaction kinetics and can be integrated into

constitutive models implemented within user-defined finite element subroutines.

6.5.2: Polymer Chain Distribution and Crosslinking Behavior

To model the evolving molecular structure under thermal aging, we consider the polymer chain

distribution using the Schulz-Zimm form, which generalizes the most probable distribution:

W (r) =

σσ
µΓ(σ)

(cid:17)

σ

(cid:16)

r
u

(cid:16)

(cid:17)

,

− σr
u

exp

(6.2)

where r is the chain length, µ is the number-average chain length, and σ characterizes distribu-

tion breadth.

The evolution of the chain weight fraction distribution due to random crosslinking is given by

Kimura’s expression:

W (r, x) =

"

r
µ2

·

σσ+1
Γ(σ + 1)

(cid:17)

σ−1

(cid:16)

r
u

∞X

+

αk

k=2

(cid:17)

(cid:16)

r
u

#

kσ+2k−3

(cid:16)

· exp

−(σ + 2µx)

(cid:17)

.

r
u

(6.3)

Crosslink density is defined as ρ = 2x, accounting for two units per linkage. The normalized

chain length distribution m(r, x) satisfies:

Z ∞

0

rm(r, x) dr = 1,

(6.4)

and relates to the weight and number fraction distributions as:

65

m(r, x) =

W (r, x)
r

=

N (r, x)
Pn

,

(6.5)

where Pn is the number-average chain length, and [Pr] is the concentration of chains of length

r.

The evolution of m(r, x) with respect to crosslinking density is governed by Saito’s integro-

differential equation:

∂m(r, x)
∂x

= −2rm(r, x) +

Z

r

0

h · (r − h) · m(h, x) · m(r − h, x) dh,

(6.6)

For an initial Schulz-Zimm distribution, the analytical solution becomes:

W (r, ρ) =

1
µ

(cid:20)

· (rµ)σ exp

−(σ + ρµ)

(cid:21)

·

r
µ

∞X

k=0

σ(k+1)(σ+1)(ρµ)k
(k + 1)!Γ[(k + 1)(σ + 1)]

(cid:19)

k(σ+2)

!

(cid:18)

r
µ

. (6.7)

6.5.3: Monte Carlo Simulation of Random Cross-Linking

To complement the analytical formulation, a Monte Carlo simulation is developed to stochastically

model the crosslinking process. This method enables the inclusion of history-dependent interac-

tions and spatial randomness not captured in deterministic models.

1. Initialization: Generate a chain population according to the Schulz-Zimm distribution.

2. Crosslink Selection: Randomly select chain ends for potential crosslinking based on prede-

fined probabilities or spatial proximity.

3. Update: Modify molecular structures to include T- or H-shaped branching after each suc-

cessful link.

4. Iteration: Repeat until a target crosslinking density ρ is reached.

This approach allows modeling of branching, gelation, and network connectivity, providing

valuable validation and insights beyond analytical solutions.

66

 
6.5.4: Evolution of Chain Length Distribution with Crosslink Density

Figure 6.10: Evolution of the chain length distribution W (r, ρ) as a function of increasing
crosslink density ρ, with parameters µ = 100, σ = 2. As ρ increases from 0 to 0.05, the peak of
the distribution shifts toward shorter chain lengths and the curve broadens, indicating the
progressive formation of branched and networked structures due to thermal aging.

Fig.6.10 illustrates how thermal aging influences the molecular architecture of a cross-linked

polymer network. At low crosslink density, the system is dominated by long linear chains. As ρ in-

creases, the average chain length decreases and the distribution broadens, indicating the emergence

of branched and interconnected structures consistent with network formation.

6.6: Constitutive Model

Constitutive models for rubber-like materials are typically classified as either phenomenological

or micro-mechanical. The distinction lies in the construction of the strain energy function W (F).

Phenomenological models express W using strain invariants (e.g., I1, I2), while micro-mechanical

models derive constitutive behavior from polymer physics [42, 106].

This section presents a physically motivated constitutive model for cross-linked polymers that

captures nonlinear inelastic behaviors, including an idealized Mullins effect, as well as chemically

induced aging. The idealized Mullins effect (also called stabilized softening) is modeled indepen-

dently of unloading/reloading paths.

67

6.6.1: Statistical Mechanics of Polymers

Experimental studies show that elastomers consist of long polymer chains oriented randomly in

three dimensions [117, 118]. Each chain comprises multiple monomers and is constrained by two

end linkages, often due to physical or chemical cross-linking (see Fig.6.11).

We assume the polymer network is a three-dimensional assembly of one-dimensional chains

distributed in all orientations (Fig.6.12). The macroscopic strain energy W (F) is calculated as the

directional average of the microscopic energies of these chains, integrated over the unit sphere.

A freely rotating chain (FRC) contains n Kuhn segments, each of length l, resulting in an end-

√

to-end distance

nl. The sub-network stretch in a given direction d is defined as:

p

λd =

dT FT Fd.

(6.8)

We use the following notation: scalars F , vectors F, second-order tensors F, and fourth-order

tensors F. A bar over a variable indicates normalization by Kuhn length l, e.g., ¯X = X
l .

Figure 6.11: Schematic of a single polymer chain consisting of n segments of length l, total
contour length nl.

Figure 6.12: (a) A single polymer chain with n segments and contour length nl; (b) a
micro-sphere structure with directional vector di ∈ R3.

68

nll3D Network Energy

The macroscopic strain energy function is computed via directional integration:

W (F) =

Z

S

1
4π

rψd dA ≈
ρd

kX

i=1

r ωiψdi.
ρdi

(6.9)

Here, ψd is the entropic energy of a single chain in direction d, and ωi are weights associated with

collocation directions on the sphere.

Single Chain Energy

The entropic energy of a single freely rotating chain is:

ˆψc(¯n) = ¯nKBT

Z

τ

0

ˆβ dˆτ + ψ0,

ˆβ =

(cid:20)

1 − 1 + τ 2

¯n

(cid:21)

β,

τ =

λd
√
¯n

,

β = L−1(τ ) ≈ τ

1 − τ

+ 2τ − 8
9

τ 2.

(6.10)

(6.11)

with:

Damage Modeling

We consider two damage mechanisms: mechanical and environmental. These are assumed to be

superimposable. Thus, chain deterioration is modeled as a function of temperature T , time t, and

deformation F.

Under mechanical loading, short chains debond from cross-links or surfaces. Following Mar-

ckmann et al. [81], damage is modeled via evolution of active chain length ¯nd and active chain

density ρd
r:

ˆnd = n0 exp(µλd

m),

r = ρr0 exp(−µλd
ρd

m),

(6.12)

where λd

m is the maximum experienced stretch in direction d, and µ is a damage sensitivity param-

eter.

For chemical aging (absent mechanical load), damage is modeled through changes in chain

69

density:

ρ0 = ρr∞ exp(µλm) + ρr0

· c ln c ·

c − µλm
µ0 exp(kt − Ea
RT )

,

(6.13)

where c is concentration, µ0 and k are fitting parameters, Ea is activation energy, and R is the

universal gas constant.

Constitutive Formulation

Assuming incompressibility:

det(F) = 1.

The first Piola-Kirchhoff stress tensor is:

P =

∂W (F)
∂F

− pF−T ,

where p enforces incompressibility.

The derivative of W is:

∂W
∂F

=

kX

j=1

ωjρdj
r

∂ψdj
t
∂λdj

· 1
2λdj

·

¯Cdj)
∂(dT
j
∂ ¯F

:

∂ ¯F
∂F

.

(6.14)

(6.15)

(6.16)

With simplification:

√

∂ψc
∂λd

=

¯nKBT ˆβ,

∂(dT ¯Cd)
∂ ¯F

:

∂ ¯F
∂F

= 2J−1/3 ¯F(d ⊗ d).

(6.17)

The final expression for the first Piola-Kirchhoff stress is:

P = KBT

kX

√

ρdi
r

i=1

¯ndiL−1

(cid:18)

(cid:19)

ωi
λdi

λdi
√
¯ndi

J−1/3 ¯F(di ⊗ di) − pF−T .

(6.18)

70

6.7: Finite Element Linearization

This section presents the algorithmic framework and variational formulation for implementing the

proposed constitutive model in a finite element (FE) environment. Both the second Piola-Kirchhoff

stress tensor and the consistent tangent modulus C are derived within a total Lagrangian finite

strain framework under the assumption of incompressibility and plane stress. The nonlinear system

is solved using the Newton-Raphson method, which requires linearization of the internal force

residual. The linearized stress-strain relation takes the form:

δS = C : δE,

(6.19)

where δS and δE are the increments of the second Piola-Kirchhoff stress and Green-Lagrange strain

tensors, respectively.

6.7.1: Stress Linearization

The second Piola-Kirchhoff stress is derived from the strain energy density function W(F), differ-

entiated with respect to the Green-Lagrange strain E = 1

2 (C − I). Under the incompressibility

constraint, the stress tensor is expressed as:

S =

∂W(F)
∂E

− pJC−1 = 2

∂W(F)
∂ ¯C

:

∂ ¯C
∂C

− pJC−1,

where ¯C = J −2/3C is the isochoric (deviatoric) part of the right Cauchy-Green tensor.

The projection tensor ∂ ¯C

∂C in the Lagrangian description is given by:

∂ ¯C
∂C

= J −2/3

(cid:18)

I − 1
3

(cid:19)

C ⊗ C−1

,

(6.20)

(6.21)

where I is the fourth-order identity tensor. Substituting Eq. (6.21) into Eq. (6.20) gives:

S = 2J −2/3 ∂W(F)
∂ ¯C

− 2
3

J −2/3

(cid:18)

∂W(F)
∂ ¯C

71

(cid:19)

: (C ⊗ C−1)

− pJC−1.

(6.22)

6.7.2: Out-of-Plane Stress Elimination

To enforce the plane stress condition, we set the out-of-plane Cauchy stress σ33 to zero. Solving

for the Lagrange multiplier p, we define:

Sn = 2

∂W
∂ ¯C

,

and obtain:

p = J −5/3Sn

33F 2
33

− 1
3

J −5/3Sn : C.

(6.23)

(6.24)

Substituting into the expression for S, the final form of the second Piola-Kirchhoff stress under

plane stress becomes:

(cid:0)

S = J −2/3

6.7.3: Tangent Modulus Derivation

Sn − Sn

33F 2

33C−1

(cid:1)

.

(6.25)

The derivation procedure for the tangent modulus is illustrated in Fig.6.13.

72

Figure 6.13: Tangent Modulus Derivation Procedure.

6.7.4: Plane Stress Simplification

For the plane stress case, the simplified derivation procedure and resulting tensor are shown in

Fig.6.14.

73

Figure 6.14: Plane Stress Simplification Procesure.

The detailed derivations of Sn and Cn are provided in the Appendix.

74

Figure 6.15: A schematic figure of the hierarchical micro-macro multi-scale approach.

Algorithm 1: Hierarchical Atomistic-Continuum Multi-Scale Method

1. Initialization

• Discretize the domain with a finite element mesh; assign an atomistic RVE to each

integration point.

• Loop over integration points:

– Set initial deformation gradient F = I.
– Characterize each atomistic RVE; store initial atomic coordinates.
– Compute first Piola-Kirchhoff stress P and tangent modulus CPF from the RVE.
– Upscale P and CPF to the macroscopic integration point.
– Compute second Piola-Kirchhoff stress S and tangent modulus CSE via

continuum mechanics.
• End loop over integration points.

2. For each increment (n → n + 1):

• Initialize variables from the last converged state.
• Apply external loads or prescribed displacements.

3. For each iteration (i → i + 1):

• Assemble global stiffness matrix K; solve for incremental displacement u.
• Update nodal positions accordingly.
• Loop over integration points:

– Recalculate deformation gradient F in current and reference configurations.
– Apply F to RVE under periodic boundary conditions.
– Recompute P and CPF at the atomistic level.
– Upscale results to the continuum integration point.
– Update S and CSE.

• End loop over integration points.
• Compute internal force vector and residual; check for convergence.

4. Check for convergence:

• If converged, proceed to the next increment.
• Otherwise, repeat the next iteration.

75

ΓuΓtInitial configurationtuFDeformed configurationC, SxxxxxxxxxxxxxxxxGaussIntegrationpointMulti-scale model6.8: Validation and Results

Fig.6.16 presents the stress-strain curve from monotonic tension-to-failure tests conducted by [35].

These experiments were carried out at temperatures of 288◦C and 316◦C with four different aging

durations: unaged (conducted at room temperature without aging), aged for 100, 250, and 1000

hours. The unaged experimental results were utilized to calibrate the parameters of the constitutive

model developed. Fig.6.16 shows the model predictions for the stress-strain curve of PMR-15

under these different aging conditions. Comparison of the predicted stress-strain curve with the

experimental data demonstrated an error of less than 10%, validating the accuracy and effectiveness

of the model.

Figure 6.16: Stress-strain curve obtained in monotonic tension to failure tests for PMR-15 aged in
Argon in (a) 288oC, and (b) 316oC [35].

Fig.6.17 compares deformation results obtained from two implementations developed in this

work: (a) a custom Abaqus subroutine and (b) a MATLAB script. The MATLAB code was origi-

nally inspired by the framework presented in [9], which modeled hydrolysis and thermo-oxidative

aging. However, the code was significantly adapted and extended to model inert aging conditions

relevant to our study. The strong agreement between both implementations confirms the correct-

ness and consistency of our computational models.

76

ExperimentalDataStress (MPa)Strain (%)01234567805101520253035Aging Temperature = 288 oC Fitted CurveUnaged Aged for 100 HourAged for 250 HourAged for 1000 HourAging Temperature = 316 oC 246810120Stress (MPa)Strain (%)05101520Experimental DataUnaged Aged for 50 HourAged for 250 HourAged for 500 HourFigure 6.17: Deformation results comparison: (a) Abaqus subroutine implementation (UMAT)
and (b) MATLAB simulation.

The UMAT subroutine (the UMAT code, written in Fortran, is provided in the Appendix) was

derived using an open-source library called Tensor Toolbox, which was implemented in the Fortran

programming language. This library provides essential tensor operations that are crucial for the

derivation of the tangent modulus in the context of nonlinear material models. The Tensor Toolbox

enables efficient and accurate handling of tensor calculations, which are integral to computing the

material’s stress-strain relationship and, consequently, the tangent modulus.

The tangent modulus is essential due to the nonlinear behavior of materials like rubber un-

der large deformations. In nonlinear finite element analysis (FEA), iterative solvers, such as the

Newton-Raphson method, rely on the tangent stiffness matrix to update the material’s behavior at

each iteration and ensure accurate convergence of the solution.

6.8.1: Numerical Simulation Results

Two sets of numerical results on different inert thermal aging temperature are presented here to

demonstrate the finite element simulation of the mechanical response of cross-linked polymers.

The objective of this simulation is to show the capability of the numerical model to accurately

predict the response of the polymer material under different inert aging conditions. These two

77

U, Magnitude+0.000e+00+3.340e+00+6.680e+00+1.002e+01+1.336e+01+1.670e+01+2.004e+01+2.338e+01+2.672e+01+3.006e+01+3.340e+01+3.674e+01+4.008e+01abexamples are a classical dog-bone specimen under prescribed displacements along −x direction.

Due to the thin film assumption, the plane stress assumption is adopted.

Figure 6.18: Evolution of Second Piola–Kirchhoff Stress in Epoxy Resin Under Inert Thermal
Aging at 60◦C and Mechanical Loading (UMAT Results).

Figure 6.19: Evolution of Second Piola–Kirchhoff Stress in Epoxy Resin Under Inert Thermal
Aging at 80◦C and Mechanical Loading (UMAT Results).

78

Viewport: 1 ODB: C:/temp/Dogbone_T333_t1.odb(Avg: 75%)S, Mises+8.356e+00+2.132e+01+3.429e+01+4.726e+01+6.022e+01+7.319e+01+8.616e+01+9.913e+01+1.121e+02+1.251e+02+1.380e+02+1.510e+02+1.640e+02Viewport: 1 ODB: C:/temp/Dogbone_T333_t30.odb(Avg: 75%)S, Mises+7.794e+00+1.992e+01+3.204e+01+4.416e+01+5.629e+01+6.841e+01+8.053e+01+9.266e+01+1.048e+02+1.169e+02+1.290e+02+1.412e+02+1.533e+02ab(Avg: 75%)S, Mises+6.589e+00+1.689e+01+2.719e+01+3.749e+01+4.779e+01+5.809e+01+6.839e+01+7.870e+01+8.900e+01+9.930e+01+1.096e+02+1.199e+02+1.302e+02cCHAPTER 7: CONSTITUTIVE

MODELING OF DIFFUSION-LIMITED

OXIDATION (DLO) COUPLED WITH A

LARGE DEFORMATION THEORY FOR

POLYMER DEGRADATION

7.1: Abstract

The impact of oxidation on polymer degradation is a crucial aging process, particularly when it’s

confined to the surface of the sample, often referred to as diffusion-limited oxidation (DLO). DLO

entails a complex interplay between oxygen absorption, diffusional transport, and oxygen’s reac-

tion within elastomers. In this study, a novel kinetics-based model for oxygen absorption is devel-

oped and rigorously validated against a range of experimental datasets. Furthermore, the diffusion-

reaction equation is extended into three dimensions and solved using the Alternating Direction

Implicit (ADI) method. Various reaction rate functions are considered to account for chain scission

and network reformation reactions, enabling the description of non-uniform degradation patterns.

The enhanced model is employed to simulate heterogeneous oxidation and to showcase how

different contributing factors influence the oxidation behavior of nitrile rubber. Notably, the out-

comes of our proposed model align closely with empirical data concerning the degree of oxidation

in polymers. Additionally, a constitutive model is formulated to incorporate the intricate interplay

between diffusion, chemical reactions, and large deformations in polymers. This model is seam-

lessly integrated into a finite element framework to illustrate its predictive capabilities concern-

ing the effects of diffusion-limited oxidation (DLO) on the mechanical properties of cross-linked

polymers. It allows for a comprehensive numerical analysis of the coupled diffusion-reaction and

mechanical behaviors exhibited by polymers undergoing DLO.

79

7.2: Introduction

Polymers employed in engineering applications, including those used in aerospace and turbine

engines, possess a finite lifespan primarily due to the effects of degradation. This degradation

predominantly results from environmental factors that are present throughout the manufacturing

and service life of these materials. The three primary factors responsible for the heterogeneous

aging of polymers are ultraviolet (UV) radiation, ozonation, and oxidation [6, 60, 67]. To ensure

the effective design of polymer materials for engineering applications in extreme conditions, it be-

comes crucial to accurately predict their long-term stiffness and strength, as well as to model the

mechanisms leading to degradation and eventual failure. In many cases, polymers experience a

substantial degradation in their performance and mechanical properties due to the effects of aging,

which becomes particularly pronounced in highly specialized applications [116]. High-temperature

oxidation or thermo-oxidative aging in polymers serves as a prominent illustration of aging mecha-

nisms. In this scenario, oxygen permeates elastomeric materials and engages in chemical reactions

with active polymer sites. Consequently, layers of oxides start to develop in proximity to the ex-

posed surface, fostering conditions conducive to the initiation and spread of micro-cracks within

the material. This phenomenon can exert a substantial influence on the mechanical characteristics

and overall structural integrity of the polymer as time progresses. [27]. These micro-cracks not only

propagate but also exacerbate the oxidation process within the material, ultimately culminating in

the failure of the material [15, 16, 25]. Upon exposure to air, elastomers interact with oxygen at

their surface boundary, a process governed by Henry’s law. The adsorbed oxygen then proceeds to

diffuse into the material, initiating polymer reactions that ultimately lead to polymer degradation.

A substantial body of recent research has delved into examining the alterations in the polymer’s

network structure stemming from thermo-oxidation [26, 32, 103]. The aging of polymers under

thermal oxidative conditions, especially under accelerated conditions, frequently encompasses com-

plex and non-uniform degradation mechanisms associated with diffusion-limited oxidation (DLO)

phenomena [29, 50]. The shape of DLO highly depends on the oxidation rate, oxygen diffusiv-

80

ity, thickness of the sample, and the presence of antioxidants [14, 48]. Heterogeneous oxidation

in elastomers is driven by diffusion reaction, resulting in changes in the mechanical properties of

elastomers [24, 78, 94, 95, 108]. Chain scission and cross-linking processes play pivotal roles in re-

shaping the polymer network, resulting in significant alterations in the material’s behavior during

oxidation. These changes induced by oxidation are extensively documented in recent literature and

include phenomena such as an increase in modulus, a decrease in glass transition temperature, and

reduced failure strain, among others [11, 39–41, 82]. The stress distribution in the elastomers can

also be affected by oxidation. Hence, considering DLO would enable us to accurately predict the

stiffness of polymers during aging. Therefore, the uncertainty in the performance prediction of the

sample can be avoided [49, 75]. Modeling heterogeneous oxidation in elastomers poses significant

challenges due to the intricate interplay of physical, chemical, and thermo-mechanical mechanisms

that evolve over time during the aging process. In the context of polymer oxidation, a comprehen-

sive understanding involves a set of six elementary reactions. Extensive research has also been

conducted to investigate the impact of chain reactions on polymer oxidation, resulting in the devel-

opment of kinetic models that aim to describe these oxidative processes in polymers [7, 8, 28, 111].

For a thick elastomer specimen, the DLO model was introduced in [96], which described the depth

and extent of DLO in various circumstances. The adsorption and diffusion of oxygen in the elas-

tomer play a major role in heterogeneous thermo-oxidative aging [56,64,102]. The absorption pro-

cesses and the subsequent diffusion of oxygen molecules into elastomers have been investigated

experimentally [57]. In this study, the paper reports findings on the solubility of oxygen and the

degree of oxidation throughout the thickness of Nitrile Rubber during the aging process. Addition-

ally, a model has been developed to effectively represent the absorption of oxygen and the extent

of oxidation in elastomers. This model takes into consideration variables such as temperature and

aging time to provide insights into these critical aspects of polymer behavior during aging [58].

In order to predict how polymers respond to diffusion-limited oxidation, researchers have in-

tegrated the balance of momentum with the diffusion-reaction equation. This approach allows for

an in-depth exploration of how non-uniformly oxidized polymers behave as they age. The concept

81

of coupling multiple physical phenomena in material constitutive modeling has been extensively

studied in the field of mechanics. For instance, a coupled theory was developed to explain high-

temperature oxidation in metallic alloys. This illustrates the interdisciplinary nature of materials

science and mechanics in understanding complex phenomena [79]. The coupled stress-diffusion re-

sponse in polymer gels has been extensively studied by [19–21,37,51]. The reaction-diffusion cou-

pling in photo-responsive gels is also noteworthy work [113]. Several works have also been done

on coupled photo-chemical reaction-large deformation theory as [33, 100, 112]. The non-uniform

degree of oxidation leads to variations in mechanical properties and results in an inhomogeneous

distribution of stress and strain within the oxidized polymer [83, 89, 92].

This study refines the oxygen absorption model and validates its ability to replicate experimen-

tal findings. We expand the diffusion-reaction equation into three dimensions and elaborate on

the solution method using the ADI numerical approach. Novel reaction rate functions for chain

scission and network reformation, accounting for oxygen concentration and oxidation levels, en-

hance predictive accuracy. Our DLO modeling, employing these reaction rate functions, exhibits

improved correlation with experimental data. Additionally, we develop a constitutive model to

explore oxidation progression under mechanical loads and its impact on polymer mechanical prop-

erties during aging. This model incorporates the interplay between diffusion, chemical reactions,

and polymer deformation. We then implement the model within a finite element framework to ana-

lyze how polymers respond to DLO. Finally, we investigate the advancement of oxygen diffusion

and the evolution of oxidation extent in cases of stress-coupled oxidation.

The paper is organized as follows. In section 7.3, the modified oxygen absorption model is pre-

sented. Section 7.4 explains the extended diffusion-reaction equation and ADI solution procedure.

The results are analyzed in detail to explain the model’s capability for the simulation of oxidation

across the thickness of polymers. In section 7.5, the developed constitutive model for coupled

stress-oxidation is introduced, and the model is further implemented in a FE setting by writing a

user element subroutine (UEL). In the end, the concluding remarks are presented.

82

Figure 7.1: Oxygen Absorption Rate vs. Temperature.

7.3: Oxygen Absorption

Oxygen is soluble in polymers, and the permeation process can be split into solubility and diffusivity

[107]. The solubility (also named oxygen uptake or absorption) describes the initial oxidation step

and is provoked by the oxygen gradient between the polymer and the ambient environment as

the driving force. Oxygen uptake is a function of temperature, ambient oxygen pressure, and the

aging duration or rather the aging progress [101]. Fig.7.1 Illustrates the varying oxygen absorption

rates of a polymer at different temperatures. This section uses the experimental data of oxygen

absorption to verify the proposed oxygen uptake model. The experiment has been conducted on

nitrile rubber, and provides the experimental results [59]. For details of the experiment procedure

and the material used, see [59]. Oxygen uptake was measured by a respirometer and Differential

Oxygen Analyzer. The samples were aged in sealed chambers at multiple elevated temperatures

before the oxygen concentration of the contained air was measured. The difference in the amount

of oxygen is measured before and after exposure, which is associated with elastomer absorption

[55, 56]. The oxygen absorption model of oxidized material is developed to analyze the oxygen

solubility in polymers. In aerobic aging, solubility is a significant factor mainly governed by the

partial pressure of oxygen pox. The oxygen absorption rate,

˙mO2, at the exterior surface can be

83

calculated using Henry’s law which describes it with respect to solubility rate coefficient S and the

partial pressure pox, as given below:

˙mO2 = S · pox,

(7.1)

where S is a material parameter describing the solubility rate coefficient and is governed by chain

scission and network reformation processes during exposure to air. The novelty of the proposed

oxygen absorption model is that aging processes are divided into chain scission and network ref-

ormation parts to see how these processes affect the oxygen uptake in the superficial layer of the

elastomer.

Another novelty of the proposed model is avoiding the complete stop of oxygen absorption.

According to the work done by [58], oxygen diffusion is not prevented even for the totally oxidized

part of the elastomer. So, even for the oxidized boundary, the oxygen is absorbed, diffusing into the

elastomer in favor of the reaction. But, in the oxygen absorption model in [58], it was assumed that

the fully aged superficial layer inhibits any uptake of oxygen molecules. Hence, in this work, the

β parameter is defined, which is the model parameter and describes the maximum percentage of

reduction of oxygen uptake by aging. Therefore, a complete oxygen uptake stop is avoided since it

can only be reduced to a certain minimum of the initial value. Oxygen adsorbing at the elastomer’s

surface declines with aging and can be described by:

˙mO2 = pox · S0 · exp(

(cid:18)

(cid:18)

1 − β

−Es
RT

)

psci + pref
2

(cid:19)(cid:19)

,

(7.2)

where S0 describes the initial solubility rate of the matrix, Es the activation energy of solubility, R

the universal molar gas constant, and T the absolute temperature. The aging parameters describe

aging progress, whereas the value 0 describes virgin material, and 1 stands for a completely aged

material. The accompanying evolution equation of chain scission and network reformation that

describe the aging progress can be expressed as follows:

˙psci = νsci

p

· exp(

(cid:0)

−Esci
p
RT

)

1 − psci

(cid:1)

2 ,

(7.3)

84

˙pref = νref

p

· exp(

−Eref
p
RT

(cid:0)

)

1 − pref

(cid:1)2

,

(7.4)

where νp represents a proportionality factor and can be imagined as a temperature-independent

velocity of the aging progress. Ep is the activation energy of the aging process. The oxygen uptake

model introduced in [58] did not give any information about the effect of the size and geometry of

the specimen on oxygen absorption. In this paper, the L parameter is defined, which is the surface-

volume ratio of the sample. The basic solubility rate is considered a function of the surface-volume

ratio in the model. An experiment was conducted on nitrile rubber, and the total mole of oxygen

per gram of the sample weight was reported as a function of aging temperature and duration [59].

The bar-shaped samples (55x20x2mm) were aged at temperatures of 100, 120oC and for periods

up to 72 h. The above model is used to describe the experimental results. The experimental data is

converted to the total mole of oxygen per surface-volume ratio to investigate the capability of the

developed model. Therefore, the evolution equation is solved by using ode45 MATLAB non-stiff

differential equations solver and fitted to the experimental data.

Figure 7.2: Oxygen consumption of the nitrile rubber as a function of aging temperature and
duration.

85

Fig. 7.2 shows that the modeling approach can replicate the experimentally obtained oxygen

absorption during aging. With ongoing exposure, the aging parameters increase according to the

evolution equations, and the solubility of the sample is reduced. Hence, the oxygen absorption rate

is decelerated during aging. The parameters used for the modeling shown in Fig. 7.2 are chosen as

declared in Table 7.1.

pox(Kpa)
S0 [mol/L/s/Pa]
Es
Esci
p

2
1.8 × e−6
5.6 × e4
47000

R[J/mol/K]
νsci
p
νref
p
Eref
p

8.3145
10
10
47000

Table 7.1: Parameters for modeling of oxygen absorption behavior.

Nevertheless, measuring the oxygen uptake provides no information about the distribution of

oxidation in the elastomer. Thus, no conclusion can be drawn about any kind of heterogeneity

of oxidation. Hence, besides solubility, diffusivity plays an important role when considering the

process of oxidative aging. Therefore, the diffusion process itself should be considered more in

detail.

7.4: Diffusion reaction equation

As stated in the previous section, the chemical aging of an elastomer by oxidation is generally

defined as a process involving several steps, including absorption, diffusion, and chemical reactions.

After oxygen molecules are adsorbed due to the affinity of the polymer’s surface for O2 molecules,

a concentration gradient occurs in the exterior layer before the oxygen molecules subsequently

diffuse into the material’s interior. The most common approach to describe this kind of diffusion

process is given by Fick‘s second law as follows.

∂C
∂t

= ∇ · (Dox∇C) ,

(7.5)

where oxygen concentration C and diffusivity are formulated as a function of the location in the

material and aging time. Oxygen molecules react chemically with elastomers while it is diffusing

86

into the sample. Hence, based on the set of reactions for oxidation [99] and the general mass balance

law for oxygen concentration, we can write,

∂C
∂t

= ∇ · (Dox∇C) − rsci

ox + rref
ox .

(7.6)

Oxygen is consumed in the propagation reaction and produced in the termination reaction. rsci

ox and

rref
ox are the rate of consumption or production of oxygen in chain scission and network reformation

reactions, respectively. The reaction rate functions are a function of aging temperature and oxygen

concentration. Both chemical aging reactions prevent oxygen from further diffusion since they

create a chemical compound with a network of elastomer molecules. Chain scission and network

reformation are treated individually and are described by an Arrhenius approach. The reaction

rate functions are a function of the elastomer’s oxygen concentration, temperature, and extent of

oxidation. If there is no oxygen, no oxidative reaction can occur. Consequently, the progression of

aging can be described by the following reaction rate functions.

rsci
ox = F (C)rsci

o exp

ox = F (C)rref
rref

o

exp

(cid:18)

(cid:18)

(cid:19)

−Esci
r
RT

(1 − psci)n,

(cid:19)

−Eref
r
RT

(1 − pref )n,

(7.7)

(7.8)

are constant factors influencing the velocity of the reactions. F (C) is a function

where rsci

o and rref

o

to describe the reaction rate dependency on oxygen availability. The exponent n is a parameter

used to adjust the effect of the extent of oxidation on reactivity. Thus, the reactivity reduces with

the progressing aging process since the amount of convenient reaction partners decreases. The

evolution equations that describe the progress of aging can be expressed as [58]:

˙psci = νsci

p

rsci
ox
rsci
ox0

,

87

(7.9)

˙pref = νref

p

rref
ox
rref
ox0

.

(7.10)

Diffusivity is a function of temperature and is also affected by aging. Here, it is assumed that

both types of oxidative reactions, chain scission, and network reformation, reduce the diffusivity

of the elastomer in the same way. Thus, the diffusion coefficient Dox is formulated as a function

of the aging parameters chain scission psci and network reformation pref . The diffusion coefficient

is written as:

Dox = D0 · e−ED/Rθ.ξ, with ξ = ιpsci · pref .

(7.11)

The model parameter ι is used to adjust the influence of the aging parameters on the diffusion

coefficient. In order to get a better understanding of heterogeneous aging, the diffusion-reaction

equation is extended for three dimensions as the following equation:

∂C
∂t

= Dox

∂C 2
∂2x

+

∂Dox
∂x

∂C
∂x

+ Dox

∂C 2
∂2y

+

∂Dox
∂y

∂C
∂y

+ Dox

∂C 2
∂2z

+

∂Dox
∂z

∂C
∂z

− rsci

ox + rref

ox . (7.12)

7.4.1: Numerical solution of the diffusion-reaction equation

The ADI method is used to solve the reaction-diffusion equation in three dimensions. In the ADI

method, every time step is subdivided into three steps. This numerical procedure works so that a

3D problem is transformed into three 1D problems which can be solved with an implicit approach.

The numerical solution can be simplified to three tridiagonal systems per time step. First, the con-

centration is calculated in the x-direction, second in the y-direction, and then in the z-direction [77].

The Douglas higher order ADI scheme is utilized for the accurate numerical solution to solve the

3D diffusion-reaction equation [36]. The second-order finite difference method is also used for dis-

cretizing the first derivative of oxygen concentration and diffusion with respect to directions. The

first step provides the oxygen concentration C in the time step m + 1/3. The following equations,

i, j, and k, are used as variables to describe location dependency in the x, y, and z directions.

88

(cid:16)

m+ 1
3
i,j,k

C

− C m
i,j,k

(cid:17)

∆t

(cid:20)

(cid:16)

1
2

δ2
x

=

D
h2

m+ 1
i,j,k + C m
3

i,j,k

C

(cid:17)

(cid:0)

(cid:1)

(cid:0)

C m

i,j,k

+ δ2
z

C m

i,j,k

+ δ2
y

(cid:21)

(cid:1)

+

Di+1,j,k − Di−1,j,k
2∆x

m+ 1
3
i+1,j,k

C

− C

+

+

Di,j+1,k − Di,j−1,k
2∆y
Di,j,k+1 − Di,j,k−1
2∆z

m+ 1
i−1,j,k + C m
3
4∆x
− C m
2∆y

i,j+1,k

C m

i+1,j,k

i,j−1,k

C m

i,j,k+1

− C m
2∆z

i,j,k−1

− C m

i−1,j,k

− rsci,m

ox,i,j,k + rref,m
ox,i,j,k.

(7.13)

In the equation above, δ2 is the second-order difference operator. The concentration C m+1/3

i,j,k

is the first step solution used by the ADI method and calculated at every grid point of the area

with i, j, k = 1, ..., N . N stands for the number of grid points in the x, y, and z directions and

h represents the constant distance between two points in the mesh grid. The Eqn. 7.13 can be

written in the tridiagonal matrix and right-hand-side vector. Subsequently, an LU-factorization of

the tridiagonal matrix is executed, which can be solved by the Thomas algorithm. The solution

of Eqn. 7.13 contains a N tridiagonal equation system. A constant oxygen concentration C0 is

assumed in the ambient environment. Hence, Dirichlet boundary conditions are used. After the

first step, oxygen concentration C m+1/3

i,j,k

is obtained, the ADI method uses a second step to calculate

the next step oxygen concentration C m+2/3

i,j,k

.

(cid:16)

m+ 2
3
i,j,k

C

− C m
i,j,k

(cid:17)

∆t

(cid:20)

(cid:16)

1
2

δ2
x

=

D
h2

m+ 1
i,j,k + C m
3

i,j,k

C

(cid:17)

(cid:16)

+ δ2
y

m+ 2
i,j,k + C m
3

i,j,k

C

(cid:17)

(cid:0)

C m

i,j,k

+ δ2
z

(cid:21)

(cid:1)

+

+

m+ 1
3
i+1,j,k

C

Di+1,j,k − Di−1,j,k
2∆x
Di,j+1,k − Di,j−1,k
2∆y
Di,j,k+1 − Di,j,k−1
2∆z

m+ 2
3
i,j+1,k

+

C

− C

− C

i+1,j,k

i,j+1,k

m+ 1
i−1,j,k + C m
3
4∆x
m+ 2
i,j−1,k + C m
3
4∆y
− C m
2∆z

i,j,k−1

C m

i,j,k+1

− C m

i−1,j,k

− C m

i,j−1,k

− r

sci,m+ 1
ox,i,j,k + r
3

ref,m+ 1
3
ox,i,j,k

.

(7.14)

Analogously, the ADI method uses the final step to calculate the final oxygen concentration

C m+1.

89

(cid:1)

(cid:0)

C m+1
i,j,k

− C m
i,j,k
∆t

=

(cid:17)

(cid:16)

+ δ2
y

m+ 2
i,j,k + C m
3

i,j,k

C

(cid:17)

(cid:0)

+ δ2
z

C m+1

i,j,k + C m

i,j,k

(cid:1)i

i,j,k

m+ 1
i,j,k + C m
3
m+ 1
3
i+1,j,k

C

− C

m+ 2
3
i,j+1,k

C

− C

h

(cid:16)

C

δ2
x

D
2h2
Di+1,j,k − Di−1,j,k
2∆x
Di,j+1,k − Di,j−1,k
2∆y
C m+1

i,j,k+1

+

+

i+1,j,k

i,j+1,k

− C m

i−1,j,k

− C m

i,j−1,k

m+ 1
i−1,j,k + C m
3
4∆x
m+ 2
i,j−1,k + C m
3
4∆y
− C m

+

Di,j,k+1 − Di,j,k−1
2∆z

− C m+1

i,j,k−1 + C m
4∆z

i,j,k+1

i,j,k−1

− r

sci,m+ 2
ox,i,j,k + r
3

ref,m+ 2
3
ox,i,j,k

.

(7.15)

After finishing the calculation of a time step, the oxygen concentration Ci,j,k with i, j, k =

1, · · ·, N is obtained, which is used to update the aging parameters psci and pref , and subsequently

the diffusion coefficient Dox at every grid point. The aging parameters are updated according to the

evolution equations. The ADI method is unconditionally stable, and its accuracy is of second order

with respect to time and space [36]. The evolution equations are solved numerically by means of

the implicit finite difference method, which reads for the chain scission as follows. The numerical

solution for the network reformation is implemented similarly and is therefore not shown here.

psci,m+1
i,j,k

− psci,m
i,j,k

∆t

= νsci
p

rsci
ox
rsci
ox0

.

(7.16)

The mentioned numerical method for solving the 3D diffusion reaction equation is used to

visualize a cross-section in the sample during aging. Three-dimensional plotting gives a compre-

hensive visual impression of thermo-oxidative aging of a bar-shaped sample’s cross-section and

DLO. Model parameters are chosen as two-dimensional equations published in [58]. For illustra-

tion, a cube is considered and is aged up to 1000 seconds at 393 Kelvin. This aging time cannot be

compared to the real aging behavior of elastomers. It is just a fictitious value and is used only to

show the basic functionality of the model. Fig. 7.3 gives an example of the spatial development of

the aging parameter psci for various depths of the cube.

90

Figure 7.3: 3D Plot of the Aging Parameter ψci Computed Using the Developed MATLAB Code
for Varying Thicknesses and Aging Durations — (a) After 200 s, (b) After 1000 s.

In the ongoing process, oxidative reactions consume oxygen, decreasing the diffusivity. After

a while, the surface layer is fully aged, and hence diffusivity is almost diminished. The diffusion

coefficient decreases so that diffusion is reduced to a minimum by the progress of aging. For such

high temperatures, DLO inhibits oxygen from reaching the core region of the sample. Therefore,

heterogeneous oxidation occurs over the thickness of the sample.

7.4.2: Reaction rate function schemes

The model introduced is applied to simulate the heterogeneous aging of peroxide-cured NBR in

an oxygen-containing environment. An experiment has been conducted on nitrile rubber, and the

carbonyl-to-nitrile ratio was reported during aging. To describe the degree of oxidation from both

the theoretical and experimental side, it was assumed that the carbonyl nitrile ratio is equalized with

the aging parameter psci [57]. The modeling approach tries to consider the experimental findings.

Therefore, the reaction terms are modified to adjust the influence of aging on reactivity. Here,

for better correlation with experimental data, some reaction rate functions for chain scission and

network reformation processes are proposed, and their capability is investigated. The reaction-

diffusion equation with reaction terms for chain-scission and network reformation is used, and the

model is adjusted and fitted to the data obtained by FTIR investigation in [57]. The reaction rate

functions consider the non-linear contribution of the aging parameters psci and pref and oxygen

91

concentration. The aim of investigating the proposed reaction rate functions is to find the best

reaction functions that can accurately describe the extent of oxidation in elastomers. Here are three

schemes for reaction rate functions for chain scission reactions based on an Arrhenius approach.

(cid:18)

rsci
ox =

C
C0

rsci
0 exp

−Esci
r
RT
(cid:18)

(cid:19)

(1 − psci)2,

(cid:19)

(1 − psci),

rsci
ox =

αC
1 + βC

rsci
0 exp

−Esci
r
RT

(7.17)

(7.18)

(cid:21) (cid:20)

(cid:20)

2αC
1 + βC

rsci
ox =

1 −

βC
2(1 + βC)

(cid:21)

(cid:18)

rsci
0 exp

(cid:19)

−Esci
r
RT

(1 − psci)1/2.

(7.19)

The reaction terms for network reformation rref

ox are formulated analogously. The α and β

are the parameters to consider the effect of oxygen concentration on the reaction rate functions.

The model can describe the oxidative profiles, whereas equating the aging parameter psci and the

carbonyl nitrile ratio must be interpreted cautiously.

92

Figure 7.4: Comparison of the experimentally measured degree of oxidation of nitrile rubber aged
at 393 K with model predictions using three different reaction rate schemes.

The first scheme of reaction rate functions was proposed by Herzig et al. [57]. Fig. 7.4 illus-

trates that the proposed functions for reaction rate have a better correlation with experimentally

measured degree of oxidation. It is clear that several factors like the type of raw material, oxygen

concentration, and presence of antioxidants or temperature influence the degree of heterogeneity.

The parameters used for modeling nitrile rubber oxidation shown in Fig. 7.4 are chosen as stated

in Table 7.2.

D0[m2/s]
rsci
0 [mol/l/s]
Esci
[J/mol]
r
νsci
p

1 × e−5
7

rref
0
2.3 × e4 Eref

r

ED[J/mol]

[mol/l/s]
[J/mol]
νref
p

0.5

3.5 × e4
7
2.4 × e4
0.5

Table 7.2: Parameters used for modeling of carbonyl-nitrile-ratio of the NBR.

93

7.4.3: Linkage of oxygen absorption and diffusion-reaction behavior

After oxygen is absorbed, there are four possibilities for subsequent progression: Oxygen is dis-

solved in the elastomer without any chemical reaction, it reacts with the elastomer (oxidative reac-

tions), it reacts with the antioxidant, or it is adsorbed by the elastomer again. The amount of ab-

sorbed oxygen can be calculated by summing up these cases during aging. The currently dissolved

oxygen in the elastomer is added to the oxygen consumption in the chain scission reaction, the

oxygen production in the network reformation reaction, and binding oxygen molecules by means

of antioxidation. Thus, the amount of oxygen absorbed from the beginning of exposure until aging

time t can be calculated as an output of the diffusion-reaction model. That allows us to compare the

calculated amount of absorbed oxygen with the experimental data obtained in Section 7.3. There-

fore, the following relationship can calculate the total amount of absorbed oxygen.

y

Z

y

t

(cid:0)

mox = mdis + mre =

Co2dv +

V

0

V

rsci
ox

− rref

ox + rant
ox

(cid:1)

dv.

(7.20)

The currently dissolved amount of oxygen is obtained by the summation of the present oxygen

at every spatial element of the whole sample. The total amount of oxygen reacted after absorption is

calculated by summing up the consumption, production, and reaction of oxygen with antioxidants

of every previous time step since the start of exposure. The three reaction rate function schemes

for chain scission and network reformation are implemented for calculating the Eqn. 7.20. The

amount of absorbed oxygen can be plotted over the aging duration and represents a by-product of

the model. The efficacy of the schemes for reaction rate functions is investigated by comparing the

model results with oxygen absorption tests, illustrated in Fig. 7.5.

94

Figure 7.5: Comparison of oxygen absorption profiles over aging time using three reaction rate
models.

The primary reason for the disparities between the experimental data and the simulation out-

comes is the existence of antioxidants. In the initial stages of exposure, a portion of the oxygen

reacts with antioxidants, an effect that is not currently accounted for in the diffusion-reaction model.

However, over time, the antioxidants will gradually lose their effectiveness. To improve the align-

ment between simulation results and experimental data, it is crucial to incorporate the reaction

between oxygen and antioxidants into the model, thereby providing a more accurate representation

of the real-world conditions.

7.5: Constitutive model for stress-coupled oxidation

In this section, the consistent chemo-mechanically coupled theory is developed for polymer oxida-

tion. The developed model will most likely help to reach an advanced level in describing elastomer

aging. For example, the effect of mechanical loads on diffusion or the extent of oxidation, chain

scission, and network reformation, in general, can be analyzed more deeply.

We identify such a macroscopically homogeneous body B with the region of space it occupies

in a fixed reference configuration and denote by X an arbitrary material point of B. A motion

95

of B is then a smooth one-to-one mapping X = χ(X, t) with deformation gradient, velocity, and

velocity gradient is given by:

F = ∇χ,

v = ˙χ,

L = gradv = ˙FF−1.

(7.21)

We use the standard notation of modern continuum mechanics [52]. Specifically: ∇ and Div denote

the gradient and divergence with respect to the material point X in the reference configuration;

grad and div denote these operators with respect to the point x=χ(X, t) in the deformed body; a

superposed dot denotes the material time-derivative. We consider a multiplicative decomposition

of the total deformation gradient based on the large deformation theory of polymers as follows:

F = FeFp,

(7.22)

where Fe is the elastic component of the deformation gradient and Fp is the irreversible plastic

component. Referring to Eqn. 7.22, the determinant J can be decomposed into,

J = JeJp,

Where, Je = det(Fe),

Jp = det(Fp).

(7.23)

Further, the velocity gradient is defined as L = ˙FF−1, which can be subdivided into the elastic part

and inelastic part as Le = ˙FeF−e and Lp = ˙FpF−p. The elastic and plastic deformation tensor can

be further defined from the standard continuum mechanics as symmetric and skew-symmetric. So,

Le = De + We and Lp = Dp + Wp. We also make another kinematic assumption that the plastic

stretch Dp can be additively decomposed into Dp = Ds + Wvp, where DS represents the inelastic

strain rate due to oxidation-induced shrinkage and Wvp represents the inelastic strain rate due to

bulk viscoplastic deformation in the polymer. The readers are suggested to refer to earlier literature

for the detailed equations of the kinematics of deformation of polymers, such as [2–4]. Among the

set of reactions in elastomers, chain scission reaction plays a major role in oxidative shrinkage [2].

Hence, we have considered the following simple form for DS as Eq.(24), where γ < 0 is a material

96

parameter determining the volume shrinkage amount.

DS = γ.rsci
ox

· I,

with, tr(DS) = 3γrsci
ox .

(7.24)

Motivated by the free energy representation in [71], a separable form of the free energy is

considered as:

ψR(C e, Cox, psci, pref ) = ψmech

R + ψchem

R + ψdif f
R .

(7.25)

In this article, the 8-chain model of Arruda-Boyce [5] is utilized as a free energy function for

the mechanical part as,

R = nc · KB · T · N ·
ψmech

(cid:20)√

N −1Λζ + ln(

ζ
sinh(ζ)

(cid:21)

)

+

K
2

ln(J)2,

(7.26)

where nc is the chain density in the elastomer, KB denotes for Boltzmann’s constant, N is the

number of segments of length l in a chain, Λ is the chain stretch, and ζ abbreviates the inverse

Langevin function. In this article, we aim to describe the influence of heterogeneous oxidation

on the mechanical properties of elastomers. To do so, it is assumed that the chain density can be

affected by both oxidative reactions, chain scission, and network reformation. The chain density in

elastomers is reduced by chain scission reactions and is increased by network reformation reactions.

Therefore, the following relationship is proposed for considering the effect of oxidation on chain

density during aging.

(cid:0)

nc = nc0

1 − psci + Γpref

(cid:1)

.

(7.27)

In the equation. (27), nc0 represents the original chain density at the beginning of the exposure.

Both aging parameters can evolve inhomogeneously, thereby leading to inhomogeneous stress. Fol-

lowing the work presented in [47], we consider that the material’s bulk modulus does not depend

on oxidation. For the reactive part of the free energy, a simple quadratic function of the extent of

97

the reaction is considered.

ψchem

R =

(cid:2)

1
2

H sci(1 − psci)2 + H ref (1 − pref )2

(cid:3)

,

(7.28)

where the parameters H sci and H ref are the chemistry modulus for the chain scission and network

reformation reactions. Here we have used Flory-Huggins theory to define the free energy for oxy-

gen diffusion [42]. Oxygen is considered the only diffusing species into elastomers, as given by,

ψdif f
R = µo2

0 C o2

R + RT C o2

R

(cid:18)

ln(

ΩC o2
R
1 + ΩC o2
R

) + χ(

1
1 + ΩC o2
R

(cid:19)

)

,

(7.29)

where µo2

0 is the reference chemical potential for oxygen, χ is the dimensionless Flory-Huggins

interaction parameter, and Ω is the volume of a mole of oxygen. Thus, the complete free energy

expression can now be written as,

ψR = nc · KB · T · N ·

(cid:20)√

N −1Λ.ζ + ln(

ζ
sinh(ζ)

(cid:21)

)

+

K
2

H sci(1 − P sci)2 + H ref (1 − P ref )2

(cid:3)

(cid:2)

1
2

+µo2

0 C o2

R +RT.C o2

R

ln(

ln(J)2+
(cid:18)

ΩC o2
R
1 + ΩC o2
R

) + χ(

1
1 + ΩC o2
R

(cid:19)

)

.

(7.30)

From the free energy and using the guidelines for thermodynamic restriction, it is possible to get the

specific constitutive equations for the second Piola stress Se, chemical potential of oxygen, µo2and

affinity of each reaction. The Second Piola-Kirchoff stress can be expressed in the following form:

Se = 2

¯S =

∂ψR
∂Ce = K · lnJ · C−e + J − 2
3.N − Λ2
nc · KB · θ
· I.
N − Λ2
3

(cid:21)

(cid:20)

3

(cid:20)

¯S − 1
3

(cid:21)

( ¯Ce : ¯S). ¯C−e

,

(7.31)

(7.32)

Chemical potential can be obtained by taking the partial derivative of the total free energy with

98

respect to the oxygen concentration as,

µo2 =

∂ψ
∂CR

(cid:18)

= µo2

0 + RT

ln(

Ω.C o2
R
1 + ΩC o2
R

) + (

1
1 + ΩC o2
R

) + χ(

1
1 + ΩC o2
R

(cid:19)

)2

.

(7.33)

This work considers the effect of oxygen concentration and the extent of reactions in the reaction

rate functions for the coupled stress-oxidation scenario. Therefore, the following relationships are

proposed to control the influences of oxygen availability, aging parameters, and applied stress on

the evolution of reactions in the elastomers in the case of stress-coupled oxidation.

rsci
ox =

rref
ox =

αC
1 + βC
αC
1 + βC

rsci
0 exp

rref
0

exp

(cid:19)

(cid:18)

(cid:18)

−Esci
r
RT
−Eref
r
RT

(1 − psci).F sci,
(cid:19)

(1 − pref ) · F ref .

(7.34)

(7.35)

The F sci and F ref are reactive forces for chain scission and network reformation reactions,

which express the dependency of reaction rate on the mechanical load. We can express the reactive

forces for the individual reaction as,

F sci = Asci + γtr(C e · Se) − µo2,

F ref = Aref + µo2.

(7.36)

(7.37)

The Asci and Aref are defined as the affinity of chain scission and network reformation reactions.

Due to the fact that both oxidation reaction effects are considered in chain density during aging, the

following relationships can be obtained from the total free energy:

99

∂psci + H sci(1 − psci)
(cid:18)
N −1Λζ + ln

(cid:20)√

Asci = − ∂ψmech

R

= nc0KBN
Aref = − ∂ψmech

R

∂pref + H ref (1 − pref )
(cid:20)√
N −1Λζ + ln

= −Γnc0KBN

(cid:19)(cid:21)

ζ
sinh(ζ)

+ H sci(1 − psci),

(7.38)

(cid:19)(cid:21)

(cid:18)

ζ
sinh(ζ)

+ H ref (1 − pref ).

(7.39)

7.5.1: Finite element implementation

In the case of stress-coupled oxidation, the balance of momentum’s partial differential equation is

coupled with the balance of oxygen concentration equation. There are two governing differential

equations required to be solved in this case:

DivT + b = 0,

˙C o2

R = Div(D∇ · C o2

R ) − rsci

ox + rref
ox .

(7.40)

(7.41)

With b a non-inertial body force and Cauchy stress T . For solving the set of coupled equations,

we need boundaries and initial conditions to complete the theory. For the local force balance, a

pair of boundary conditions are considered, in which the displacement u and the surface traction

tensor are specified on sub-surfaces of the boundary ∂B of the body B in a sense ∂B = su ∪ st and

su ∩ st = ∅. Similarly, a pair of boundary conditions are considered for the oxygen concentration

equation, such that oxygen concentration and oxygen flux are known on the sub-surfaces of the

boundary. The initial condition data is taken as follows:

C o2

R (X, 0) = C o2

R (X).

(7.42)

U (X, 0) = U0(X),

100

Thus, the coupled set of Eqn.7.41, together with the boundary and initial condition, pose an

initial boundary value problem to be solved for the displacement u(X, t) and oxygen concentration

C o2 simultaneously. The finite element method is used for solving the coupled set of equations, in

which displacement and oxygen concentration are approximated in each element by:

X

X

u =

C o2

R =

uAN A,

(Co2

R )AN A,

(7.43)

(7.44)

with the index A = 1, 2, 3, · · ·, M denoting the nodes of the element, u, and C o2

R denoting nodal

displacements, and oxygen concentration, and N A the shape functions. We employ a standard

Galerkin approach. In the absence of body forces, the corresponding weak forms of the coupled

equation are:

Z

B

Z

T.

∂N A
∂x

Z

dv =

N A.(T.n)dA,

∂B
Z

˙Co2

R .N A.dv = −
Z

B

Div(D.∇.C o2
R ).
Z

B

∂N A
∂x

.dv +

Z

∂B

N A.(D.∇.Co2

R ).n.dv

(7.45)

(7.46)

−

N A.rox,sci.dv +

N A.rox,ref .dv.

B

B

This system of coupled equations is solved using a Newton procedure by defining the following

element-level residuals for the displacement and oxygen concentration,

Ru = −

Z

B

Z

T.

∂N A
∂x

Z

dv +

N A.(T.n)dA,

∂B
Z

Rc = −
Z

˙C o2

R .N A.dv −
Z

B

B

Div(D.∇.C o2

R ).

−

N A.rox,sci.dv +

N A.rox,ref .dv.

B

B

∂N A
∂x

.dv +

Z

∂B

N A.(D.∇.C o2

R ).n.dv

(7.47)

(7.48)

101

In addition to the residuals, the following tangents are also required for the iterative Newton solver:

(Kuu)AB = − ∂(Ru)A
∂uB ,
(Kuc)AB = − ∂(Ru)A
∂cB ,
(Kcu)AB = − ∂(Rc)A
∂uB ,
(Kcc)AB = − ∂(Rc)A
∂cB .

(7.49)

(7.50)

(7.51)

(7.52)

The residuals and elemental level stiffness matrix integrals are evaluated numerically using

Gaussian quadrature. Since this is a standard method in the finite element literature, we do not

present details here. In this work, the set of equations is solved numerically for each element by

writing a user-element subroutine in FORTRAN language programming. In its notation, ABAQUS

/ standard (2020) requires certain matrices denoted as RHS and AMATRX to be evaluated and/or

updated by the user element subroutine UEL. The work presented by Chester et.al [21] is referred

to for more details on the solution method.

7.6: Concluding remark

The work presented aims at providing a comprehensive view of the heterogeneous oxidation of elas-

tomers. Knowledge about oxygen uptake and penetration is essential to understand the mechanical

property changes of the elastomer during aging. Heterogeneous oxidation is the simultaneous pro-

cess of oxygen uptake and diffusion into elastomers. The introduced equations have described

the absorption of oxygen phenomena. The results of the modified model are in good agreement

with the experimentally reported data. Oxygen absorption has been considered as a function of

the aging temperature, the progress of aging, and the surface-volume ratio. The diffusion-reaction

equation has been extended in 3 dimensions to simulate the oxidation of complex-shaped elastomer

components. The ADI numerical solution is utilized for the sake of illustration of heterogeneous

oxidation over the thickness of the material. For better correlation with experimental data of the

degree of oxidation, some reaction rate functions have been proposed for chain scission and net-

102

work reformation reactions. The capability of functions for predicting the extent of reactions in

elastomers has been investigated by comparing the results with experimental data. In addition, this

article presents a modified constitutive model for chemo-mechanically coupled large deformation

of polymers. The model incorporates the coupled effect of the diffusion chemical reactions and

mechanical loading of the continuum-level constitutive response of the material. The model’s ca-

pability has been analyzed by implementing it into the Abaqus finite element software using the

user element subroutine. The model can predict the oxidation rates in polymers under DLO and

general stress-coupled oxidation scenarios.

103

CHAPTER 8: GRAPHIC USER

INTERFACE

8.1: Introduction

The integration of graphical user interfaces (GUIs) into finite element workflows significantly en-

hances user accessibility and efficiency, particularly for complex simulations. Abaqus, a powerful

finite element analysis (FEA) software, provides a Python scripting interface that enables advanced

customization and automation. Tkinter, the standard GUI library for Python, can be seamlessly used

within the Abaqus environment to create intuitive interfaces for tasks such as model input, material

parameter fitting, and postprocessing. This section outlines the core knowledge, procedural steps,

and best practices for integrating Tkinter-based GUIs into Abaqus.

This chapter describes the rationale, design, architecture and implementation of the GUI, with

a focus on how it integrates with the existing software ecosystem to enhance usability and function-

ality.

8.2: Motivation and Role of the GUI

8.2.1: Limitations of the Commercial Software

While the commercial software provides robust simulation capabilities (e.g., thermomechanical

solvers, structural analysis, transient simulations), it lacks support for:

• Advanced data preprocessing or augmentation workflows

• Customized model reduction techniques or post-processing routines

• User-developed features tailored to novel use cases or academic research

• Automation and integration across multiple simulation steps

104

8.2.2: Importance of Using a GUI in FEM Simulation

• User-Friendly Interaction A graphical user interface (GUI) provides an intuitive way for

users, especially those who are not programming experts, to input complex model parame-

ters, select materials, define boundary conditions, and configure simulation settings without

writing code.

• Error reduction GUIs can incorporate input validation, drop-down menus, and predefined

templates to reduce user errors in setting up simulations, leading to more reliable and consis-

tent results.

• Visual feedback and real-time updates With a GUI, users can visualize geometry, mesh

quality, boundary conditions, and intermediate results interactively. This immediate feedback

helps in detecting and correcting setup issues early.

• Streamlined Workflow GUIs integrate multiple steps of the simulation process, pre-processing,

solving, and post-processing, in a single environment, improving productivity and reducing

the learning curve.

• Accessibility and Collaboration By abstracting technical complexities, GUIs enable multi-

disciplinary teams (engineers, designers, analysts) to work collaboratively, sharing simula-

tion setups and results easily.

• Customization and Automation Modern FEM GUIs often allow scripting or batch pro-

cessing within a user-friendly interface, balancing ease of use with advanced control and

automation.

• Improved Decision Making With interactive plots, graphs, and visualization embedded in

GUIs, users can better interpret simulation outputs and make informed design or analysis

decisions.

8.2.3: Purpose of the Custom GUI

The custom-developed GUI serves the following purposes:

• Improved Accessibility Graphical User Interfaces (GUIs) eliminate the need for users to

105

Figure 8.1: Overview of the UMAT Workflow.

interact with complex command line instructions or script syntax. This makes it easier for

non-programmers to run simulations using intuitive inputs such as forms, drop-down menus,

and buttons, rather than modifying code manually.

• Enhanced Workflow Efficiency GUIs simplify multistep or repetitive tasks by integrating

them into single-button operations. This not only speeds up workflow, but also promotes

consistency and reduces the likelihood of user error, especially valuable in industrial or team-

based environments.

• Improved User ExperienceBy offering visual cues, progress indicators, and clear error mes-

sages, GUIs create a more engaging and user-friendly environment. This structured feedback

is especially helpful during lengthy or computation-heavy simulations, minimizing confusion

and frustration.

8.3: UMAT Derivation

8.3.1: UMAT Subroutine Workflow

This section illustrates the workflow for running a UMAT subroutine in Abaqus without utilizing the

dedicated graphical user interface (GUI), emphasizing the potential for errors and the complexity

often faced by first-time users. The complete process, from initial setup to execution, is depicted

in Fig.8.1. Procedures of how to run an Abaqus subroutine without using the GUI is listed in the

following:

106

1. Step 1: Adjust the variables in .for file (Fortran code)(see Fig.8.2);

2. Step 2: Open Abaqus and Select Job in Module list (see Fig.8.3);

3. Step 3: Choose Input File from Source drop-down list (see Fig.8.4);

4. Step 4: Select Input file (see Fig.8.5);

5. Step 5: In job editor, use General, upload .for file under User Subroutine file, select .ODB

in Result format (see Fig.8.6).

Figure 8.2: Variables in .for file.

Figure 8.3: Select Job in Module list.

107

Figure 8.4: Choose Input file from Source drop-down list.

Figure 8.5: Select Input file.

Figure 8.6: In job editor, use General, upload .for file under User subroutine file, select ODB
in Result format.

108

8.4: Objective of Integration

• Provide a custom user interface for input setup, including geometry, boundary conditions,

and material data.

• Automate the generation of input files (.inp) for Abaqus simulations.

• Incorporate custom pre-/post-processing routines (e.g., ML-driven preprocessing).

• To trigger and monitor Abaqus simulations from a single interface.

• To parse and visualize output (ODB) results post-simulation.

8.5: Background

The Abaqus Scripting Interface is a powerful tool that enables users to programmatically interact

with the model database ( mdb ), simulation jobs and result files. However, it traditionally requires

knowledge of Python scripting and familiarity with the Abaqus data structures, which may not

always be accessible to all users.

Tkinter, being a lightweight and widely available GUI package, can be embedded within Abaqus

Python scripts. This combination allows users to interact with Abaqus through buttons, text inputs,

and file selectors, thereby reducing the need for direct coding interaction.

8.6: Methodology

The methodology involves embedding a Tkinter-based GUI within a Python script that utilizes the

Abaqus Scripting Interface. The general flow of the integration process is as follows:

Tkinter GUI Initialization: A simple Tkinter window is created with user input fields and

buttons.

1. User Input Collection: Users provide the necessary simulation parameters through the GUI,

such as material properties, boundary conditions, or job settings.

2. Simulation Setup via Abaqus Scripting: Upon user interaction (e.g. clicking a button), the

script reads the input values and uses the Abaqus Scripting Interface to generate the model,

109

apply loads, mesh, and prepare the job.

3. Execution and Feedback: The script either submits the job automatically or allows the user

to manually trigger the execution. Feedback or status messages can be displayed within the

GUI or console.

8.7: Running the GUI Script in Abaqus

To streamline the simulation workflow within the Abaqus environment, a custom GUI was devel-

oped using Python’s tkinter library. This user interface offers a visual alternative to command-

line job submission, improving accessibility for users without scripting experience. The GUI inte-

grates model management, job configuration, and boundary condition setup, all from a centralized

interface.

The derived user interface (see Fig.8.7) allows users to specify or browse for an input file

(.inp), assign a job name, and choose or create a working model. The model selection dropdown

dynamically lists all available models in the current Abaqus database, while the optional model

creation entry enables the user to define a new model interactively.

In addition to job submission, the GUI supports the creation of basic static steps and the appli-

cation of default boundary conditions. Users can enable these features by filling in the Step Name

field or toggling the ”Add default boundary condition” checkbox. This boundary condition is

applied to a default region named ALL, if it exists within the model’s part definitions.

Upon submitting a job, the GUI invokes the Abaqus mdb.Job() constructor, linking the pro-

vided input file to the selected or newly created model. The job is then submitted and monitored

through Abaqus’ internal job manager. Completion status and any issues are reported via message

boxes and console logging, allowing users to verify successful execution directly from the GUI.

This extensible tool not only enhances productivity but also serves as a customizable framework

for automating simulation setup and execution workflows. With minor extensions, it can be adapted

to include additional step types, material definitions, and load conditions - all accessible through

a GUI framework consistent with Abaqus’ internal data structures. The modular structure of the

110

code also allows for easy extension and customization based on project needs. The script can be

executed within Abaqus by running:

a b a q u s c a e

s c r i p t = a b a q u s _ g u i . py

The full source code is presented in the Appendix, and the layout of the GUI is visualized in the

accompanying figure.

Figure 8.7: Tkinter-based GUI for Abaqus Job Submission.

Fig.8.7 shows the extended version of the Abaqus job submission GUI, which incorporates

interactive elements for model-driven simulations. This graphic user interface introduces three

major enhancements:

• Model Creation: Users can specify a new model name to be created on the fly. If the name

does not exist in the current session, a new model object will be added to the Abaqus model

database (mdb.models).

• Step Definition: A dedicated field allows users to enter a static step name, which will be

added to the selected model. This enables users to configure basic analysis sequences directly

from the GUI, avoiding manual scripting.

• Boundary Condition Control: An optional checkbox allows users to automatically assign

a fixed displacement boundary condition to a predefined set named ALL. This is useful for

running basic simulations without further pre-processing in the CAE interface.

111

The GUI reflects a modular design where each input field corresponds to a standard Abaqus

scripting operation. Once all fields are populated and the user clicks Run Job, the script processes

the inputs in sequence:

it verifies file paths, creates or selects the model, defines the step and

boundary condition (if applicable), then submits and monitors the job. The GUI provides immediate

feedback via pop-up message boxes and printed console messages.

This tool is particularly valuable for iterative simulation tasks or instructional environments,

where students or engineers can perform simulation runs without needing to understand the under-

lying scripting syntax. It can also be expanded to support material assignments, meshing operations,

or post-processing routines in future iterations.

8.8: Implementation Workflow:

Integrating Tkinter with Abaqus Scripting To successfully implement a graphical user interface

(GUI) using Tkinter within the Abaqus environment, a structured process must be followed. This

section outlines the detailed workflow for integrating Tkinter into Abaqus scripts and describes

how to run and use the resulting GUI.

8.8.1: Workflow Overview

The overall workflow for creating and executing a Tkinter-Abaqus integrated application consists

of the following steps:

112

Step Description

1

2

3

Design the GUI layout using Tkinter widgets (e.g., labels, entry fields,
buttons) to capture input parameters such as material properties, geom-
etry specifications, and analysis settings.
Define callback functions that are triggered by GUI interactions (e.g.,
clicking a ”Run” or ”Submit” button). These functions will gather inputs
from the GUI and use them to execute Abaqus commands.
Import Abaqus scripting modules (from abaqus import *, from
abaqusConstants import *) to enable programmatic control of mod-
els, parts, materials, and jobs.

4 Within the callback functions, use the Abaqus scripting interface to:

• Create or access an existing model (mdb.Model).
• Define geometry, material, and section properties.
• Assign boundary conditions, loads, and mesh controls.
• Generate input files and submit analysis jobs.

5

6

7

8

Add user feedback elements in the GUI (e.g., status messages, message
boxes, progress bars) to guide and inform users during model setup and
simulation execution.
Save the script (e.g., GUI_name.py). Test it inside Abaqus by:
• Running it through Abaqus/CAE: File > Run Script.
• Or

(batch mode):

abaqus cae

terminal

from the
noGUI=GUI_name.py.

Validate the script by ensuring the GUI inputs correctly translate into
Abaqus model features, the job submission runs without error, and the
output database (.odb) can be postprocessed as needed.
Optionally extend the GUI to include:

• File dialog access for user-defined data or geometry files.
• Postprocessing

read ODB results

routines

(e.g.,

via

session.openOdb).

• Export options for simulation results (e.g., CSV, plots).

Table 8.1: Workflow for Integrating Tkinter with Abaqus Scripting.

8.8.2: Detailed Steps

Step 1: Designing the GUI Layout

• Use Tkinter classes such as Tk(), Label, Entry, Button, and Frame to construct the graphical

layout.

113

• Arrange the widgets logically using .pack(), .grid(), or .place() layout managers.

r o o t = t k . Tk ( )

t k . L a b e l ( r o o t ,

t e x t =” E n t e r Young ' s Modulus : ” ) . p a c k ( )

e n t r y _ E = t k . E n t r y ( r o o t )

e n t r y _ E . p a c k ( )

Step 2: Defining Callback Functions

• Define Python functions that will be triggered by user actions (e.g., clicking a button).

• These functions should extract user input from the GUI and use it to drive the generation or

manipulation of Abaqus models.

Step 3: Embedding Abaqus Commands

• Abaqus scripting commands such as mdb.Model, mdb.Job, or session.openOdb can be em-

bedded directly inside the callback functions.

• This allows user inputs from the GUI to dynamically control the simulation parameters.

Step 4:

Initializing and Running the Tkinter Main Loop

• After defining the GUI elements and callback functions, the Tkinter event loop must be

started by calling root.mainloop().

• This command will keep the window open and responsive to user inputs.

r o o t . m a i n l o o p ( )

Step 5: Executing the Script Inside Abaqus

• The script must be run inside the Abaqus Python interpreter to have access to both Tkinter

and Abaqus modules (abaqus, abaqusConstants, etc.).

• There are two common ways to run the script:

1. Inside Abaqus/CAE:

– Open Abaqus/CAE.

– Go to File → Run Script, then select and execute the Python file.

114

2. From the Command Line:

– Use the following command to run the script in no-GUI (batch) mode:

a b a q u s c a e noGUI=GUI_name . py

3. Note: If complex GUI interactions are required, running within Abaqus/CAE is gener-

ally preferred.

8.8.3: Visual Summary of the Workflow

The flowchart (Fig.8.8) illustrates the step-by-step process of integrating a Tkinter-based graphical

user interface (GUI) with Abaqus scripting functionalities. Each stage of the workflow plays a cru-

cial role in enabling user-friendly model creation and job submission within Abaqus. The process

is as follows.

1. Design GUI Layout The first step involves creating the user interface using Tkinter widgets

such as labels, entry fields, buttons, and frames. This layout defines how users will input

parameters, interact with the script, and visualize the outputs.

2. Define Callback Functions Next, callback functions are programmed to respond to user inter-

actions (e.g., button clicks). These functions act as bridges between the GUI and the Abaqus

Python scripting environment.

3. Abaqus commands embedded within the callback functions are inserted in Abaqus scripting

commands (such as part creation, material definition, meshing, and job submission). This

ensures that GUI actions directly trigger model-building and analysis operations in Abaqus.

4. Run Tkinter Main Loop After setting up the GUI and functions, the Tkinter Main Event

Loop (mainloop()) is executed. This keeps the GUI active, allowing users to interact with it

dynamically.

5. Execute Script Inside Abaqus Finally, the entire Python script is run inside the Abaqus/-

CAE or Abaqus/Viewer environment. This ensures that all necessary Abaqus modules and

libraries are properly loaded and that the GUI can control Abaqus operations seamlessly.

115

This structured approach enables users to build custom, accessible interfaces for complex sim-

ulation workflows, thereby improving efficiency and reducing the likelihood of user errors.

Figure 8.8: Flowchart of Designing Process.

8.9: Designed User Interface

This section presents the designed GUI and its core functionalities. The primary purpose of the

GUI is to streamline operations by bypassing predefined, lengthy subroutines that are prone to user

error.

Figure 8.9: Graphical User Interface (GUI) Layout.

In addition to Fig.8.9, a simplified version (shown in Fig.8.10) is also developed. This sim-

plified interface is designed specifically for executing user subroutines, without requiring manual

selection or interaction with the full simulation setup. This version bypasses the UMAT selection

step described in the previous section, enabling quicker deployment and testing of subroutine be-

havior in a controlled environment. The original selection interface is intended to provide detailed

explanations of each environmental aging concept along with references to the corresponding orig-

inal research papers.

116

DesignGUI LayoutDefineCallbackFunctionsEmbedAbaqusCommandRunTkinterMain LoopExecuteScriptInsideAbaqusEndMaterial ModelSelectionInput Option 1:Copy, paste Exp. data from.txt fileInput Option 2:Load .csv dataPrint Experimental DataRun SubroutineRun Matlab CodeFigure 8.10: Layout of Graphical User Interface for Subroutine Execution Only.

Fig.8.11 shows the interface for selecting a sample material model. It is a click-based selec-

tion tool that navigates the user to the interface on the right-hand side. Once the user inputs the

experimental stress-strain data into the designated fields and clicks the Material Model button, the

tool automatically runs a Genetic Algorithm to identify the best-fit material parameters. These

parameters can then be used for further correction or subsequent calculations.

Figure 8.11: Material Model Selection Tool and Parameter Fitting Interface.

117

The implemented GUI provides two user-friendly methods for inputting stress-strain data: di-

rect copy-paste from a .txt file or loading from a .csv file. As shown in Fig. 8.9, users can paste

space-separated stress-strain pairs into a multi-line text box or click the Load CSV File button to im-

port a file with two columns—stress and strain. In both cases, clicking the Plot Data button parses

the input and generates a stress-strain plot, enabling quick visual validation before proceeding to

model fitting or optimization. The code is provided in Appendix.

Fig. 8.12 illustrates the interface design options. The first option is a simple selection button,

while the second is a drop-down menu that provides a list of predefined choices.

Figure 8.12: Selection options available in the Genetic Algorithm user interface.

The Graphic User Interface (GUI) of the Genetic Algorithm (GA) is illustrated in Fig.8.13.

The purpose of this function is to simplify the process of determining the parameters of the best-

fit material model. Instead of manually entering the material model code and its parameter range,

users only need to enter the experimental stress and strain data. Once the data are provided, users

can simply click the ”Run Genetic Algorithm” button to obtain the optimized parameters that best

fit the provided experimental results.

Figure 8.13: Genetic Algorithm User Interface.

118

Option 1: selectOption 2: drop-down8.10: Conclusion

Integrating Tkinter-based GUIs with Abaqus scripting offers an accessible and efficient solution

for automating simulation setup, reducing the learning curve for users less familiar with scripting.

This approach facilitates the development of intuitive interfaces for defining materials, configur-

ing models, submitting jobs, and managing results—thereby streamlining the overall workflow in

computational simulations. All Tkinter code implementations referenced in the previous sections

are provided in the Appendix.

119

BIBLIOGRAPHY

[1] S Abouzahr, GL Wilkes, and Z Ophir. Structure-property behaviour of segmented polyether-
mdi-butanediol based urethanes: effect of composition ratio. Polymer, 23(7):1077–1086,
1982.

[2] Nicoli M Ames, Vikas Srivastava, Shawn A Chester, and Lallit Anand. A thermo-
mechanically coupled theory for large deformations of amorphous polymers. part ii: Ap-
plications. International Journal of Plasticity, 25(8):1495–1539, 2009.

[3] L Anand and C Su. A theory for amorphous viscoplastic materials undergoing finite de-
formations, with application to metallic glasses. Journal of the Mechanics and Physics of
Solids, 53(6):1362–1396, 2005.

[4] Lallit Anand, Nicoli M Ames, Vikas Srivastava, and Shawn A Chester. A thermo-
mechanically coupled theory for large deformations of amorphous polymers. part i: For-
mulation. International Journal of Plasticity, 25(8):1474–1494, 2009.

[5] Ellen M Arruda and Mary C Boyce. A three-dimensional constitutive model for the large
stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids,
41(2):389–412, 1993.

[6] Roger A Assink, Mathew Celina, Julie M Skutnik, and Douglas J Harris. Use of a respirom-
eter to measure oxidation rates of polymeric materials at ambient temperatures. Polymer,
46(25):11648–11654, 2005.

[7] L Audouin, V Gueguen, A Tcharkhtchi, and J Verdu. “close loop” mechanistic schemes for
hydrocarbon polymer oxidation. Journal of Polymer Science Part A: Polymer Chemistry,
33(6):921–927, 1995.

[8] Ludmila Audouin, Laurence Achimsky, and Jacques Verdu. Kinetic modeling of low-
In Handbook of polymer degradation,

temperature oxidation of hydrocarbon polymers.
pages 753–764. CRC Press, 2000.

[9] Amir Bahrololoumi, Aref Ghaderi, Mamoon Shaafaey, and Roozbeh Dargazany. A micro-
mechanical constitutive model to predict hygrothermal aging of cross-linked polymers. In
ASME International Mechanical Engineering Congress and Exposition, volume 85680, page
V012T12A039. American Society of Mechanical Engineers, 2021.

[10] 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.

120

[11] Kenneth J Bowles and Gregory Nowak. Thermo-oxidative stability studies of celion
6000/pmr-15 unidirectional composites, pmr-15, and celion 6000 fiber. Journal of Com-
posite Materials, 22(10):966–985, 1988.

[12] P Budrugeac. Accelerated thermal ageing of nitrile-butadiene rubber under air pressure.

Polymer Degradation and stability, 47(1):129–132, 1995.

[13] F Bueche. Molecular basis for the mullins effect. Journal of Applied Polymer Science,

4(10):107–114, 1960.

[14] 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.

[15] M Celina, J Wise, DK Ottesen, KT Gillen, and RL Clough. Correlation of chemical and me-
chanical property changes during oxidative degradation of neoprene. Polymer degradation
and Stability, 68(2):171–184, 2000.

[16] 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.

[17] J Chen, W Wang, and J Li. Effect of thermal vacuum environment on the performance of

adhesive bonding materials, aerosp. Mater. Technol, 42(6):92–96, 2012.

[18] Jun Chen, Nengwen Ding, Zhifeng Li, and Wei Wang. Organic polymer materials in the

space environment. Progress in Aerospace Sciences, 83:37–56, 2016.

[19] Shawn A Chester and Lallit Anand. A coupled theory of fluid permeation and large deforma-
tions for elastomeric materials. Journal of the Mechanics and Physics of Solids, 58(11):1879–
1906, 2010.

[20] Shawn A Chester and Lallit Anand. A thermo-mechanically coupled theory for fluid per-
meation in elastomeric materials: application to thermally responsive gels. Journal of the
Mechanics and Physics of Solids, 59(10):1978–2006, 2011.

[21] Shawn A Chester, Claudio V Di Leo, and Lallit Anand. A finite element implementation of
a coupled diffusion-deformation theory for elastomeric gels. International Journal of Solids
and Structures, 52:1–18, 2015.

[22] Sung-Seen Choi and Jong-Chul Kim. Lifetime prediction and thermal aging behaviors of sbr
and nbr composites using crosslink density changes. Journal of Industrial and Engineering
Chemistry, 18(3):1166–1170, 2012.

[23] Hsoung-Wei Chou, Jong-Shin Huang, and Shih-Ting Lin. Effects of thermal aging on fatigue

121

of carbon black–reinforced epdm rubber. Journal of applied polymer science, 103(2):1244–
1251, 2007.

[24] X Colin, L Audouin, and J Verdu. Kinetic modelling of the thermal oxidation of polyisoprene
elastomers. part 3: Oxidation induced changes of elastic properties. Polymer degradation
and stability, 92(5):906–914, 2007.

[25] X Colin, C Marais, and J Verdu. A new method for predicting the thermal oxidation of
thermoset matrices: application to an amine crosslinked epoxy. Polymer Testing, 20(7):795–
803, 2001.

[26] X Colin and J Verdu. Thermal ageing and lifetime prediction for organic matrix composites.

Plastics, rubber and composites, 32(8-9):349–356, 2003.

[27] X Colin and J Verdu. Strategy for studying thermal oxidation of organic matrix composites.

Composites Science and technology, 65(3-4):411–419, 2005.

[28] Xavier Colin and Jacques Verdu. Mechanisms and kinetics of organic matrix thermal oxi-
dation. In Long-term durability of polymeric matrix composites, pages 311–344. Springer,
2011.

[29] AV Cunliffe and A Davis. Photo-oxidation of thick polymer samples—part ii: The influ-
ence of oxygen diffusion on the natural and artificial weathering of polyolefins. Polymer
Degradation and Stability, 4(1):17–37, 1982.

[30] Hüsnü Dal, Osman Gültekin, and Kemal Açıkgöz. An extended eight-chain model for hy-
perelastic and finite viscoelastic response of rubberlike materials: Theory, experiments and
numerical aspects. Journal of the Mechanics and Physics of Solids, 145:104159, 2020.

[31] Roozbeh Dargazany, Vu Ngoc Khiêm, Uwe Navrath, and Mikhail Itskov. Network evolu-
tion model of anisotropic stress softening in filled rubber-like materials: Parameter identifi-
cation and finite element implementation. Journal of Mechanics of Materials and Structures,
7(8):861–885, 2013.

[32] J Decelle, N Huet, and V Bellenger. Oxidation induced shrinkage for thermally aged epoxy

networks. Polymer Degradation and Stability, 81(2):239–248, 2003.

[33] Mohammad Dehghany, Haohui Zhang, Reza Naghdabadi, and Yuhang Hu.

A
thermodynamically-consistent large deformation theory coupling photochemical reaction
and electrochemistry for light-responsive gels. Journal of the Mechanics and Physics of
Solids, 116:239–266, 2018.

[34] Julie Diani, Bruno Fayolle, and Pierre Gilormini. A review on the mullins effect. European

Polymer Journal, 45(3):601–612, 2009.

122

[35] Bradley K Diedrick. Effectsof prior aging at 260° c in argon on inelastic deformation behav-

ior of pmr-15 polymer at 260° c: Experiment and modeling. 2010.

[36] Jim Douglas. Alternating direction methods for three space variables. Numerische Mathe-

matik, 4:41–63, 1962.

[37] Fernando P Duda, Angela C Souza, and Eliot Fried. A theory for species migration in a
finitely strained solid with application to polymer network swelling. Journal of the Mechan-
ics and Physics of Solids, 58(4):515–529, 2010.

[38] AE 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(2):189–206, 2010.

[39] Esteve Ernault, Emmanuel Richaud, and Bruno Fayolle. Thermal-oxidation of epoxy/amine
followed by glass transition temperature changes. Polymer Degradation and Stability,
138:82–90, 2017.

[40] BCXALVJ Fayolle, X Colin, L Audouin, and J Verdu. Mechanism of degradation induced

embrittlement in polyethylene. Polymer degradation and stability, 92(2):231–238, 2007.

[41] Bruno Fayolle, Emmanuel Richaud, Jacques Verdu, and Fabienne Farcas. Embrittlement of
polypropylene fibre during thermal oxidation. Journal of materials science, 43:1026–1032,
2008.

[42] Paul J Flory. Principles of polymer chemistry. Cornell university press, 1953.

[43] 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.

[44] Aref Ghaderi, Vahid Morovati, Amir Bahrololoumi, and Roozbeh Dargazany. A physics-
informed neural network constitutive model for cross-linked polymers.
In ASME In-
ternational Mechanical Engineering Congress and Exposition, volume 84607, page
V012T12A007. American Society of Mechanical Engineers, 2020.

[45] 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.

[46] 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.

123

[47] Marco Gigliotti, Loic Olivier, Dinh Quy Vu, Jean-Claude Grandidier, and Marie Christine
Lafarie-Frenot. Local shrinkage and stress induced by thermo-oxidation in composite mate-
rials at high temperatures. Journal of the Mechanics and Physics of Solids, 59(3):696–712,
2011.

[48] Pieter Gijsman, Weifu Dong, Adam Quintana, and Mathew Celina. Influence of temperature
and stabilization on oxygen diffusion limited oxidation profiles of polyamide 6. Polymer
Degradation and Stability, 130:83–96, 2016.

[49] Kenneth T Gillen, Mathias Celina, and Michael R Keenan. Methods for predicting more con-
fident lifetimes of seals in air environments. Rubber chemistry and technology, 73(2):265–
283, 2000.

[50] Kenneth T Gillen and Roger L Clough. Rigorous experimental confirmation of a theoretical

model for diffusion-limited oxidation. Polymer, 33(20):4358–4365, 1992.

[51] Sanjay Govindjee and Juan C Simo. Coupled stress-diffusion: Case ii. Journal of the Me-

chanics and Physics of Solids, 41(5):863–887, 1993.

[52] Morton E Gurtin, Eliot Fried, and Lallit Anand. The mechanics and thermodynamics of

continua. Cambridge University Press, 2010.

[53] David E Hanson, Marilyn Hawley, Robert Houlton, Kiran Chitanvis, Philip Rae, E Bruce Or-
ler, and Debra A Wrobleski. Stress softening experiments in silica-filled polydimethylsilox-
ane provide insight into a mechanism for the mullins effect. Polymer, 46(24):10989–10995,
2005.

[54] Leila Haroonabadi, Ali Dashti, and Mahsa Najipour. Investigation of the effect of thermal
aging on rapid gas decompression (rgd) resistance of nitrile rubber. Polymer Testing, 67:37–
45, 2018.

[55] A Herzig, M Johlitz, and A Lion. Ageing experimental investigation on the consumption
of oxygen and its diffusion into elastomers during the process of ageing. In Constitutive
Models for Rubber IX, pages 39–44. CRC Press, 2015.

[56] A Herzig, M Johlitz, and A Lion. Consumption and diffusion of oxygen during the ther-
moxidative ageing process of elastomers: Absorption und diffusion von sauerstoff während
der thermooxidativen alterung von elastomeren. Materialwissenschaft und Werkstofftechnik,
47(5-6):376–387, 2016.

[57] A Herzig, M Johlitz, RJ Pazur, I Lopez-Carreon, and CG Porter. Thermal ageing of peroxide-
cured nbr, part ii: Diffusion-limited oxidation and constitutive modelling. In Constitutive
Models for Rubber XI, pages 587–592. CRC Press, 2019.

124

[58] A Herzig, L Sekerakova, M Johlitz, and A Lion. A modelling approach for the heterogeneous
oxidation of elastomers. Continuum Mechanics and Thermodynamics, 29:1149–1161, 2017.

[59] Alexander Herzig, Michael Johlitz, and Alexander Lion. An experimental set-up to analyse
the oxygen consumption of elastomers during ageing by using a differential oxygen analyser.
Continuum Mechanics and Thermodynamics, 27:1009–1017, 2015.

[60] E Hoffman, W Daugherty, E Skidmore, K Dunn, and D Fisher. Aging performance of model
9975 package fluoroelastomer o-rings. Technical report, Savannah River Site (SRS), Aiken,
SC (United States), 2011.

[61] Kurt Hornik, Maxwell Stinchcombe, and Halbert White. Multilayer feedforward networks

are universal approximators. Neural networks, 2(5):359–366, 1989.

[62] Rubber Houwink. Slipping of molecules during the deformation of reinforced rubber. Rub-

ber Chemistry and Technology, 29(3):888–893, 1956.

[63] Yani Ioannou, Duncan Robertson, Darko Zikic, Peter Kontschieder, Jamie Shotton, Matthew
Brown, and Antonio Criminisi. Decision forests, convolutional networks and the models in-
between. arXiv preprint arXiv:1603.01250, 2016.

[64] Michael Johlitz. On the representation of ageing phenomena. The Journal of Adhesion,

88(7):620–648, 2012.

[65] 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.

[66] Michael Johlitz and Alexander Lion. Chemo-thermomechanical ageing of elastomers based
on multiphase continuum mechanics. Continuum Mechanics and Thermodynamics, 25:605–
624, 2013.

[67] S Kamarudin, PY Le Gac, Y Marco, and AH Muhr. Formation of crust on natural rubber
after ageing. In Proceedings of the 7th European Conference on Constitutive Models for
Rubber, ECCMR, pages 197–202, 2011.

[68] Leif Kari. Dynamic stiffness of chemically and physically ageing rubber vibration isolators
in the audible frequency range: Part 1: constitutive equations. Continuum Mechanics and
Thermodynamics, 29(5):1027–1046, 2017.

[69] Yoichi Kodera and Benjamin J McCoy. Distribution kinetics of radical mechanisms: re-

versible polymer decomposition. AIChE journal, 43(12):3205–3214, 1997.

[70] Evgenia Kollia, Antonios Vavouliotis, and Vassilios Kostopoulos. Thermal ageing of carbon

125

fiber-reinforced cyanate ester composites under inert and oxidative environment. Polymer
Composites, 40(S2):E1388–E1396, 2019.

[71] 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, 137:103858, 2020.

[72] 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.

[73] Shin-ichi Kuroda, Koichiro Terauchi, Kyohei Nogami, and Itaru Mita. Degradation of aro-
matic polymers—i. rates of crosslinking and chain scission during thermal degradation of
several soluble aromatic polymers. European polymer journal, 25(1):1–7, 1989.

[74] Julie Lambert-Diani and Christian Rey. New phenomenological behavior laws for rubbers
and thermoplastic elastomers. European Journal of Mechanics-A/Solids, 18(6):1027–1043,
1999.

[75] Pierre Yves Le Gac, Mathew Celina, Gérard Roux, Jacques Verdu, Peter Davies, and Bruno
Fayolle. Predictive ageing of elastomers: Oxidation driven modulus changes for polychloro-
prene. Polymer Degradation and Stability, 130:348–355, 2016.

[76] Pierre-Yves Le Gac, V Le Saux, M Paris, and Yann Marco. Ageing mechanism and me-
chanical degradation behaviour of polychloroprene rubber in a marine environment: Com-
parison of accelerated ageing and long term exposure. Polymer degradation and stability,
97(3):288–296, 2012.

[77] Jichun Li and Yi-Tung Chen. Computational partial differential equations using MATLAB®.

Crc Press, 2019.

[78] 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.

[79] Kaspar Loeffel and Lallit Anand. A chemo-thermo-mechanically coupled theory for elastic–
viscoplastic deformation, diffusion, and volumetric swelling due to a chemical reaction. In-
ternational Journal of Plasticity, 27(9):1409–1431, 2011.

[80] Thanh-Tam Mai, Yoshihiro Morishita, and Kenji Urayama. Novel features of the mullins
effect in filled elastomers revealed by stretching measurements in various geometries. Soft
matter, 13(10):1966–1977, 2017.

[81] Gilles Marckmann, Erwan Verron, Laurent Gornet, Grégory Chagnon, Pierre Charrier, and
Patrice Fort. A theory of network alteration for the mullins effect. Journal of the Mechanics

126

and Physics of Solids, 50(9):2011–2028, 2002.

[82] M Minervino, M Gigliotti, MC Lafarie-Frenot, and JC Grandidier. The effect of thermo-
oxidation on the mechanical behaviour of polymer epoxy materials. Polymer testing,
32(6):1020–1028, 2013.

[83] M Minervino, M Gigliotti, MC Lafarie-Frenot, and JC Grandidier. A coupled experimen-
tal/numerical approach for the modelling of the local mechanical behaviour of epoxy poly-
mer materials. Journal of the Mechanics and Physics of Solids, 67:129–151, 2014.

[84] 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.

[85] Leonard Mullins. Softening of rubber by deformation. Rubber chemistry and technology,

42(1):339–362, 1969.

[86] 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.

[87] 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:369–385, 2020.

[88] 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.

[89] L Olivier, NQ Ho, JC Grandidier, and MC Lafarie-Frenot. Characterization by ultra-micro
indentation of an oxidized epoxy polymer: correlation with the predictions of a kinetic model
of oxidation. Polymer Degradation and Stability, 93(2):489–497, 2008.

[90] Ozgur Ozmen. Effects of prior aging at 316 c in argon on inelastic deformation behavior of

pmr-15 polymer at 316 c: Experiment and modeling. 2009.

[91] Abhishek Panigrahi, Abhishek Shetty, and Navin Goyal. Effect of activation functions on
the training of overparametrized neural nets. arXiv preprint arXiv:1908.05660, 2019.

[92] Marina Pecora, Yannick Pannier, Marie-Christine Lafarie-Frenot, Marco Gigliotti, and
Camille Guigon. Effect of thermo-oxidation on the failure properties of an epoxy resin.
Polymer Testing, 52:209–217, 2016.

127

[93] Vorapong Pimolsiriphol, Pongdhorn Saeoui, and Chakrit Sirisinha. Relationship among ther-
mal ageing degradation, dynamic properties, cure systems, and antioxidants in natural rubber
vulcanisates. Polymer-Plastics Technology and Engineering, 46(2):113–121, 2007.

[94] Kishore V Pochiraju. Modeling thermo-oxidative aging and degradation of composites. In
Long-term durability of polymeric matrix composites, pages 383–425. Springer, 2011.

[95] KV Pochiraju, GP Tandon, and GA Schoeppner. Evolution of stress and deformations in
high-temperature polymer matrix composites during thermo-oxidative aging. Mechanics of
Time-Dependent Materials, 12:45–68, 2008.

[96] Adam Quintana and Mathew C Celina. Overview of dlo modeling and approaches to predict
heterogeneous oxidative polymer degradation. Polymer Degradation and Stability, 149:173–
191, 2018.

[97] 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:1011–1023, 2015.

[98] Sudip Ray and Ralph P Cooney. Thermal degradation of polymer and polymer composites.
In Handbook of environmental degradation of materials, pages 185–206. Elsevier, 2018.

[99] LM Rincon-Rubio, B Fayolle, L Audouin, and J Verdu. A general solution of the closed-
loop kinetic scheme for the thermal oxidation of polypropylene. Polymer Degradation and
Stability, 74(1):177–188, 2001.

[100] Trisha Sain, Kaspar Loeffel, and Shawn Chester. A thermo–chemo–mechanically coupled
constitutive model for curing of glassy polymers. Journal of the Mechanics and Physics of
Solids, 116:267–289, 2018.

[101] John Scheirs, Stephen W Bigger, and Norman C Billingham. A review of oxygen uptake
techniques for measuring polyolefin oxidation. Polymer testing, 14(3):211–241, 1995.

[102] Lars Steinke. Ein Beitrag zur Simulation der thermo-oxidativen Alterung von Elastomeren.

VDI-Verlag, 2013.

[103] GP Tandon, KV Pochiraju, and GA Schoeppner. Modeling of oxidative development in

pmr-15 resin. Polymer Degradation and Stability, 91(8):1861–1869, 2006.

[104] 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.

[105] Namrata Singh Tomer, Florence Delor-Jestin, RP Singh, and Jacques Lacoste. Cross-linking

128

assessment after accelerated ageing of ethylene propylene diene monomer rubber. Polymer
Degradation and Stability, 92(3):457–463, 2007.

[106] Leslie Ronald George Treloar. The physics of rubber elasticity. Oxford University Press,

USA, 1975.

[107] GJ Van Amerongen.

Influence of structure of elastomers on their permeability to gases.

Journal of polymer science, 5(3):307–332, 1950.

[108] Arnaud Vieyres, Roberto Pérez-Aparicio, Pierre-Antoine Albouy, Olivier Sanseau, Kay
Saalwächter, Didier R Long, and Paul Sotta. Sulfur-cured natural rubber elastomer networks:
correlating cross-link density, chain orientation, and mechanical response by combined tech-
niques. Macromolecules, 46(3):889–899, 2013.

[109] Haixiao Wan, Ke Gao, Sai Li, Liqun Zhang, Xiaohui Wu, Xiaodong Wang, and Jun Liu.
Chemical bond scission and physical slippage in the mullins effect and fatigue behavior of
elastomers. Macromolecules, 52(11):4209–4221, 2019.

[110] Elias Wang, Atli Kosson, and Tong Mu. Deep action conditional neural network for frame

prediction in atari games. Technical report, Technical Report, Stanford University, 2017.

[111] J Wise, KT Gillen, and RL Clough. Quantitative model for the time development of diffusion-

limited oxidation profiles. Polymer, 38(8):1929–1944, 1997.

[112] Jiangtao Wu, Zeang Zhao, Craig M Hamel, Xiaoming Mu, Xiao Kuang, Zaoyang Guo, and
H Jerry Qi. Evolution of material properties during free radical photopolymerization. Jour-
nal of the Mechanics and Physics of Solids, 112:25–49, 2018.

[113] Chen Xuan and Lihua Jin. Concurrent reaction and diffusion in photo-responsive hydrogels.

Journal of the Mechanics and Physics of Solids, 124:599–611, 2019.

[114] YS Rohana Yahya, AR Azura, and Z Ahmad. Effect of curing systems on thermal degra-
dation behaviour of natural rubber (smr cv 60). Journal of Physical Science, 22(2):1–14,
2011.

[115] Yang. Paper-related code release, 2025. https://github.com/chyg-code/appendix/tree/main.

[116] C Joshua Yeh, Matt Dowland, Randall G Schmidt, and Kenneth R Shull. Fracture and ther-
mal aging of resin-filled silicone elastomers. Journal of Polymer Science Part B: Polymer
Physics, 54(2):263–273, 2016.

[117] Oon H Yeoh. Characterization of elastic properties of carbon-black-filled rubber vulcan-

izates. Rubber chemistry and technology, 63(5):792–805, 1990.

129

[118] Oon H Yeoh and PD Fleming. A new attempt to reconcile the statistical and phenomeno-
logical theories of rubber elasticity. Journal of Polymer Science Part B: Polymer Physics,
35(12):1919–1931, 1997.

[119] 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.

[120] Danming Zhong, Yuhai Xiang, Tenghao Yin, Honghui Yu, Shaoxing Qu, and Wei Yang.
A physically-based damage model for soft elastomeric materials with anisotropic mullins
effect. International Journal of Solids and Structures, 176:121–134, 2019.

130

APPENDIX

All Graphic-User Interface source code has been uploaded and can be viewed at [115]. For

integration point of micro-sphere approach, the integration points and weighting factors of the unit-

sphere is shown in Table 8.2

i
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

di(1)
0.0
0.0
1.0
0.0
0.0
0.707106781187
0.707106781187
0.707106781187
0.707106781187
0.836095596749
0.836095596749
0.836095596749
0.836095596749
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.387907304067

di(2)
0.0
1.0
0.0
0.707106781187
0.707106781187
0.0
0.0
0.707106781187
0.707106781187
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.836095596749
0.836095596749
0.836095596749
0.836095596749
0.387907304067
0.387907304067
0.387907304067
-0.387907304067

di(3)
1.0
0.0
0.0
0.707106781187
0.707106781187
0.707106781187
0.707106781187
0.0
0.0
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.387907304067
0.836095596749
0.836095596749
0.836095596749
0.836095596749

wi
0.0265214244093
0.0265214244093
0.0265214244093
0.0199301476312
0.0199301476312
0.0199301476312
0.0199301476312
0.0199301476312
0.0199301476312
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487
0.0250712367487

Table 8.2: Integration points and weighting factors of the unit-sphere

131