Part 1: Mathematical model, analysis, and analytical solutions

Authored by

C. Xu
Graduate Student, University of Houston

M. K. Mudunuru
Graduate Student, University of Houston

K. B. Nakshatrala
Department of Civil & Environmental Engineering
University of Houston, Houston, Texas 77204–4003
phone: +1-713-743-4418, e-mail: knakshatrala@uh.edu
website: http://www.cive.uh.edu/faculty/nakshatrala
3

D/D0

arXiv:1511.05538v2 [cond-mat.mtrl-sci] 21 Jan 2016

Material degradation due to moisture and temperature

2

1
0

Experimental data for tension
Proposed model for tension
Experimental data for compression
Proposed model for compression
Experimental data for shear
Proposed model for shear

0.001
0.002
Trace of strain

0.003

This figure shows a good agreement between the experimental data and the proposed constitutive model for the diffusivity under tension, compression and shear. Analysis is performed
for various extracted sample sizes, and the coefficient of determination is found to be close
to 1. This calibration study provides confidence in the proposed model to apply for brittle
materials.

2016
Computational & Applied Mechanics Laboratory

Material degradation due to moisture and temperature
Part 1: Mathematical model, analysis, and analytical solutions

C. Xu, M. K. Mudunuru, and K. B. Nakshatrala
Department of Civil and Environmental Engineering, University of Houston.
Correspondence to: knakshatrala@uh.edu
Abstract. The mechanical response, serviceability, and load bearing capacity of materials and
structural components can be adversely affected due to external stimuli, which include exposure to
a corrosive chemical species, high temperatures, temperature fluctuations (i.e., freezing-thawing),
cyclic mechanical loading, just to name a few. It is, therefore, of paramount importance in several
branches of engineering – ranging from aerospace engineering, civil engineering to biomedical engineering – to have a fundamental understanding of degradation of materials, as the materials in
these applications are often subjected to adverse environments. As a result of recent advancements
in material science, new materials like fiber-reinforced polymers and multi-functional materials that
exhibit high ductility have been developed and widely used; for example, as infrastructural materials or in medical devices (e.g., stents). The traditional small-strain approaches of modeling these
materials will not be adequate. In this paper, we study degradation of materials due to an exposure to chemical species and temperature under large-strain and large-deformations. In the first
part of our research work, we present a consistent mathematical model with firm thermodynamic
underpinning. We then obtain semi-analytical solutions of several canonical problems to illustrate
the nature of the quasi-static and unsteady behaviors of degrading hyperelastic solids.

NOMENCLATURE
ρ

A
ζ
ψ
λ, µ
κ
u
v
ϑ
c
Rs
cp

density of solid in deformed conﬁguration [kg m−3 ]
speciﬁc Helmholtz potential [J kg−1 ]
dissipation functional [J kg−1 s−1 ]
strain energy density functional [J m−3 ]
Lamé parameters [Pa]
bulk modulus [Pa]
displacement [m]
velocity [m s−1 ]
temperature [K]
concentration [1]
speciﬁc vapor constant [J kg−1 K−1 ]
heat capacity [J kg−1 K−1 ]

Key words and phrases. Degradation; aging; continuum damage mechanics; coupled chemo-thermo-mechano
analysis; semi-analytical solutions; constitutive modeling; hyperelasticity.
1

MϑE
McE
dϑc
κ
η
Dϑϑ
Dκκ
Dϑκ , Dκϑ
T
h
q
h
q

thermal expansion tensor [J m−3 K−1 ]
chemical expansion tensor [J m−3 ]
thermo-chemo coupled parameter [J kg−1 K−1 ]
speciﬁc chemical potential [J kg−1 ]
speciﬁc entropy [J kg−1 K−1 ]
thermal diﬀusion tensor [m2 s−1 ]
diﬀusivity tensor [m2 s−1 ]
Dufour-Soret eﬀect tensors [m2 s−1 ]
Cauchy stress [Pa]
diﬀusive ﬂux vector [kg m−2 s−1 ]
heat ﬂux vector [J m−2 s−1 ]
volumetric source [kg m−3 s−1 ]
volumetric heat source [J m−3 s−1 ]

1. INTRODUCTION AND MOTIVATION
aterial and structural degradation is a major problem in infrastructure and various other
real-life applications. Most of the well-known manifestations, such as “wear out”, “fracture”,
“spalling”, and “section loss”, are related to the phenomenon of degradation [Batchelor et al., 2003].
Virtually, every material degrades when subjected to hostile environment and external stimuli.
Importance of this phenomena has triggered a surge in research to develop more resistible materials. Consequently, understanding the general behavior of degrading materials has attracted
the interest of researchers. A fundamental study of degradation is crucial to several branches of
engineering: aerospace, mechanical, civil, and biomedical. Moreover, some new materials, such as
ﬁber-reinforced polymers and multi-functional materials that exhibit high ductility have been widely
used recently; for example, as infrastructural materials or in medical devices (e.g., stents). In order
to model these materials, the traditional small-strain assumption will not be suﬃcient anymore.
In a nutshell, degradation means the loss in either serviceability or functionality. To be precise,
a material is said to be undergoing thermal degradation at a spatial point x ∈ Ω if the available
isothermal density is lower than the reference available isothermal power at that particular point.
That is,

M

dA
dt

6
ϑ>ϑref

dA
dt

for x ∈ Ω

(1.1)

ϑ=ϑref

Similarly, the chemical/moisture degradation can be deﬁned as follows:
dA
dt

6
c>cref

dA
dt

for x ∈ Ω

(1.2)

c=cref

where A denotes the speciﬁc Helmholtz potential of the material. Ω is the degrading body under
consideration, t is the time of interest, ϑref and cref are the speciﬁed reference temperature and
reference concentration. Note that degradation not only reduces the durability of materials but
also alters material properties. For instance, material damage can induce anisotropy in thermal
conductivity and diﬀusivity [Peng and Landel, 1975; Venerus et al., 2004; Zheng et al., 2011].
Herein, we develop a coupled continuum mathematical model for thermal and chemical-induced
degradation of solids, which are initially hyperelastic. We now outline three main reasons for such
a need.
2

• There is irrefutable experimental evidence that many modern infrastructural materials used
in repair and retroﬁtting applications exhibit large deformations. For example, the popular
high-early-strength Engineered Cementitious Composites (ECC) are capable of delivering
a compressive strength of 21 MPa within 4 hours after placement. Moreover, the long-term
tensile strain capacity of ECC members is more than 2% [Li, 2006; Wang and Li, 2006].
• In order to understand degradation mechanisms due to moisture, chemical, and temperature, coupling at various levels is needed (which is due to balance laws, material parameters, boundary conditions, and initial conditions). With existing and popular multi-physics
packages such ABAQUS [Aba, 2014], ANSYS [Ans, 2015], and COMSOL [Com, 2014], it is
possible to couple certain degradation mechanisms to some extent at material parameters,
boundary conditions, and initial conditions. However, such packages do not oﬀer ﬂexibility to couple important heat and mass transfer terms in balance laws. This is of utmost
importance in capturing the eﬀects of chemo-thermo-mechano degradation.
• Finally, when a new model or a thermodynamic framework is developed, stability of the
solutions for the corresponding initial boundary value problem needs to be shown. However, such an analysis is rarely performed when a new degradation model/framework is
developed in literature. Herein, for the proposed degradation framework we shall perform
stability analysis in the sense of Lyapunov. Subsequently, this methodology shall be used
to construct a robust computational framework in the part-II of the paper.
Hence, due to the above reasons small strain assumptions to model degradation and healing behavior
of these infrastructural systems are rarely valid. The proposed framework takes in to account the
underlying degradation mechanisms. Correspondingly, the respective parameters have a physical
meaning and can be calibrated through experiments.
It should be emphasized that elasticity is an idealization. There is no material whose response
is perfectly elastic. But there are situations in which the response of certain materials under
normal conditions can be idealized to be hyperelastic. For example, large blood arteries and rock.
Many of these materials function in hostile environments, and are constantly subjected to adverse
external stimuli. One often is interested in the unsteady response of the bodies made of hyperelastic
materials subjected to degradation/healing. The application areas in mind are the response of high
performance cementitious materials (which undergo large strains and large deformations) and several
important coupled deformation-thermal-transport processes in biomechanics and biomedicine. In
the next couple of subsections, we shall discuss various degradation mechanisms and the deﬁciencies
in the existing frameworks in modeling chemo-thermo-mechano degradation.
1.1. Degradation mechanisms. There are many mechanisms that can result in the degradation of materials and structures. In general, the degradation mechanisms can be divided into four
categories: mechanical processes, chemical reactions, biological degradation [Gu et al., 1998], and
radiation [Kaplan, 1989]. For mechanical processes, the performance of materials can be aﬀected
adversely by fatigue [Jung et al., 2000], pressure loading [Rajagopal et al., 2007], and swelling of
solid mixtures [Buonsanti et al., 2011]. Examples of chemical degradation include humid and alkaline eﬀects [Björk et al., 2003], exposure to chlorides and carbon-dioxide [Glasser et al., 2008], and
calcium leaching [Gawin et al., 2009]. Biological degradation refers to the dissolution of materials
by bacteria or other microorganisms. Degradation induced by radiation includes radiation damage
as well as other mechanical and chemical processes triggered by radiation.
3

Table 1. Various degradation mechanics and their primary manifestation. Many other
factors can be found in [Naus, 2007].

Degradation factor
Primary manifestation
cracking
reduced durability
vibration
cracking
Physical processes
freezing and thawing
cracking/scaling/disintegration
abrasion/erosion/cavitation
section loss
thermal exposure/thermal cycling cracking/spalling/strength loss
eﬄorescence/leaching
increased porosity
phosphate
surface deposits
Chemical processes
sulfate attack
volume change/cracking
acids/bases
disintegration/spalling/leaching
alkali-aggregate reactions
disintegration/cracking
The coupling eﬀects between these mechanisms can have a signiﬁcant impact on the rate of
deterioration of materials and structures. For instance, see Table 1 for some important factors that
aﬀect the degradation modeling in infrastructural materials such as concrete. Therefore, developing
an appropriate and general model for material degradation is useful to predict the life span of a
given structure. A comprehensive understanding of chemo-thermo-mechano degradation not only
plays a pivotal role in improving the quality and reliability of existing infrastructure, but also
has a tremendous impact on the economy [Herrmann, 2013]. In this paper, we shall assume that
predominant degradation mechanisms are moisture and temperature. We propose a general threeway strongly coupled degradation model based on a thermodynamic framework. This three-way
coupling is between mechanical, thermal, and transport processes.
1.2. Thermodynamics of chemo-thermo-mechano degradation. Herein, we shall provide a brief review and current status of thermal and chemical degradation. In the literature,
thermal degradation is modelled based on variants of thermoelasticity by incorporating damage
variables. Some popular research works in this direction are [Willam et al., 2005] for modeling
thermo-mechanical damage processes in heterogeneous cementitious materials and [Allam et al.,
2013] on the behavior of reinforced concrete slabs exposed to ﬁre. On the other hand, some popular
research works for the chemical degradation are [Björk et al., 2003] on the environmental eﬀects of
alkalinity and humidity on concrete slabs, [Cho and Kim, 2010] on moisture damage mechanisms occurring within asphaltic materials and pavements, [Bouadi and Sun, 1990] on thermal and moisture
eﬀects on structural stiﬀness and damping of laminated composites, and [Weitsman and Guo, 2002;
Weitsman, 2006] on ﬂuid-induced damage and absorption in polymeric composites. However, none
of the above mentioned papers on thermal or chemical degradation have a proper thermodynamic
basis.
There are two popular approaches to construct thermodynamically-consistent degradation models. The ﬁrst approach is based on the theory of the internal variable, wherein a scalar (or a tensor)
variable is introduced to model the degree of damage [Weitsman, 1987; Grasberger and Meschke,
2004; Springman and Bassani, 2009; Rajagopal et al., 2007]. For instance, the damage variable may
represent the measure of the fraction of broken cross-links or micro-cracks in a representative volume
element of the body [Kachanov, 1986; Lemaitre and Desmorat, 2005; Voyiadjis and Kattan, 2005].
4

The main disadvantage of this approach is that it is diﬃcult (or sometimes impossible) to measure
the internal variables through experiments or associate them to physical quantities/parameters.
The second approach is to build a thermodynamic framework by modeling all the relevant coupled processes. This achieved by taking into account the dependence of material properties on the
deformation of the solid, temperature, and concentration of chemical species. The degradation parameters under this approach have physical basis and can be calibrated using experiments (for example, see Section 5 of this paper). Herein, we shall employ the second approach to develop a thermodynamically consistent degradation model. It should be noted that certain research works exist in literature wherein the degradation models using the second approach. For example, see [Muliana et al.,
2009; Darbha and Rajagopal, 2009; Karra and Rajagopal, 2012; Klepach and Zohdi, 2014]. However, it appears that the above cited works suﬀer from the main drawback that they considered
thermodynamics of chemo-thermo-mechano degradation in the context of a closed system as opposed to an open system, which is the approach taken in this paper. Moreover, the models are not
as comprehensive as the one proposed in this paper.

1.3. Scope of the paper. In this paper we set out to achieve the following objectives:
(i) We derive a general chemo-thermo-mechano degradation model by appealing to the maximization of rate of dissipation. It will also be shown that many popular models are special
cases of the proposed mathematical model. For example, we will show that the small-strain
moisture degradation model proposed in [Mudunuru and Nakshatrala, 2012] is a special case
of the proposed model.
(ii) We will calibrate the proposed degradation model with existing experimental data sets. This
calibration study should provide conﬁdence in employing the proposed constitutive model to
model degradation of various brittle and quasi-brittle materials like ceramics, glass ﬁbers, and
concrete.
(iii) A systematic mathematical analysis is presented for the proposed model under large/ﬁnite
deformations. In particular, we shall show that the unsteady solutions under the proposed
degradation model are bounded and are stable in the sense of Lyapunov.
(iv) Last but not the least, semi-analytical solutions to several canonical problems are presented,
which provide insights into the behavior of degrading structural members. This will be valuable
for developing better design and safety codes.
The rest of the paper is organized as follows. Section 2 introduces the notation, mathematical
preliminaries, and the relevant balance laws. Section 3 presents a mathematical model for degradation of materials due to moisture and temperature, which is valid even under ﬁnite deformations
and large strains. The constitutive relations are obtained by appealing to the maximization of
rate of dissipation hypothesis, which ensures that the constitutive model satisﬁes the second law
of thermodynamics a prior. In Section 4, the proposed model is calibrated with an experimental
dataset. The coupled initial boundary value problem arising from the proposed degradation model
is presented in Section 5. We also show the solutions of the proposed mathematical model are
bounded and stable. In Section 6, solutions to several canonical problems are presented to illustrate
the predictive capabilities of the proposed model, and to highlight the eﬀects of degradation on the
structural behavior. Finally, conclusions are drawn in Section 7.
A list of the main symbols used in the paper are provided in the Nomenclature.
5

2. NOTATION, PRELIMINARIES, AND BALANCE LAWS
Let us consider a body B. The body occupies a reference conﬁguration Ω0 (B) ⊂ Rnd , where
“nd” denotes the number of spatial dimensions. A point in the reference conﬁguration is denoted
by p ∈ Ω0 (B). We shall denote the time by t ∈ [0, T ], where T is the length of the time interval
of interest. Due to motion, the body occupies diﬀerent spatial conﬁgurations with time. We shall
denote the conﬁguration occupied by the body at time t as Ωt (B) ⊂ Rnd . A corresponding spatial
point will be denoted as x ∈ Ωt (B). The gradient and divergence operators with respect to p are,
respectively, denoted by Grad[•] and Div[•]. Similarly, the gradient and divergence operators with
respect to x are, respectively, denoted by grad[•] and div[•].
The motion of the body is mathematically described by the following invertible mapping:
x = ϕ(p, t)

(2.1)

The displacement vector ﬁeld can then be written as:
u = x − p = ϕ(p, t) − p

(2.2)

The velocity vector ﬁeld is deﬁned as:
∂ϕ(p, t)
(2.3)
∂t
where a superposed dot indicates the material/total time derivative, which is the derivative with
respect to time holding the reference coordinates ﬁxed. The gradient of motion (which is also
referred to as the deformation gradient) is deﬁned as:
v = ẋ :=

∂ϕ(p, t)
= I + Grad[u]
(2.4)
∂p
where I denotes the second-order identity tensor. The corresponding right Cauchy-Green tensor is
denoted by:
F = Grad[x] ≡

C = FT F

(2.5)

where (•)T denotes the transpose of a second-order tensor. The velocity gradient with respect to x
and the symmetric part of the velocity gradient are, respectively, deﬁned as follows:
L := grad[v] ≡ ḞF−1

1
D :=
L + LT
2
The Green-St. Venant strain tensor is deﬁned as:

1
1
Grad[u] + Grad[u]T + Grad[u]T Grad[u]
E = (C − I) =
2
2
In situations the following assumption holds:
p
max
kϕ(p, t) − pk2 + kGrad[u]k2 ≪ 1
p∈Ω0 (B),t∈[0,T ]

(2.6)
(2.7)

(2.8)

(2.9)

one is justiﬁed to employ the following linearized strain tensor:
 1

1
El =
Grad[u] + Grad[u]T ≈
grad[u] + grad[u]T
(2.10)
2
2
where k • k denotes the Frobenius norm [Antman, 1995].
Since we will be dealing with processes in addition to the mechanical deformation, we need to
introduce quantities other than the ones that are associated with the kinematics. We will denote the
6

temperature by ϑ and the speciﬁc entropy by η. The mass fraction of the chemical species is denoted
by c and the corresponding chemical potential is denoted by κ. The temperature, mass fraction
of chemical species, entropy, and chemical potential are all scalar ﬁelds, while the displacement,
velocity, and acceleration are vector ﬁelds. In some situations, it may be needed to explicitly
indicate the functional dependence of these quantities. We employ a standard notation, which will
be illustrated through the temperature ﬁeld. The temperature in terms of reference coordinates and
spatial coordinates will be denoted as follows:
ϑ = ϑ̃(p, t) = ϑ̂(x, t)

(2.11)

2.1. Balance laws. For our study, we shall consider the thermodynamic system to be the
entire degrading body. Moreover, we shall assume this thermodynamic system to be an open
system. That is, heat and mass transfers can occur across the boundary of the system. We now
present the balance laws that govern the evolution of the chosen system.
The balance of mass of the solid in the degrading body takes the following form:
ρ̇ + ρ div[v] = 0

(2.12)

where ρ is the density of the solid in the deformed conﬁguration Ωt (B). The balance of a chemical
species, which is being transported in the degrading body, can be mathematically written as:
ρċ + div[h] = h

(2.13)

where h is the mass transfer ﬂux vector in the deformed conﬁguration, and h is the volumetric source
of the chemical species in the deformed conﬁguration. We assume that the chemical species cannot
take partial stresses, which is a reasonable assumption in the degradation of materials due to small
concentrations of moisture. One can handle large moisture contents by introducing partial stresses
and using the theory of interacting continua (which is often referred to mixture theory) [Bowen,
1976]. We do not address such issues, as our focus is degradation due to small concentrations of
moisture or chemicals. The balance of linear momentum of the solid can be written as:
ρv̇ = div[T] + ρb

(2.14)

where b is the speciﬁc body force, and T is the Cauchy stress in the solid. Assuming that there is
no internal couples, the balance of angular momentum of the solid reads:
T = TT

(2.15)

Assuming that the balance of linear momentum (i.e., equation (2.14)) holds, the balance of energy
of the system (i.e., the ﬁrst law of thermodynamics) can be written as:
d
(A + ϑη) = T • D − div[κh] + κh − div[q] + q
(2.16)
dt
where A is the speciﬁc Helmholtz potential, q is the heat ﬂux vector in the deformed conﬁguration,
and q is the volumetric heat source in the deformed conﬁguration. In our study, we assume that
the Helmholtz potential A to depend on F, c, and ϑ. We also have the following relations for the
chemical potential and speciﬁc entropy:
ρ

∂A
∂c
∂A
η := −
∂ϑ

κ := +

7

(2.17)
(2.18)

Assuming the balance of chemical species to hold, we then have the following convenient form for
the balance of energy:
ρ



∂A T
F • D + ϑη̇
∂F



= T • D − div[q] − grad[κ] • h + q

(2.19)

The localized version of the second law of thermodynamics in the deformed conﬁguration (by assuming that all the aforementioned balance laws to hold) takes the following form:
ρ




∂A T
1
F • D = T • D − grad[ϑ] • q − grad[κ] • h − ρζ
∂F
ϑ

(2.20)

where ζ is the speciﬁc rate of dissipation functional, which is non-negative. The above equation is
a stronger version than the second law of thermodynamics, which is a global law and not a local
one. The second law of thermodynamics does not assert that the rate of entropy production be nondecreasing at each and every point in the system/body. Strictly speaking, equation (2.20) should be
referred to as the reduced local dissipation equality. Another point to highlight is that the second
law of thermodynamics, in its original form, is in the form of an inequality. The introduction of the
non-negative dissipation functional, which acts as a slack variable, converts the inequality into an
equality, as provided in equation (2.20).

2.2. The maximization of rate of dissipation. Among the various methodologies to derive
constitutive relations (e.g., see [Maugin, 1998]), the axiom of maximization of rate of dissipation
put-forth by Ziegler [Ziegler, 1983] is an attractive procedure. Herein, we extend this procedure to
the open thermodynamic system under consideration. We obtain the constitutive relations using
the maximization of rate of dissipation hypothesis, which needs the prescription of two functionals –
the Helmholtz potential and the dissipation functional. We assume the functional dependence of the
Helmholtz potential and the dissipation functional to be Â(F, c, ϑ) and ζ̂(D, grad[ϑ], grad[κ]; F, ϑ, c).
The mathematical statement of maximization of rate of dissipation can be written as follows:
maximize

D,grad[ϑ],grad[κ]

subject to

ρζ = ρζ̂(D, grad[ϑ], grad[κ]; F, ϑ, c)


∂A T
1
ρ
F • D = T • D − grad[ϑ] • q − grad[κ] • h − ρζ
∂F
ϑ

(2.21a)
(2.21b)

Note that ρζ is maximized with respect to arguments to the left side of ‘;’. Using the method
of Lagrange multipliers, the above constrained optimization problem is equivalent to the following
unconstrained optimization problem:
extremize

ρζ̂(D, grad[ϑ], grad[κ]; F, ϑ, c)


 
1
∂A T
F • D − T • D + grad[ϑ] • q + grad[κ] • h + ρζ
+ Λt ρ
∂F
ϑ

D,grad[ϑ],grad[κ],Λt

8

(2.22)

where Λt is the Lagrange multiplier enforcing the constraint (2.21b). The ﬁrst-order optimal conditions give rise to the following relations:


∂A T
1 + Λt
∂ζ
T=ρ
F +
(2.23a)
ρ
∂F
Λt
∂D


1
1 + Λt
∂ζ
q=−
(2.23b)
ρ
ϑ
Λt
∂grad[ϑ]


∂ζ
1 + Λt
(2.23c)
ρ
h=−
Λt
∂grad[κ]


∂A T
1
ρ
F • D − T • D + grad[ϑ] • q + grad[κ] • h + ρζ = 0
(2.23d)
∂F
ϑ
The above equations can be obtained by taking (Gâteaux) variation of the objective function in
equation (2.22) with respect to D, grad[ϑ], grad[κ] and Λt , respectively. By straightforward manipulations on equations (2.23a)–(2.23d), the Lagrange multiplier Λt can be explicitly calculated as
follows:
"
#−1
ζ
Λt = ∂ζ
−1
(2.24)
∂ζ
∂ζ
∂D • D + ∂grad[ϑ] • grad[ϑ] + ∂grad[κ] • grad[κ]
If the rate of dissipation functional ζ is a homogeneous functional of order 2 with respect to D,
grad[ϑ] and grad[κ], we then have
∂ζ
∂ζ
∂ζ
•D+
• grad[ϑ] +
• grad[κ] = 2ζ
∂D
∂grad[ϑ]
∂grad[κ]

(2.25)

which further implies that Λt = −2. The constitutive relations under Λt = −2 will simplify to:
∂A T 1 ∂ζ
F + ρ
∂F
2 ∂D
ϑ
∂ζ
q=− ρ
2 ∂grad[ϑ]
∂ζ
1
h=− ρ
2 ∂grad[κ]

T=ρ

(2.26a)
(2.26b)
(2.26c)

Remark 2.1. It should be emphasized that the dissipation functional need not be a homogeneous
functional of order two in terms of F, c and ϑ. The maximization of the rate of dissipation certainly
does not require such an assumption. However, we shall make such an assumption, as it is convenient
and the resulting constitutive relations can still model the desired degradation mechanisms.
2.3. Governing equations in the reference configuration. Since we are also interested in
developing a computational framework and obtaining numerical solutions, it will be convenient to
write the balance laws in the reference conﬁguration. To this end, we introduce:
J ≡ det[F]

(2.27)

where det[•] denotes the determinant. The balance of mass in the reference conﬁguration can be
written as:
ρ0 = Jρ
9

(2.28)

where ρ0 is the density of the undeformed solid. The balance of chemical species in the reference
conﬁguration can be rewritten as:
ρ0 ċ + Div[h0 ] = h0

(2.29)

where h0 = JF−1 h is the diﬀusive ﬂux vector in the reference conﬁguration and h0 = Jh is the
volumetric source in the reference conﬁguration. The balance of linear momentum in the reference
conﬁguration takes the following form:
ρ0 v̇ = Div[P] + ρ0 b

(2.30)

where P = JTF−T is the ﬁrst Piola-Kirchhoﬀ stress. The balance of angular momentum in the
reference conﬁguration takes the following form:
PFT = FPT

(2.31)

In the reference conﬁguration, the balance of energy can be written as:
ρ0



∂A
• Ḟ + ϑη̇
∂F



= P • Ḟ − Div[q0 ] − Grad[κ] • h0 + q0

(2.32)

where q0 = JF−1 q is the heat ﬂux vector in the reference conﬁguration and q0 = Jq is the volumetric
heat source in the reference conﬁguration. In the reference conﬁguration, the second law can be
rewritten as:
ρ0




1
∂A
• Ḟ = P • Ḟ − Grad[ϑ] • q0 − Grad[κ] • h0 − ρ0 ζ0
∂F
ϑ

(2.33)

where ζ0 = ζ is the non-negative rate of dissipation functional in the reference conﬁguration.
2.3.1. Maximization of rate of dissipation in the reference conﬁguration. The mathematical
statement of maximization of rate of dissipation can be written as follows:
maximize

ρ0 ζ0 = ρ0 ζ̃(Ḟ, Grad[ϑ], Grad[κ]; F, ϑ, c)

Ḟ,Grad[ϑ],Grad[κ]

subject to

ρ0



∂A
• Ḟ
∂F



= P • Ḟ −

1
Grad[ϑ] • q0 − Grad[κ] • h0 − ρ0 ζ0
ϑ

(2.34a)
(2.34b)

Using the method of Lagrange multipliers, one can obtain the following equivalent unconstrained
optimization problem:
extremize

ρ0 ζ̃(Ḟ, Grad[ϑ], Grad[κ]; F, ϑ, c)

Ḟ,Grad[ϑ],Grad[κ],Λ0



+ Λ0 ρ0





1
∂A
• Ḟ − P • Ḟ + Grad[ϑ] • q0 + Grad[κ] • h0 + ρ0 ζ0
∂F
ϑ
10

(2.35)

where Λ0 is the Lagrange multiplier enforcing the constraint given by equation (2.34b). The ﬁrstorder optimality conditions give rise to the following constitutive relations:


∂A
1 + Λ0
∂ζ0
P = ρ0
(2.36a)
+
ρ0
∂F
Λ0
∂ Ḟ


1 + Λ0
∂ζ0
1
q0 = −
(2.36b)
ρ0
ϑ
Λ0
∂Grad[ϑ]


1 + Λ0
∂ζ0
h0 = −
(2.36c)
ρ0
Λ0
∂Grad[κ]


∂A
1
ρ0
• Ḟ − P • Ḟ + Grad[ϑ] • q0 + Grad[κ] • h0 + ρ0 ζ0 = 0
(2.36d)
∂F
ϑ
Similar to the derivation presented earlier in the context of current conﬁguration, the Lagrange
multiplier Λ0 can be explicitly calculated as follows:
"
#−1
ζ0
Λ0 = ∂ζ0
−1
(2.37)
∂ζ0
∂ζ0
•
Grad[ϑ]
+
•
Grad[κ]
•
Ḟ
+
∂Grad[ϑ]
∂Grad[κ]
∂ Ḟ
If the rate of dissipation functional in the reference conﬁguration ζ0 is a homogeneous functional of
order 2, we have
∂ζ0
∂ζ0
∂ζ0
• Ḟ +
• Grad[ϑ] +
• Grad[κ] = 2ζ0
∂Grad[ϑ]
∂Grad[κ]
∂ Ḟ

(2.38)

which further implies that Λ0 = −2. The constitutive relations under Λ0 = −2 take the following
form:
∂A 1 ∂ζ0
+ ρ0
∂F 2 ∂ Ḟ
ϑ
∂ζ0
q0 = − ρ0
2 ∂Grad[ϑ]
∂ζ0
1
h0 = − ρ 0
2 ∂Grad[κ]

P = ρ0

(2.39a)
(2.39b)
(2.39c)

The overarching idea behind the proposed chemo-thermo-mechano degradation model is shown
in Figure 1. In the next section, we will develop the proposed constitutive model by appealing to
the maximization of rate of dissipation. Solving the coupled balance laws (even including the evolution equation for internal variable), we can get the displacement, temperature, and concentration
(internal variable, if needed).
3. A GENERAL CONSTITUTIVE MODEL FOR CHEMO-THERMO-MECHANO
DEGRADATION
Under the maximization of rate of dissipation hypothesis, the constitutive relations can be
obtained by prescribing two functionals – the Helmholtz potential and the dissipation functional.
Philosophically, the Helmholtz potential quantiﬁes the way in which the material stores energy,
whereas the dissipation functional quantiﬁes the way in which the material dissipates energy. For our
proposed chemo-thermo-mechano degradation model, we prescribe the following functional forms
11

Rate of dissipation
functional

Helmholtz potential

Maximization of rate
of dissipation

Balance of mass

Stress

Heat flux

Diffusion flux

Balance of linear
momentum

Balance of
energy

Balance of
chemical species

Evolution
equation for α

Solve coupled equations
Displacement, temperature,
concentration, and internal variable

Figure 1. Overarching idea of the proposed degradation framework: This flowchart shows
the overarching idea behind the proposed framework. We shall appeal to the axiom of
maximization of rate of dissipation to obtain constitutive relations for stress, heat flux,
diffusion flux, and evolution equation for internal variable (if required). Solving the coupled
equations, we get the solution for displacement, temperature, concentration, and internal
variable (if required).

for the speciﬁc Helmholtz potential and the rate of dissipation functional:
1 cp
1
1
{ϑ − ϑref }2 −
{ϑ − ϑref } MϑE • E + dϑc {ϑ − ϑref } {c − cref }
A = Â(F, c, ϑ) = ψ −
ρ0
2 ϑref
ρ0
Rs ϑref
1
{c − cref }2
(3.1)
{c − cref } McE • E +
−
ρ0
2
cp
1
grad[ϑ] • Dϑϑ grad[ϑ] + grad[ϑ] • Dϑκ grad[κ]
ϑ
ϑ
1
1
+ grad[κ] • Dκϑ grad[ϑ] +
grad[κ] • Dκκ grad[κ] (3.2)
ϑ
Rs ϑref

ζ = ζ̂(D, grad[ϑ], grad[κ]; F, ϑ, c) =

where Rs = R/M . Rs and R denote the speciﬁc vapor constant and the universal vapor constant
respectively, M is the molecular mass of chemical species. ϑref and cref are the speciﬁed reference
temperature and reference mass concentration, which depend on the underlying boundary value
problem. We denote cp as the coeﬃcient of heat capacity, dϑc as the thermo-chemo coupled parameter [Sih et al., 1986], MϑE as the anisotropic coeﬃcient of thermal expansion (which is assumed to
be independent of temperature, concentration, and strain), and McE as the anisotropic coeﬃcient
of chemical expansion due to concentration (which is also assumed to be independent of temperature, concentration, and strain). Both MϑE and McE are assumed to be symmetric. Dϑϑ is the
anisotropic thermal conductivity tensor and Dκκ is the anisotropic diﬀusivity tensor. Dϑκ corresponds to the anisotropic Soret eﬀect tensor, which characterizes the transport of chemical species
caused by temperature gradient. Similarly, Dκϑ is the Dufour eﬀect tensor, which represents the
heat ﬂow caused by a concentration gradient.
Remark 3.1. In chemo-thermo-elasticity and in modeling degradation of materials due to transport and reaction of chemical species, coeﬃcient of chemical expansion McE and thermo-chemo
coupling parameter dϑc play a vital role (see [Sih et al., 1986, Chapter-5] and references therein).
12

Induced-strains due to chemical expansivity will be signiﬁcant in harsh environmental conditions and
cannot be neglected [Sih et al., 1986]. Considerable inquest has been made in literature to experimentally measure McE in ceramics [Adler, 2001; Morozovska et al., 2011; Blond and Richet, 2008], laminated and polymer composites [Sih et al., 1986; Bouadi and Sun, 1989; Cai and Weitsman, 1994],
elastomers and biological materials [Harper, 2002; Myers et al., 1984; Lai et al., 1991], and concrete structures [Ulm et al., 2000; Černy and Rovnaníková, 2002; Swamy, 2002]. However, adequate
progress has not been made yet to develop constitutive models and computational frameworks for such
chemo-thermo-elastic materials or materials undergoing chemical-induced degradation. Herein, we
shall take a step forward to address this issue.
Remark 3.2. It should be noted that in the absence of electrical and magnetic ﬁelds, all of the
above tensors are symmetric [Bowen, 1976; Coussy, 2004; Jarkova et al., 2001]. Moreover, from
the Onsager reciprocal relations (which was put-forth by Onsager in 1930s [Onsager, 1931a,b]) we
have the following relationship between the Soret eﬀect tensor and the Dufour eﬀect tensor.
Dϑκ = Dκϑ

(3.3)

Additionally, physics demands that the tensors Dϑϑ and Dκκ are positive deﬁnite.
Remark 3.3. Note that the speciﬁc Helmholtz potential and correspondingly the dissipation
functional for diﬀusion can also be modelled using the following expressions:
Ac = Rs ϑref c{ln[c] − 1}
(3.4)
c
grad[κ] • Dκκ grad[κ]
(3.5)
ζc =
Rs ϑref
Both equations (3.1)–(3.2) and (3.4)–(3.5) result in similar partial diﬀerential equation structure for
modeling Fickian diﬀusion.
Under the proposed model, the speciﬁc entropy and chemical potential take the following form:
∂A
1 ∂ψ
cp
1
η=−
=−
+
{ϑ − ϑref } + MϑE • E − dϑc {c − cref }
(3.6)
∂ϑ
ρ0 ∂ϑ ϑref
ρ0
1 ∂ψ
1
∂A
=
+ Rs ϑref {c − cref } − McE • E + dϑc {ϑ − ϑref }
(3.7)
κ=
∂c
ρ0 ∂c
ρ0
From equations (2.26a)–(2.26c), we have the constitutive relations in deformed conﬁguration as:
1 ∂ψ T 1
1
∂A T
F =
F − {ϑ − ϑref } FMϑE FT − {c − cref } FMcE FT
(3.8a)
T=ρ
∂F
J ∂F
J
J
∂ ζ̂
ρ
ρ
ϑ
= −ρcp Dϑϑ grad[ϑ] − Dϑκ grad[κ] − Dκϑ grad[κ]
(3.8b)
q=− ρ
2 ∂grad[ϑ]
2
2
1
∂ ζ̂
ρ
ρ
ρ
h=− ρ
=−
Dκκ grad[κ] −
Dϑκ grad[ϑ] −
Dκϑ grad[ϑ]
(3.8c)
2 ∂grad[κ]
Rs ϑref
2ϑ
2ϑ
The rate of dissipation functional for the degradation model in the reference conﬁguration is taken
as follows:
ζ = ζ̃(Ḟ, Grad[ϑ], Grad[κ]; F, ϑ, c)
cp
1
= Grad[ϑ] • Dϑϑ Grad[ϑ] + Grad[ϑ] • Dϑκ Grad[κ]
ϑ
ϑ
1
1
Grad[κ] • Dκκ Grad[κ]
+ Grad[κ] • Dκϑ Grad[ϑ] +
ϑ
Rs ϑref
13

(3.9)

where Dαβ = F−1 Dαβ F−T , α and β represent ϑ or κ. Correspondingly, the constitutive relations
in the reference conﬁguration take the following form:
∂ψ
∂A
=
− {ϑ − ϑref } FMϑE − {c − cref } FMcE
∂F
∂F
ϑ
∂ ζ̃
ρ0
ρ0
q0 = − ρ0
= −ρ0 cp Dϑϑ Grad[ϑ] − Dϑκ Grad[κ] − Dκϑ Grad[κ]
2 ∂Grad[ϑ]
2
2
P = ρ0

∂ ζ̃
ρ0
1
ρ0
ρ0
=−
h0 = − ρ0
Dκκ Grad[κ] −
Dϑκ Grad[ϑ] −
Dκϑ Grad[ϑ]
2 ∂Grad[κ]
Rs ϑref
2ϑ
2ϑ

(3.10a)
(3.10b)
(3.10c)

3.1. Coupling terms for the degradation model. The following hyperelastic material models will be employed in this paper:
λ
(tr[E])2 + µE • E
St. Venant-Kirchhoﬀ model
(3.11a)
2
κ
ψ = (ln[J])2 + µE • E
Modiﬁed St. Venant-Kirchhoﬀ model
(3.11b)
2
λ
ψ = µtr[E] + µln[J] + (ln[J])2 Neo-Hookean model
(3.11c)
2
where ψ is the stored strain energy density functional, λ and µ are the Lamé parameters, and
κ = λ + 2µ
3 is the bulk modulus. Recall that J = det[F]. The Lamé parameters in the degrading
model are given by the following expressions:
ψ=

ϑ
cref
ϑref
ϑ
c
− µ2 (x)
µ(x, c) = µ0 (x) − µ1 (x)
cref
ϑref
λ(x, c) = λ0 (x) − λ1 (x)

c

− λ2 (x)

(3.12a)
(3.12b)

where λ0 and µ0 are the Lamé parameters of the virgin material. λ1 and µ1 are the parameters that
account for the eﬀect of concentration of chemical species on degradation of solid. λ2 and µ2 are
the parameters that account for the temperature eﬀect on the degrading solid. It should be noted
that λ1 , µ1 , λ2 , and µ2 are all positive. Furthermore, these parameters are constrained such that
the bulk modulus and shear modulus are strictly positive.
3.1.1. Deformation dependent diﬀusivity. The eﬀect of deformation on diﬀusivity is modeled as
follows: When tensile and shear strains are predominant, we have the following constitutive model
(exp[ηS IIE ] − 1)
(exp[ηT IE ] − 1)
+ (DS − D0 )
(exp[ηT ErefT ] − 1)
(exp[ηS ErefS ] − 1)
(exp[ηM S IIIE ] − 1)
− D0 )
(exp[ηM S ErefM S ] − 1)

Dκκ = D0 + (DT − D0 )
+ (DM S

(3.13)

where IE , IIE , and IIIE are the ﬁrst, second, and third invariants of the Green-St. Venant strain
tensor. These are deﬁned as follows:
IE := tr[E]
p

2 dev[E] • dev[E] =


1
dev[E]
IIIE := det
IIE
IIE :=

14

r

(3.14a)
2
(3tr[E2 ] − (tr[E])2 )
3

(3.14b)
(3.14c)

where dev[E] := E − 31 tr[E]I is the deviatoric part of E. These invariants are used to model the
eﬀect of dilation, magnitude of distortion, and mode of distortion on the diﬀusivity of the solid.
ηT , ηS , and ηM S are non-negative parameters. ErefT , ErefS , and ErefM S are reference measures of
the tensile strain, shear strain, and mode of shear strain respectively. D0 , DT , DS , and DM S are,
respectively, the reference diﬀusivity tensors under no strain, tensile strain, and shear strain.
When compression and shear strains are predominant, deformation dependent diﬀusivity is
modeled as follows:
(exp[ηS IIE ] − 1)
(exp[ηT IE ] − 1)
+ (DS − D0 )
(exp[ηT ErefT ] − 1)
(exp[ηS ErefS ] − 1)
(exp[ηM S IIIE ] − 1)
− D0 )
(exp[ηM S ErefM S ] − 1)

Dκκ = D0 + (D0 − DC )
+ (DM S

(3.15)

where ηC is a non-negative parameter, ErefC is a reference measure of the compression strain, and
DC is the reference diﬀusivity tensor under compressive strain.
Remark 3.4. Note that deformation dependent diﬀusivity given by equations (3.13) and (3.15)
can be constructed using a diﬀerent set of invariants of a given strain tensor. This invariants can
be either principal or Hencky type [Lurie, 1990; sek and Kruisová, 2006; Criscione et al., 2000]
based on the nature of material and associated experimental data. The proposed framework can
accommodate such models with minor modiﬁcations.
In case of transversely isotropic materials with ﬁbers running along the direction Mtf , the following invariants are needed to model deformation dependent diﬀusivity in addition to the invariant
set given by equations (3.14a)–(3.14c)
IVE := Mtf • EMtf
VE := Mtf • E2 Mtf

(3.16a)
(3.16b)

For more details on selection of invariants for transversely isotropic or orthotropic materials see
[Lurie, 1990; Holzapfel, 2000; Ogden, 1997].
3.1.2. Deformation dependent thermal conductivity. The eﬀect of deformation of the solid on
thermal conductivity is modeled as follows [Bhowmick and Shenoy, 2006]:
Dϑϑ = K0ϑ (1 + IE )−δ

(3.17)

where δ is a non-negative parameter. K0ϑ is the reference conductivity tensor under no strain.
Based on molecular dynamics simulations, Bhowmick and Shenoy [Bhowmick and Shenoy, 2006]
suggested δ to be 9.59 and K0ϑ = 4.61ϑ−1.45 (for certain brittle-type materials). For various other
ductile or brittle-type materials, these parameters can be determined by experiments or can be
constructed using Lennard-Jones potential in molecular dynamics.
Remark 3.5. Due to the lack of experimental data, we assume that Dufour and Soret tensors do
not depend on the deformation of solid. However, it should be noted that the proposed thermodynamic
and computational framework can accommodate deformation dependent Dufour and Soret tensors
with minor modiﬁcations (whenever such an experimental evidence is available).
3.2. Special cases of the general degradation model and their thermodynamic status. The following popular non-degradation constitutive models can be shown as special cases of
the proposed degradation model, as shown in Figure 2, when the material parameters are assumed
15

to be independent of concentration, temperature, and deformation of the solid. That is, the Lamé
parameters and Dαβ (α and β represent either ϑ or κ) are independent of c, ϑ, and E.
(1) Fourier and Fickian models: The standard heat conduction constitutive model is obtained by assuming the solid to be rigid and mass concentration of diﬀusing chemical
species to be equal to zero. Similarly, to recover the standard Fickian model we assume
the solid to be rigid and temperature of the homogenized body to be constant.
(2) Dufour-Soret model: This model is obtained by assuming the solid to be rigid. Furthermore, the thermo-chemo coupling parameter dϑc is neglected.
(3) Linearized elasticity and hyperelasticity: To obtain hyperelastic constitutive models, we assume isothermal conditions and mass concentration of diﬀusing chemical species
to be equal to zero. The linearized elasticity model can be recovered from any given hyperelastic model by assuming that the small strains assumption given by equation (2.9)
holds.
(4) Thermoelasticity: The standard thermoelasticity model can be recovered by assuming
mass concentration of diﬀusing chemical species to be equal to zero. The material parameters are assumed to be independent of temperature and deformation.
(5) Chemoelasticity: Similarly, the standard chemoelasticity model can be recovered by assuming isothermal conditions. The material parameters are assumed to be independent of
concentration and deformation.
(6) Chemo-Thermoelasticity: Herein, we assume that the material parameters are independent of concentration, temperature, and deformation. In addition, thermo-chemo coupling
parameter dϑc , Dufour tensor, and Soret tensor are neglected.
One can also derive specialized (thermo-mechano and chemo-mechano) degradation models:
(1) Thermo-mechano degradation model: This model is obtained from the thermoelasticity
model by relaxing the assumption that material parameters are independent of temperature
and deformation.
(2) Chemo-mechano degradation model: Similar to thermo-mechano degradation model, this
degradation model is obtained from the chemoelasticity model by relaxing the assumption
that material parameters are independent of concentration and deformation.
3.2.1. Status of the degradation model in [Mudunuru and Nakshatrala, 2012]. The small-strain
chemo-mechano degradation model proposed in [Mudunuru and Nakshatrala, 2012] is a special case
of the proposed chemo-thermo-mechano degradation, and can be obtained under a plethora of
assumptions. These assumptions include steady-state response, small strains, and isothermal conditions with negative volumetric heat source in the entire degrading body. One also needs to neglect
chemo-thermo, chemo-mechano, and thermo-mechano couplings. Moreover, the functional forms of
the speciﬁc Helmholtz potential and rate of dissipation functional need to be:
Rs ϑref
1
ψ+
{c − cref }2
ρ
2
1
grad[κ] • Dκκ grad[κ]
ζ=
Rs ϑref

A=

(3.18)
(3.19)

where the stored strain energy density functional is given by:
ψ = ψ̂(El , c) =

λ(c)
(tr[El ])2 + µ(c)El • El
2
16

(3.20)

Proposed chemo-thermo-mechano
degradation model

Assumption
Model

No degradation

No diffusion

Thermal degradation

Chemo-thermoelasticity

Isothermal

Chemical degradation

No degradation

No diffusion,
isothermal

No deformation

No degradation

Thermo-elasticity

Hyper-elasticity

Dufour-Soret
model

Chemo-elasticity

No deformation

Small strain

No deformation

Fourier model

Linear elasticity

Fickian model

Figure 2. Special cases of the proposed chemo-thermo-mechano degradation model: Many
existing degrading and non-degrading constitutive models are special cases of the proposed
hierarchical model, with appropriate assumptions.

Under the small strain assumption given by equation (2.9), the Cauchy stress, chemical potential,
and mass transfer ﬂux vector can be written as:
∂A
= λ(c)tr[El ]I + 2µ(c)El
(3.21)
T=ρ
∂El
∂A
κ=
= Rs ϑref {c − cref }
(3.22)
∂c
1
∂ ζ̂
ρ
h=− ρ
=−
Dκκ grad[κ]
(3.23)
2 ∂grad[κ]
Rs ϑref
The balance of chemical species and the balance of linear momentum for the solid are given by
equations (2.13) and (2.14). Under the isothermal condition, the balance of energy simpliﬁes to the
following expression:
ρ
grad[κ] • Dκκ grad[κ]
(3.24)
q=−
Rs ϑref
which means that q needs to be non-positive in order to maintain the isothermal condition. The
deformation dependent diﬀusivity Dκκ is based on the small strain assumption, which is obtained
by linearizing the equations (3.13) and (3.15). Note that this model is developed based on the
experimental evidence that the relative diﬀusion rate varies exponentially with respect to the trace
of strain [McAfee, 1958a,b]. In this paper, we have taken a step further to calibrate these materials parameters according to the experimental data for ﬁnite strains based on the model given by
equations (3.13) and (3.15).
4. CALIBRATION WITH EXPERIMENTAL DATA
In this section, we will calibrate the proposed model for the diﬀusivity using the experimental
data set reported in [McAfee, 1958a,b]. These experiments were conducted on spherical shells made
of glass, which is a brittle material. These papers report the variation of diﬀusivity under various
17

0.01

pi
po = 0

Figure 3. Calibration with experimental data: A pictorial description of the degrading shell
used for calibrating the proposed model with the experimental data.

deformation modes: tension, compression, and shear. The calibration study presented below, which
also includes a statistical analysis of the ﬁt, will be valuable in two ways. First, it demonstrates
the predictive capabilities of the proposed constitutive model, and provides conﬁdence in the model
to be able to apply to other brittle materials like ceramics and even to quasi-brittle materials
like concrete with appropriate modiﬁcations. Second, it provides order-of-magnitude estimates for
various parameters in the diﬀusivity model for realistic materials. This will guide in the selection
of values for these parameters in the subsequent numerical studies. Figure 3 provides the geometry
and the loading on a spherical shell. The inner and outer radii are, respectively, a = 0.99 and
b = 1.0. The boundary conditions for the deformation subproblem is that the pressure within the
sphere is varied from pi = 0 to pi = 0.68947 MPa (100 psi) and the external surface is traction free.
For the diﬀusion subproblem, c(a) = 0 and c(b) = 1. In this scenario, it can be assumed that the
tensile strain is predominant. The diﬀusion can be assumed to be isotropic. Hence, equation (3.13)
is simpliﬁed as follows:
D = D0 + (DT − D0 )

(exp[ηT IE ] − 1)
(exp[ηT Eref T ] − 1)

(4.1)

The sample size to estimate the parameters in the proposed deformation-diﬀusivity model has been
taken to be 3. It has been reported that D0 = 7.26 × 10−13 m2 /sec−1 for glass ﬁbers by [McAfee,
1958a]. Based on the chosen sample size and value of D0 , the estimated diﬀusivity parameters are
given as follows:
ηT = 1.43 × 104 , DT = 23.39 × 10−13 , Eref T = 1.833 × 10−3

(4.2)

Using the experimental data reported in [McAfee, 1958b] under compressive and shear strains, and
following a similar procedure as before, the following diﬀusivity parameters are obtained:
ηC = 401.19, DC = 8.66 × 10−13 , Eref C = 1.0 × 10−3

(4.3a)

ηS = −239.61, DS = 8.65 × 10−13 , Eref S = 3.0 × 10−3

(4.3b)

We then compared the proposed model (which is obtained based on sample size of 3 points) with
the experimental data set of 10 points. Figure 4 shows the relation between the relative diﬀusion
coeﬃcient D/D0 and various strain invariants. From this ﬁgure, it is evident that the proposed
model is in a good agreement with the experimental data. Table 2 provides a statistical analysis on
the ﬁt of the experimental data with the proposed model. The coeﬃcient of determination is close
to 1. This means that the proposed model based on parameter set given by equations (4.2)–(4.3b)
is a good ﬁt to the set of experimental data of various sample sizes. To calibrate DM S , ηM S , and
18

D/D0

3
Experimental data for tension
Proposed model for tension
Experimental data for compression
Proposed model for compression
Experimental data for shear
Proposed model for shear

2

1
0

0.001
0.002
Trace of strain

0.003

Figure 4. Calibration with experimental data: This figure compares the experimental data
reported in [McAfee, 1958a,b] with the proposed constitutive model. The sample size is taken
to be 3. The strain invariants are given by equations (3.14a)–(3.14b). A good agreement
has been observed between the experimental data and the proposed constitutive model for
the diffusivity under tensile, compressive, and shear strains.
Table 2. A statistical analysis of the fit: This table provides the goodness of fit of the
proposed model with the experimental data set reported in [McAfee, 1958a,b]. Analysis
is performed for various extracted sample sizes, and under tension, compression and shear
strains. It is observed that the coefficient of determination is close to 1.
Sample size

10
25
50
75

Mean data

Standard deviation data

Coefficient of determination

Tension

Compression

Shear

Tension

Compression

Shear

Tension

Compression

Shear

1.507
1.505
1.521
1.391

1.093
1.094
1.095
1.097

1.114
1.108
1.107
1.115

0.526
0.511
0.515
0.468

0.073
0.062
0.062
0.062

0.068
0.065
0.062
0.059

0.988
0.986
0.987
0.984

0.999
0.999
0.999
0.999

0.997
0.997
0.997
0.997

ErefM S , we need additional experimental data related to mode of shear. However, such a data set to
calibrate the eﬀect of distortion due to shear on the diﬀusivity of glass is not available in literature
yet. Hence, we could not calibrate DM S , ηM S , and ErefM S . Whenever such an experimental data
is available, one can easily calibrate these parameters in a similar fashion.
5. INITIAL BOUNDARY VALUE PROBLEM AND MATHEMATICAL ANALYSIS
From the above statements, the governing equations for the proposed chemo-thermo-mechano
degrading model are stated as follows. Let the boundary of Ω0 (B) be denoted as ∂Ω0 and the
b 0 (p). Similarly, ∂Ωt denotes
corresponding unit outward normal to this boundary be denoted by n
the boundary of Ωt (B) and the corresponding unit outward normal to this boundary is denoted by
b (x, t). For the deformation subproblem, the boundary is divided into two complementary parts:
n
19

N
D
N
D
N
D
ΓD
u and Γu such that Γu ∪ Γu = ∂Ω0 and Γu ∩ Γu = ∅. Γu is the part of the boundary on which
displacement is prescribed and ΓN
u is the part of the boundary on which traction is prescribed.
Similarly, for the transport and thermal subproblem, the boundary is divided into complemenN
D
N
D
N
D
N
D
N
tary parts: ΓD
c and Γc and Γϑ and Γϑ such that Γc ∪ Γc = ∂Ω0 , Γϑ ∪ Γϑ = ∂Ω, Γc ∩ Γc = ∅, and
D
N
N
D
Γϑ ∩ Γϑ = ∅. Γc is the part of the boundary on which concentration is prescribed. Γc is the part of
the boundary on which total/diﬀusive ﬂux is prescribed. ΓD
ϑ is the part of the boundary on which
is
the
part
of
the
boundary
on which thermal ﬂux is prescribed.
temperature is prescribed. ΓN
ϑ
In case of steady-state analysis, it should be noted that the meas ΓD
> 0, meas ΓD
c
ϑ > 0, and
D
meas Γu > 0. However, such a condition is not required for studying transient problems.

5.1. Governing equations of the proposed model. The governing equations for the deformation sub-problem take the following form:
ρ0 v̇(p, t) = Div[P] + ρ0 b(p, t) in Ω0 ×]0, I[
p

(5.1a)

u(p, t) = u (p, t)

on

ΓD
u ×]0, I[

(5.1b)

Pb
n0 (p) = tp (p, t)

on ΓN
u ×]0, I[

(5.1c)

u(p, t = 0) = ui (p)

in Ω0

(5.1d)

v(p, t = 0) = vi (p)

in Ω0

(5.1e)

where up (p, t) denotes the prescribed displacement on the boundary and tp (p, t) is the prescribed
traction on the boundary. ui (p) and vi (p) are the initial conditions for the displacement and
velocity, respectively.
The governing equations for the transport sub-problem take the following form:
ρ0 ċ(p, t) + Div[h0 ] = h0 (p, t) in Ω0 ×]0, I[

(5.2a)

c(p, t) = cp (p, t)

on ΓD
c ×]0, I[

(5.2b)

b 0 (p) = hp (p, t)
h0 • n

on ΓN
c ×]0, I[

(5.2c)

in Ω0

(5.2d)

i

c(p, t = 0) = c (p)

where cp (p, t) denotes the prescribed concentration on the boundary, hp (p, t) is the prescribed
total/diﬀusive ﬂux on the boundary, and ci (p) is the initial condition for the concentration ﬁeld.
The governing equations for the thermal sub-problem take the following form:
ρ0 ϑ(p, t)η̇ = −Div[q0 ] − Grad[κ] • h0 + q0 (p, t) in Ω0 ×]0, I[

(5.3a)

ϑ(p, t) = ϑp (p, t)

on ΓD
ϑ ×]0, I[

(5.3b)

b 0 (p) = q p (p, t)
q0 • n

on ΓN
c ×]0, I[

(5.3c)

in Ω0

(5.3d)

ϑ(p, t = 0) = ϑi (p)

where ϑp (p, t) denotes the prescribed temperature on the boundary, q p (p, t) is the prescribed heat
ﬂux on the boundary, and ϑi (p) is the initial condition for the temperature ﬁeld.
5.2. On the stability of unsteady solutions. We now show that the unsteady solutions
under the proposed mathematical model for degradation are stable in the sense of a dynamical
system. There are diﬀerent notions of stability, and herein we shall establish the stability in the sense
of Lyapunov [Dym, 2002]. For the entire analysis presented in this section, we assume homogeneous
20

Dirichlet boundary conditions on the entire boundary for the diﬀusion and thermal sub-problems.
Let


ϕ 





v
χ :=

ϑ 




c

(5.4)

Consider the following functional, which is deﬁned on the reference conﬁguration:
V(χ) :=

Z

Ω0 (B)

ρ0



1
A + ϑη + v • v
2



dΩ0 + Πmech,ext (ϕ)

(5.5)

where Πmech,ext (ϕ) is the potential energy due to external mechanical loading, which is assumed to
be conservative. This implies the following
d
Πmech,ext (ϕ) = −
dt

Z

ρ0 b • v dΩ0 −
Ω0 (B)

Z

ΓN
u

tp • v dΓ0

(5.6)

In the literature, the above functional has been shown to be a Lyapunov functional for linearized
thermoelasticity and for themo-hyperelasticity. For example, see [Ericksen, 1966; Coleman and Dill,
1973; Gurtin, 1975] and references therein. Herein, we shall show that the above functional is
a legitimate Lyapunov functional for the proposed degradation model, and speciﬁcally use the
Lyapunov’s second method for stability (which is a classical result in the theory of dynamical
systems; for example, see [Hale and Kocak, 1991; Strogatz, 2001; Wiggins, 2003]) to establish the
stability of the solutions under the proposed degradation model.
To this end, we shall take the reference or equilibrium state as:

χeq


ϕeq



0
:=

0


0









(5.7)

where ϕeq is the static equilibrium deformation. The above functional is a candidate for Lyapunov
functional, as it satisﬁes:
V(χ = χeq ) = 0 and V(χ 6= χeq ) > 0

(5.8)

dV
≤0
dt

(5.9)

We now show that

21

dV
=
dt

Z

ρ0

Ω0 (B)

=

Z



∂A
∂A
∂A
• Ḟ +
ċ + ϑ̇η + ϑη̇ + v • v̇
ϑ̇ +
∂F
∂ϑ
∂c

ρ0 (κ ċ + ϑη̇) dΩ0 +

=

Z

κDiv[h0 ] dΩ0 −

Z

Grad[κ] • h0 dΩ0 −

Ω0 (B)

=

Z

Z


Z
∂A
+ η ϑ̇ dΩ0 +
∂ϑ

Ω0 (B)

ϑref
ϑ

Ω0 (B)



Z

ϑ − ϑref
Div[q0 ] dΩ0 −
ϑ

Ω0 (B)

Ω0 (B)

=



dΩ0 +

d
Πmech,ext (ϕ)
dt



d
ρ0 v • v̇ + P • Ḟ dΩ0 + Πmech,ext (ϕ)
dt

ρ0 (κ ċ + ϑη̇) dΩ0

Ω0 (B)

=−

ρ0

Ω0 (B)

Ω0 (B)

Z

Z



Z

Ω0 (B)



ϑref
1−
ϑ



1
Grad[ϑ] • q0 + Grad[κ] • h0
ϑ

ϑ − ϑref
Grad[κ] • h0 dΩ0
ϑ

Ω0 (B)

Div[q0 ] dΩ0 −



dΩ0 = −

Z

Z

Ω0 (B)



ϑref
1−
ϑ



Grad[κ] • h0 dΩ0

ϑref
ζ0 dΩ0
ϑ

(5.10)

Ω0 (B)

Since ζ0 > 0 if χ 6= χeq , ϑ, ϑref > 0, one can conclude that
dV
<0
dt

(5.11)

From the Lyapunov stability of continuous systems [Dym, 2002; Hale and Kocak, 1991], one can
conclude that χ = χeq is asymptotically stable.
6. SEMI-ANALYTICAL SOLUTIONS TO CANONICAL PROBLEMS
In this section, we shall appeal to semi-inverse methods to obtain solutions to some popular
canonical boundary value problems [Ogden, 1997]. Incompressible neo–Hookean chemo-thermomechano degradation model is considered here. Similar analysis can be performed for other compressible and incompressible chemo-mechano, thermo-mechano, and chemo-thermo-mechano degradation models. Coordinate system under consideration is either spherical or cylindrical. In all the
problems discussed below, we assume concentration and temperature to be functions of time t and
radius r (which is a current conﬁguration variable). This assumption is often made because the underlying problem has either cylindrical or spherical symmetry. We also assume that the volumetric
sources corresponding to temperature and concentration are equal to zero. In this paper, as we are
mainly interested in degradation of solid due to temperature and transport of chemical species, we
shall neglect Dufour eﬀect, Soret eﬀect, thermo-chemo coupling parameter dϑc , and anisotropic coeﬃcient of thermal and chemical expansions. In order to reduce the complexity of ﬁnding solutions
based on semi-inverse method for deformation sub-problem, we shall neglect the inertial eﬀects and
body forces.
Based on the assumptions provided here, the governing equations for the transport sub-problem
in cylindrical coordinates reduce to:
ρ

∂c 1 ∂rhr
+
= 0, c(r = ri , t) = ci , c(r = ro , t) = co , c(r, t = 0) = c0
∂t r ∂r
22

(6.1)

where hr is the mass transfer ﬂux in the radial direction. Similarly, the governing equations for the
thermal sub-problem in cylindrical coordinates can be written as:
ρϑ

∂η 1 ∂rqr
∂κ
+
=−
hr , ϑ(r = ri , t) = ϑi , ϑ(r = ro , t) = ϑo , ϑ(r, t = 0) = ϑ0
∂t
r ∂r
∂r

(6.2)

where qr is the heat ﬂux in the radial direction. In spherical coordinates, the governing equations
for the transport sub-problem are:
ρ

∂c
1 ∂r 2 hr
+ 2
= 0, c(r = ri , t) = ci , c(r = ro , t) = co , c(r, t = 0) = c0
∂t r ∂r

(6.3)

The governing equations for the thermal sub-problem in spherical coordinates are:
ρϑ

1 ∂r 2 qr
∂κ
∂η
+ 2
=−
hr , ϑ(r = ri , t) = ϑi , ϑ(r = ro , t) = ϑo , ϑ(r, t = 0) = ϑ0
∂t
r ∂r
∂r

(6.4)

Another quantity of interest in material degradation is the extent of damage at a particular
location or along the cross-section of the degrading body. In case of incompressible neo-Hookean
chemo-thermo-mechano degradation model, this quantity can be deﬁned as follows:

 

µ1 c
µ2 ϑ
µ
=1−
−
(6.5)
Dµ (x, t) :=
µ0
µ0 cref
µ0 ϑref
For virgin material, Dµ = 1. If Dµ approaches zero, then the material has degraded the most. In
addition, equation (6.5) also provides the following information:
◮ Amount of degradation at a given location and time,
◮ The parts of the body that suﬀered extensive damage, and
◮ The eﬀect of temperature and moisture (or concentration of chemical species) on the mechanical properties of materials.
6.1. Inflation of a degrading spherical shell. We now study the behavior of a degrading
(thick) spherical shell subjected to pressure loading. Figure 5 provides a pictorial description of
the boundary value problem. In addition to the obvious theoretical signiﬁcance, this problem has
relevance to safety, reliability and defect monitoring of degrading spherical structures (such as a
tank shell and a bearing structure) due to pressure loading.
Due to the spherical symmetric associated with the problem, spherical coordinates are used to
analyze the inﬂation of degrading spherical shell. Consider a spherical body of inner radius Ri and
outer radius Ro deﬁned in the reference conﬁguration as follows:
Ri ≤ R ≤ Ro ,

0 ≤ Θ ≤ π,

0 ≤ Φ ≤ 2π

(6.6)

where (R, Θ, Φ) are the spherical polar coordinates in the reference conﬁguration. The inner and
outer surfaces R = Ri and R = Ro are, respectively, subjected to pressures pi and po with pi ≥ po .
That is, the thick cylinder is inﬂated with pressure. The deformation in the current conﬁguration
can be described as follows:
ri ≤ r = m(R) ≤ ro ,

θ = Θ,

φ=Φ

(6.7)

where (r, θ, φ) are the spherical polar coordinates in the current conﬁguration, and ri and ro are,
respectively, the inner and outer radii of the shell in the current (deformed) conﬁguration. The
23

deformation gradient, the left Cauchy-Green tensor, and the right Cauchy-Green tensor have the
following matrix representations:



 dm

dm 2
0
0
0
0
dR

 dR
m2
{F} =  0 m
(6.8)
0  , {C} = {B} = 
0
0 
R
R2
2
m
m
0
0 R
0
0 R2
Incompressibility implies that

r = m(R) =
where ro =
Trr

q

q
3

R3 + ri3 − Ri3

ri ≤ r ≤ ro

Ro3 + ri3 − Ri3 . The non-zero components of the Cauchy stress are:


R4
dm 2
r2
= −p + µ(c, ϑ) 4 , Tθθ = Tφφ = −p + µ(c, ϑ) 2
= −p + µ(c, ϑ)
dR
r
R
3

(6.9)

(6.10)

The governing equations for the balance of linear momentum in the spherical polar coordinates
(e.g., see [Sadd, 2014]) reduce to:
2Trr − Tθθ − Tφφ
∂Trr
∂p
∂p
+
= 0,
= 0,
=0
∂r
r
∂θ
∂φ
The above equations imply that p is independent of θ and φ. That is,
p = p(r, t)

(6.11)

(6.12)

From equation (3.6) and (3.7), the speciﬁc chemical potential and speciﬁc entropy for the degrading
spherical shell are given as follows:

 4
1 ∂ψ
µ1
r2
R
κ=
+ Rs ϑref {c − cref } = −
+
2
−
3
+ Rs ϑref {c − cref }
(6.13a)
ρ0 ∂c
2ρ0 cref r 4
R2

 4
cp
cp
R
µ1
r2
1 ∂ψ
+
{ϑ − ϑref } =
+2 2 −3 +
{ϑ − ϑref }
(6.13b)
η=−
4
ρ0 ∂ϑ ϑref
2ρ0 ϑref r
R
ϑref

Before deriving the governing equations for the degrading shell problem, we shall do the nondimensionalization by choosing primary variables and associated reference quantities that are convenient for studying this problem. To distinguish, we shall denote all the non-dimensional quantities
using a superposed bar. We shall take µ0 , Ro , ϑref , cref , and D0 as the reference quantities, which
give rise to the following non-dimensional quantities:
r
R
Dκκ
Dϑϑ
, R=
, D κκ =
, D ϑϑ =
(6.14)
r=
Ro
Ro
D0
D0
µ1
µ2
c
ϑ
D0 t
µ1 =
, µ2 =
, c=
, ϑ=
, t= 2
(6.15)
µ0
µ0
cref
ϑref
Ro
With the stress ﬁeld in equation (6.10), we shall integrate equation (6.11) and then have the following
non-linear equation in deformation sub-problem after non-dimensionalization:

 

6
3
3 2
3
R
o
Z 2µ(c, ϑ) R − R + ri − Ri
T rr (R = Ri , t) − T rr (R = Ro , t) = po − pi =
dR
(6.16)

7
3
3 3
3
R + r i − Ri
Ri

∂ϑ
In order to reduce the complexity in ﬁnding semi-analytical solutions, we shall assume ∂r
∂t ≪ ∂t .
Substituting equation (6.13a) and (6.13b) into the constitutive relations of the proposed model,
24

the governing equations of these two sub-problems (6.3), (6.4) can be written as follows after nondimensionalization:
!
!


4
4
∂c
∂D κκ ∂c
2D κκ
r
R
∂2c
∂D κκ R
1
+
− D κκ 2 = 2ω
(6.17)
−
− 2 − 6ωDκκ
2 + r6
r
∂r
∂r
∂r
∂t
∂r
r5
R
R
!!2


4
D ϑϑ ∂D ϑϑ ∂ϑ
r
R
∂ϑ
∂2ϑ
∂c
(6.18)
− 2
ϑ
+
− D ϑϑ 2 = τ D κκ
− 2ω
2 − r5
r
∂r
∂r
∂r
∂r
∂t
R
where ω and τ are two non-dimensional parameters, which have the following expressions:
ω=

µ1
,
ρ0 Rs ϑref c2ref

τ=

Rs c2ref
cp

(6.19)

The non-linear equation (6.16) enables us to ﬁnd r i at various t for given c(R, t) and ϑ(R, t).
However, it should be noted that c(R, t) and ϑ(R, t) are also a function of r i in case of strong
coupling. This is because diﬀusivity and thermal conductivity depend on the invariants of strain E.
Hence, the integral equation (6.16) and partial diﬀerential equations (6.17) and (6.18) are strongly
coupled.
6.1.1. Steady-state analysis for shell degradation. For steady-state, we have hr r 2 = C1 and
qr r 2 + κhr r 2 = C2 , where C1 and C2 are integration constants. This implies that c and ϑ are the
solutions of the following ODEs:
!
4
3
r
R
dc
+ C1 = 0
D κκ r 2 − 2D κκ ω
(6.20a)
2 − r3
dr
R
!
!
4
2
w
R
r
dϑ
+τ
+ 2 2 − 3 − c + 1 C1 + C2 = 0
D ϑϑ r 2
(6.20b)
dr
2 r4
R
where the integration constants C 1 and C 2 are determined from the boundary conditions for the
transport and thermal sub-problems. Under weak coupling (i.e., Dϑϑ and Dκκ are constants), a
simpliﬁed form of the analytical solutions for c and ϑ can be obtained as follows:
!
4
r2
R
τ B12 D κκ
B1
Z1
+
+
A
,
+ Y1
(6.21)
+
c=ω
ϑ
=
−
+
1
2
4
2
r
r
2r
2D ϑϑ r
R
where A1 , B1 , Y1 , and Z1 are constants, which are given in terms of the boundary conditions ci , co ,
ϑi , and ϑo as follows:
!
4
8R
r2i
B1
A1 = ci −
−ω
(6.22a)
2 + r4
ri
R
i
!!
4
4
8R
8R
r i ro
r2o
r 2i
B1 =
co − ci − ω
(6.22b)
2 + r4 −
2 − r4
ri − r o
R
R
o
i
τ B12 D κκ
Z1
−
2
ri
2D ϑϑ ri



1
ri − r o
τ B12 D κκ 1
Z1 =
ϑo − ϑi −
− 2
r i ro
r2i
ro
2D ϑϑ
Y1 = ϑi +

25

(6.22c)
(6.22d)

6.1.2. Unsteady analysis for shell degradation. Herein, we shall use the method of horizontal
lines [Rothe, 1930; Picard and Leis, 1980] and shooting method [Heath, 2005] to obtain numerical
solutions to equations (6.17) and (6.18). In the method of horizontal lines, the time is discretized
ﬁrst followed by spatial discretization. The time interval of interest [0, I] is divided into N nonI
and tn = n∆t. tn is called the integral time level, where
overlapping subintervals such that ∆t = N
n = 0, · · · , N . ∆t is the time-step, which is assumed to be uniform. Employing the method of
horizontal lines with backward Euler time-stepping scheme, we obtain the following ODEs at each
time-level for equations (6.17) and (6.18):

 

!
!
(n) 4
R
d2 c(n+1)
1
1
dc(n+1)
c(n+1)
2
dD κκ


+
+
−
=
6ω
+



6 
2
2
(n)
(n)
(n)
(n)
dr
dr
dr
r (n)
r
D κκ
Dκκ ∆t
R
t=tn


!  (n) 4
!
(n)
(n)
R
c
r
2ω
dD κκ


−
−
− (n)
(6.23)

2 

5
(n)
(n)
(n)
dr
r
D κκ ∆t
Dκκ
R
t=tn
(n+1)
d2 ϑ

dr 2

+

−

2
r (n)
(n)
τ D κκ
(n)
D ϑϑ

1

+


(n)
Dϑϑ

 dc

dr
(n)

!

dD ϑϑ
dr
t=tn


!

−

dr

ϑ

(n) (n+1)

ϑ

(n)
D ϑϑ ∆t
4 2


(n)
R


− 2ω  
 −
5 
(n)
(n) 2
r
R
r (n)

t=tn

dϑ

(n+1)



=−



(n) 2
ϑ
(n)

Dϑϑ ∆t

(6.24)

where c(n) = c(r, t = tn ) and ϑ = ϑ(r, t = tn ). Algorithm 1 describes a procedure to determine
c(r, t), ϑ(r, t), and ri at various times using an iterative non-linear numerical solution strategy. The
following values are assumed for the non-dimensional parameters in the strong coupling simulations:
Ro = 1, Ri = 0.5, ∆t = 0.01, t = 2, ω = 0.05, τ = 0.2, ci = 0, ϑ0 = 0.5
co = 1, ϑi = 0.5, ϑo = 1, µ0 = 1, µ1 = 0.3, µ2 = 0.4, D 0 = 1, D T = 1.5,
D S = 1.2, ηT = ηS = 1, ErefT = ErefS = 1, K 0 = 1, δ = 10
(n)

(6.25)

(n)

In weakly coupling problem, we use D 0 , K 0 as D κκ and D ϑϑ , respectively. It should be noted that
these values are constructed based on the (brittle-type) material parameters such as glass, ceramics,
and concrete.
The numerical results are shown in ﬁgures 7–11, which reveal the following conclusions on the
overall behavior of degrading spherical shells under inﬂation:
(i) Degradation vs. non-degradation: After degradation, a spherical shell which is initially homogeneous is not homogeneous anymore.
(ii) Due to degradation, creep-like behavior is observed. Therefore, as time progresses, hoop
stresses increase. We need to note that the shell ceases to creep after a certain period of time,
which is the moment when the transport of chemical species and heat conduction are close to
steady-states.
(iii) As the pressure loading increases, the hoop stress is increasing in a non-linear fashion, which
is signiﬁcantly diﬀerent from the non-degradation shell.
26

(iv) For non-degrading shell, the chemical potential is unchanged with respect to pressure loading.
However, for strong coupling, it increases with pi in a non-linear fashion when ω is small
enough. This is because for small ω, diﬀusion takes the dominance in the coupling eﬀect.
When pressure loading increases, the diﬀusivity is increasing due to the growing strain. For
large ω, the deformation is dominant in the coupling, which is −I E term in chemical potential.
Since the ﬁrst invariant I E is always positive in this problem, chemical potential is decreasing
when the pressure loading increases.
(v) Thermo-dominated vs. chemo-dominated degradation: Weak coupling over-predicts the amount
of degradation compared to the full (or strong) coupling when thermal degradation dominates. This is because when the thermal degradation dominants, the thermal conductivity
decreases due to the increase in strain (note that the ﬁrst invariant of strain is always positive
in this problem). However, in chemo-dominated degradation, weak coupling under-predicts
the amount of degradation compared to the strong coupling case.
(vi) In case of strong coupling, healing-like behavior is observed at early time steps in thermodominated degradation (but still remains below that of the virgin material). This is because
variable heat sinks exist in the entire body (due to which temperature gets lower than the
initial condition). Hence, the material damage is less than that of at t = 0. However, this
heal-like behavior becomes less distinct (or even doesn’t exist) when the chemo-degradation
achieves the dominance.
(vii) Strong vs. weak coupling: Quantitatively and qualitatively, extent of damage for both strong
and weak coupling are considerably diﬀerent.
6.2. Bending of a degrading beam. Herein, we shall consider pure bending of a degrading
beam. At time t = 0, a ﬁnite degrading beam is suddenly bent by an action of pure end moments.
For t > 0, the centerline of the beam becomes a sector of a circle of radius rc . This centerline
is held ﬁxed for all the time. Subsequently, the stresses in the degrading beam are allowed to
relax. In addition, it is assumed that the material remains isotropic with respect to the reference
conﬁguration throughout the degradation process. These assumptions enable us to employ the
counterpart of universal deformations (also known as semi-inverse method) [Ogden, 1997] to study
such degrading beams.
A pictorial description of the initial boundary value problem is shown in Figure 12. The degrading beam is deﬁned as follows:
−L ≤ X ≤ L,

−W ≤ Y ≤ W,

−H ≤ Z ≤ H

(6.26)

where (X, Y, Z) are the Cartesian coordinates in the reference conﬁguration. We assume that the
deformation to be as follows:
r
Y
2X
+ β, θ = , z = Z
(6.27)
r=
α
γ
where (r, θ, z) are the cylindrical polar coordinates in the current conﬁguration. When X = 0,
we have β = rc2 . It should be noted that α and γ are all unknown time-dependent parameters.
These unknowns are evaluated from the incompressibility constraint, traction boundary conditions,
and pure end moments. To reduce the complexity in ﬁnding semi-analytical solutions, we shall
assume rc is given. The faces X = −L and X = L are subjected to ambient atmospheric pressure
‘patm ’. Upon deformation, the corresponding deformed faces ri and ro are maintained at patm , where
27

Algorithm 1 Inﬂation of a degrading spherical shell (numerical methodology to ﬁnd r i , c, and ϑ)
1:
2:
3:
4:
5:
6:
7:
8:
9:

INPUT: Non-dimensional material parameters, non-dimensional boundary conditions, and non(r) (c)
(ϑ)
dimensional initial conditions, MaxIters, tolerances ǫtol , ǫtol , and ǫtol .
Evaluate r i at t = 0 based on equation (6.16).
for n = 1, 2, · · · , N do
for j = 1, 2, · · · do
if j > MaxIters then
Solution did not converge in speciﬁed maximum number of iterations. EXIT.
end if
(j)
Diffusion sub-problem: Given ri , solve equation (6.23) to obtain c(j+1) . Herein, we use
shooting method to solve the ODEs.
(j)
Heat conduction sub-problem: Given r i and c(j+1) , solve equation (6.24) to obtain
(j+1)

10:
11:
12:
13:
14:
15:
16:
17:

ϑ
. Similarly, we use shooting method to solve the non-linear ODEs.
(j+1)
(j+1)
, solve for ri
given by equation
Deformation sub-problem: Given c(j+1) and ϑ
(6.16) using bisection method.
(j+1)
(j)
(j+1)
(j)
(r)
(c)
(ϑ)
− r i k < ǫtol , kc(j+1) − c(j) k < ǫtol , and kϑ
− ϑ k < ǫtol then
if kr i
(j+1)

(j+1)

OUTPUT: ri
, c(j+1) , and ϑ
else
(j)
(j+1)
.
Update the guess: r i ← r i
end if
end for
end for

. EXIT.

p
p
ri = rc2 − 2γL and ro = rc2 + 2γL are the inner and outer radius of the degrading beam. This
gives the following traction boundary conditions:
Trr (X = −L, t) = Trr (X = L, t) = patm

(6.28)

The deformation gradient F, right Cauchy-Green tensor C, and left Cauchy-Green tensor B for the
degrading beam are given as follows:




1
1
0
0
0
0
2
2
 αr


 αr
r2
0 
{F} =  0 γr 0 
(6.29)
{C} = {B} =  0
γ2
0 0 1
0
0 1
For incompressible degrading neo-Hookean material, we have αγ = 1 and the non-zero components
of the Cauchy stress tensor are given as follows:

µ(c, ϑ) 2γX + rc2
µ(c, ϑ)γ 2
, Tθθ = −p +
, Tzz = −p + µ(c, ϑ)
(6.30)
Trr = −p +
2γX + rc2
γ2
The balance of linear momentum in the cylindrical polar coordinates reduces to the following:
∂Trr
Trr − Tθθ
+
= 0,
∂r
r
28

∂p
= 0,
∂θ

∂p
=0
∂z

(6.31)

The bending moment in the deformation sub-problem can be evaluated based on the following
formula:
Z
Mbeam (t) =
Tθθ (r − rneu )dA
Across

= 2H

ZL

Tθθ (−

−L

p

rc2 + 2γXneu +

p
rc2 + 2γX) p

γ
rc2

+ 2γX

dX

(6.32)

p
where dA = 2Hdr, rneu = rc2 + 2γXneu is the neutral axis location, and Xneu is the value at which
Tθθ = 0. The chemical potential, speciﬁc entropy for the degrading beam are given as follows:

 2
1 ∂ψ
γ
µ1
r2
κ=
+ Rs ϑref {c − cref } = −
+ 2 − 2 + Rs ϑref {c − cref }
(6.33a)
ρ0 ∂c
2ρ0 cref r 2
γ

 2
1 ∂ψ
cp
µ1
r2
cp
γ
η=−
+
{ϑ − ϑref } =
+ 2 −2 +
{ϑ − ϑref }
(6.33b)
2
ρ0 ∂ϑ ϑref
2ρ0 ϑref r
γ
ϑref
Most of the non-dimensional quantities are same as that of the degrading shell problem except for
the following:
r=

r
X
γ
D0 t
, X= , γ= , t= 2
rc
rc
rc
rc

Using equations (6.27)–(6.31), we have the following non-linear equation in γ


L/r
Z c µ(c(X, t), ϑ(X, t)) γ 4 − 2γX + 12
dX = 0
2
γ 2γX + 1

(6.34)

(6.35)

−L/rc

From (6.35), γ|t=0 is given as follows:

γ|t=0

1
=
rc

q

−2L2 +

p

4L4 + rc4

(6.36)

which is the case for homogeneous neo-Hookean material. As rc is given, the parameter γ is bounded
above and below as follows:
−rc
rc
<γ<
2L
2L

(6.37)

which can be used in ﬁnding the solution for the non-linear equation given by (6.35). It should be
noted that γ|t=0 satisﬁes the inequality given by (6.37).
From equations (6.1) and (6.2), the ﬁnal form for the governing equations for transport and
thermal sub-problems for degrading beam are given as follows:






D κκ
r
γ2
∂2c
∂D κκ γ 2
1
∂c
∂D κκ ∂c
−
− 2 − 2ωD κκ
+ 4
(6.38)
+
− D κκ 2 = ω
r
∂r
∂r
∂r
∂t
∂r
r3
γ
γ2
r

 2

2

∂ϑ
∂2ϑ
∂c
D ϑϑ ∂D ϑϑ ∂ϑ
γ
r
+
− Dϑϑ 2 = τ D κκ
+ω
ϑ
−
− 2
(6.39)
r
∂r
∂r
∂r
∂t
∂r
r3
γ
29

6.2.1. Steady-state and unsteady analysis for beam degradation. In case of steady-state, we have
hr r = C1 and qr r + κhr r = C2 , where C1 and C2 are integration constants. Equations (6.38) and
(6.39) imply that c and ϑ are the solutions of the following ODEs:

 2
γ
r2
dc
− 2 + C1 = 0
D κκ r − D κκ ω
(6.40a)
dr
r2
γ

  2

w γ
r2
dϑ
+τ
+ 2 − 2 − c + 1 C1 + C2 = 0
D ϑϑ r
(6.40b)
dr
2 r2
γ

In case of weak coupling (where Dϑϑ and Dκκ are constants), the solutions for c and ϑ take the
following simpliﬁed form:


ω γ2
r2
τ B22 D κκ
+
+
B
ln[r]
c=−
+
A
,
ϑ
=
−
ln[r]2 + Z2 ln[r] + Y2
(6.41)
2
2
2
2 r2
γ
2D ϑϑ
where the constants A2 , B2 , Y2 , and Z2 (which depend on the boundary conditions) are as follows:


ω γ2
r2i
A2 = ci − B2 ln[r i ] +
+ 2
(6.42a)
2 r 2i
γ



1
ω γ2
r 2i
γ2
r 2o
+ 2− 2 − 2
(6.42b)
co − ci −
B2 =
ln[ro ] − ln[r i ]
2 r 2i
γ
ro
γ
τ B22 Dκκ
ln[ri ]2 − Z2 ln[ri ]
2D ϑϑ



τ B22 D κκ
1
2
2
ln[r i ] − ln[r o ]
Z2 =
ϑo − ϑi −
ln[ro ] − ln[r i ]
2D ϑϑ

(6.42c)

Y2 = ϑi +

(6.42d)

For unsteady analysis, we employ method of horizontal lines with backward Euler time-stepping
scheme. This gives the following ODEs at each time-level for equations (6.38) and (6.39):
!
!
2 !
γ (n)
d2 c(n+1)
c(n+1)
dc(n+1)
1
dD κκ
1
1
− (n)
= 2ω
+
+
2 + (n) 4
(n)
dr
dr
dr2
r(n)
γ (n)
r
D κκ
D κκ ∆t
t=tn
!
!
!

2
γ (n)
c(n)
ω
r(n)
dD κκ
− (n)
(6.43)
−
−

2
3
(n)
dr
r (n)
γ (n)
D ∆t
D
κκ

(n+1)
d2 ϑ

dr 2

+

1

+

r(n)

(n)

−

2τ ωD κκ
(n)

D ϑϑ

t=tn

κκ

1
(n)
D ϑϑ

!

dD ϑϑ
dr

t=tn

!

dϑ

(n+1)

−

dr

!
2
γ (n)
r (n)
dc
−


3
2
dr
r(n)
γ (n)

ϑ

(n) (n+1)

ϑ

(n)
Dϑϑ ∆t
(n)

−
t=tn

τ D κκ ω 2
(n)

D ϑϑ

=

(n)
τ D κκ
− (n)
D ϑϑ

dc
dr

!2

t=tn
 !2

−



(n) 2
ϑ
(n)

D ϑϑ ∆t

2
γ (n)
r(n)
−

2
3
r (n)
γ (n)

(6.44)

Algorithm 2 describes a procedure to determine c(r, t), ϑ(r, t), and γ at various times using an
iterative non-linear numerical solution strategy. The boundary conditions for diﬀusion and thermal
subproblems are the same as the degrading shell problem. The other parameters are assumed in
the strongly coupling simulations as follows:
L = 1, rc = 1, ∆t = 0.1, t = 2, ω = 0.05, τ = 0.5, µ0 = 1, µ1 = µ2 = 0.4, D 0 = 1,
D T = 2.0, DS = 1.5, ηT = ηS = 1, ErefT = ErefS = 1, K 0 = 1, δ = 10
30

(6.45)

(n)

(n)

In case of weak coupling, we have D0 as D κκ and K 0 as D ϑϑ , respectively.
The numerical results are shown in ﬁgures 13–16, which reveal the following conclusions on the
overall behavior of bending of degrading beams:
(i) Degradation vs. non-degradation: The main observation is that the neutral axis shifts further
to the left, similar to the phenomenon observed in viscoelastic solids [Kolberg and Wineman,
1997]. Moreover, in case of weak coupling for some instants of time the maximum stress does
not occur at either tensile or compressive sides of the beam after the onset of degradation. This
is of primal importance in regards to the calculation of failure loads/moments due to material
damage. Hence, a simple approach based on strength of materials or a more complex ﬁnite
elasticity theory to calculate stresses without accounting for degradation will lead to erroneous
results.
(ii) Initially at t = 0 and when there is no degradation, the response is that of a homogeneous
neo-Hookean material. On the onset of degradation, the material ceases to be homogeneous.
(iii) Moment relaxation is observed for weak and strong coupling degradation. Note that the
moment is a constant without degradation. Moreover, although diﬀusion is dominant in the
coupling eﬀect for chemical potential, one can still observe the deformation eﬀect on κ as
compared with no degradation case.
(iv) Strong vs. weak coupling: One can see that T θθ for strong coupling is considerably diﬀerent from
the weak coupling. This is because the degradation progress is dependent on the deformation,
concentration of the diﬀusing chemical species, and temperature of the body.
(v) The extent of damage is monotonic for weak coupling, which is not the case for strong coupling
(which helps in identifying regions that need retroﬁtting).
Remark 6.1. Based on a semi-inverse approach, under degradation, [Rajagopal et al., 2007]
have shown that there is a shift in the neutral axis for pure bending of a polymer beam. However,
their model is based on internal variables, which is diﬃcult to calibrate experimentally. On the other
hand, the proposed (and calibrated) chemo-thermo-mechano degradation model is able to predict the
shift of neutral axis without appealing to internal variable framework.
6.3. Torsional shear of a degrading cylinder. A pictorial description of the degrading
cylindrical annulus of ﬁnite length is shown in Figure 17. The bottom of the cylinder is ﬁxed
and just after time t = 0, a twisting moment is applied. We analyze the material degradation and
corresponding structural response due to the torsional shear for a prescribed angle of twist. Initially,
the body is a homogeneous neo-Hookean material and there is no transport of chemical species in
the body. For time t > 0, the outer boundary of the cylinder is always exposed to moisture (or
a diﬀusing chemical species). The inner surface of the degrading annular cylinder is held at zero
concentration. This can be achieved by constructing a mechanism which continuously removes the
moisture (or diﬀusing chemical species) from the inner boundary of the degrading cylinder. Hence,
one can control the concentration of the moisture at both inner and outer surfaces. Similar type of
initial and boundary conditions are enforced for the thermal counter part.
Consider a closed cylindrical body of inner radius Ri , outer radius Ro , and height L deﬁned as
follows:
Ri ≤ R ≤ Ro ,

0 ≤ Θ ≤ 2π,
31

0≤Z ≤L

(6.46)

Algorithm 2 Pure bending of degrading beam (numerical methodology to ﬁnd γ, c, and ϑ)
1:
2:
3:
4:
5:
6:
7:
8:
9:

INPUT: Non-dimensional material parameters, non-dimensional boundary conditions, and non(γ) (c)
(ϑ)
dimensional initial conditions, MaxIters, tolerances ǫtol , ǫtol , and ǫtol .
Evaluate γ at t = 0 based on equation (6.36). Use this as an initial guess for solving nonlinear
equation given by (6.35) or guess γ based on equation (6.37).
for n = 1, 2, · · · , N do
for i = 1, 2, · · · do
if i > MaxIters then
Solution did not converge in speciﬁed maximum number of iterations. EXIT.
end if
Diffusion sub-problem: Given γ (i) , solve equation (6.43) to obtain c(i+1) . Herein, we use
shooting method to solve the ODEs.
Heat conduction sub-problem: Given γ (i) and c(i+1) , solve equation (6.44) to obtain
(i+1)

10:
11:
12:
13:
14:
15:
16:
17:

ϑ
. Similarly, we use shooting method to solve the ODEs.
(i+1)
Deformation sub-problem: Given c(i+1) and ϑ
, solve for γ (i+1) given by equation
(6.35) using bisection method.
(i+1)
(i)
(γ)
(c)
(ϑ)
if kγ (i+1) − γ (i) k < ǫtol , kc(i+1) − c(i) k < ǫtol , and kϑ
− ϑ k < ǫtol then
(i+1)

OUTPUT: γ (i+1) , c(i+1) , and ϑ
else
Update the guess: γ (i) ← γ (i+1) .
end if
end for
end for

. EXIT the inner loop.

where (R, Θ, Z) are the cylindrical polar coordinates in the reference conﬁguration. Under torsional
shear, the deformation can be described as follows:
r = R,

θ = Θ + g(Z, t),

z = ΛZ

The components of the deformation gradient F can be written as:


1 0 0
∂g(Z, t)
{F} =  0 1 rg′  where g′ :=
∂Z
0 0 Λ

(6.47)

(6.48)

Incompressibility implied that Λ = 1. The components of the right Cauchy-Green tensor C and the
left Cauchy-Green tensor B can be written as:




1 0
0
1
0
0
 {B} =  0 1 + (rg′ )2 rg′ 
rg ′
{C} =  0 1
(6.49)
2
′
′
′
0 rg 1 + (rg )
0
rg
1
The non-zero components of the Cauchy stress T are given as follows:

2 
Trr = −p + µ(c, ϑ), Tθθ = −p + µ(c, ϑ) 1 + rg′
Tzz = −p + µ(c, ϑ),

Tθz = Tzθ = µ(c, ϑ)rg′
32

(6.50)

The balance of linear momentum in the cylindrical polar coordinates reduces to the following:
−

2
∂p
+ µ(c, ϑ)r g′ = 0,
∂r
∂p
∂θ

Symmetry in the problem implies that
takes the following form:

−

1 ∂p
+ µ(c, ϑ)rg′′ = 0,
r ∂θ

−

∂p
=0
∂z

(6.51)

= 0, which further implies that g ′′ = 0. Hence, g(Z, t)
(6.52)

g(Z, t) = Ψ1 (t)Z + Ψ2 (t)

where Ψ1 and Ψ2 are evaluated based on the input data. As the bottom of the cylinder is ﬁxed, we
have g(Z = 0, t) = 0, which implies Ψ2 (t) = 0.
The chemical potential and speciﬁc entropy are given as follows:
1 ∂ψ
µ1 r 2 Ψ21
+ Rs ϑref {c − cref } = −
+ Rs ϑref {c − cref }
ρ0 ∂c
2ρ0 cref
1 ∂ψ
cp
µ1 r 2 Ψ21
cp
η=−
+
{ϑ − ϑref } =
+
{ϑ − ϑref }
ρ0 ∂ϑ ϑref
2ρ0 ϑref
ϑref

κ=

(6.53a)
(6.53b)

Most of the non-dimensional quantities remain the same as that of the previous initial boundary
value problems except for the following:
R=

R
D0 t
, ψ = ψRo , t = 2
Ro
Ro

(6.54)

The non-dimensional twisting moment M (t) satisﬁes:

M (t) = 2π

ZRo

3

µ(c(R, t), ϑ(R, t))Ψ1 R dR

(6.55)

Ri

The Poynting eﬀect for hyperelastic materials shall be also studied. It implies the axial length
change for a cylinder under shear. The non-dimensional normal force required to keep the length
unchanged can be written as follows:
N (t) = π

ZRo

2

3

µ(c(R, t), ϑ(R, t))Ψ1 R dR

(6.56)

Ri

From equations (6.1) and (6.2), the ﬁnal form of the governing equations for transport and thermal
sub-problems can be written as:




∂2c
∂D κκ
∂D κκ ∂c
∂c
D κκ
2
+
− D κκ 2 = −ωΨ1 2D κκ + r
(6.57)
−
r
∂r
∂r
∂r
∂t
∂r


 2
D ϑϑ ∂D ϑϑ ∂ϑ
∂ϑ
∂2ϑ
∂c
2 ∂c
4
− 2τ ωD κκ rΨ1
+
− Dϑϑ 2 = τ D κκ
+ τ D κκ ω 2 r 2 Ψ1 (6.58)
−
ϑ
r
∂r
∂r
∂r
∂r
∂t
∂r
One needs to solve equations (6.56)–(6.58) to obtain c(r, t), ϑ(r, t), and M (t). Algorithm 3 describes
a numerical solution procedure to solve these equations at various times for a given angle of twist
per unit length.
33

6.3.1. Steady-state and unsteady response of degrading cylinder under torsional shear. In the
case of steady-state, c and ϑ are the solutions of the following ODEs:
dc
2
− D κκ ωrΨ1 + C 1 = 0
dr


ω 2 2
dϑ
D ϑϑ r
r Ψ1 − c + 1 C 1 + C 2 = 0
+τ
dr
2
D κκ r2

(6.59a)
(6.59b)

where C 1 and C 2 are integration constants. Under weak coupling (where Dϑϑ and Dκκ are constants), a simpliﬁed form of the analytical solutions for c and ϑ is given as follows:
c=

ω 2 2
r Ψ1 + B3 ln[r] + A3 ,
2

ϑ=−

τ B32 D κκ
ln[r]2 + Z3 ln[r] + Y3
2D ϑϑ

(6.60)

where A3 , B3 , Y3 , and Z3 are constants, which are obtained by the corresponding boundary conditions for thermal and diﬀusion sub-problem. These are given as follows:
ω
2
A3 = ci − B3 ln[r i ] − r2i Ψ1
2

1
ω 2 2
2 2
B3 =
co − ci −
r Ψ − ri Ψ1
ln[ro ] − ln[r i ]
2 o 1

(6.61a)
(6.61b)

τ B32 Dκκ
ln[ri ]2 − Z3 ln[ri ]
2D ϑϑ



τ B32 D κκ
1
2
2
ln[r i ] − ln[r o ]
Z3 =
ϑo − ϑi −
ln[ro ] − ln[r i ]
2D ϑϑ

(6.61c)

Y3 = ϑi +

(6.61d)

For unsteady analysis, method of horizontal lines with backward Euler time-stepping scheme is
employed. This gives the following ODEs at each time-level:
!
!


c(n+1)
dc(n+1)
1
dD κκ
1
d2 c(n+1)
(n) 2
−
=
2ω
+
Ψ
+
1
(n)
(n)
dr
dr
dr2
r(n)
D κκ
D κκ ∆t
t=tn
!
 (n) dD

c(n)
(n) 2 r
κκ
+ ω Ψ1
−
(6.62)
(n)
(n)
dr
D
D ∆t
t=tn

κκ

(n+1)
d2 ϑ

dr 2

+

1

+

1

!

dD ϑϑ
dr

!

dϑ

(n+1)

−

κκ

ϑ

(n) (n+1)

ϑ

=

(n)
τ D κκ
− (n)
D ϑϑ

(n)
(n)
dr
D ϑϑ
Dϑϑ ∆t
t=tn
(n)
(n)
τ Dκκ ω 2  (n) 2  (n) 4
2τ ωD κκ (n)  (n) 2 dc
+
r
Ψ
r
Ψ1
−
1
(n)
(n)
dr
D ϑϑ
D
ϑϑ
t=tn

r(n)

dc
dr

t=tn

!2

−



(n) 2
ϑ
(n)

D ϑϑ ∆t

(6.63)

The boundary conditions for diﬀusion and thermal sub-problems are the same as that of the previous
boundary value problems.
The following non-dimensional parameters are assumed in the numerical numerical simulations:
Ro = 1, Ri = 0.5, ∆t = 0.1, t = 2, ω = 0.05, τ = 0.8, µ0 = 1, µ1 = 0.5, µ2 = 0.2,
D 0 = 1, D T = 1.5, D S = 1.2, ηT = ηS = 0.1, ErefT = ErefS = 1, K 0 = 1, δ = 10

(6.64)

The numerical results are shown in Figure 18 and 19, which reveals the following important conclusions on the overall behavior of degrading structural members under torsional shear:
34

(i) The numerical results reveal that there is relaxation of moment for ﬁxed deformation. In
addition, the twisting moment required to maintain a ﬁxed angle of twist decreases with
increase in µ1 . Similar type of behavior is observed when µ1 is kept constant and µ2 is varied.
(ii) We observe moment relaxation due to material degradation when both the transport and
thermal sub-problems are close to steady states. Moreover, one can see that moment relaxation
depends on the geometry of the specimen. These aspects diﬀerentiate the stress relaxation
due to degradation from the stress relaxation due to viscoelasticity.
(iii) We observe that the normal force due to Poynting eﬀect is decreasing over time as a result
of degradation. Without degradation, the normal force is a constant (which is the case for
hyperelastic materials).
Algorithm 3 Torsional shear of a degrading cylinder (numerical methodology to ﬁnd M , c, and ϑ)
1:
2:
3:
4:
5:
6:

INPUT: Non-dimensional material parameters, non-dimensional boundary conditions, and nondimensional initial conditions.
for n = 1, 2, · · · , N do
Diffusion sub-problem: Given Ψ1 , solve equation (6.62) to obtain c(n) . Herein, we use
shooting method to solve the ODEs.
(n)
Heat conduction sub-problem: Given Ψ1 and c(n) , solve equation (6.63) to obtain ϑ .
Similar to diﬀusion sub-problem, we use shooting method to solve the non-linear ODEs.
(n)
(n)
Deformation sub-problem: Given c(n) and ϑ , solve for M
given by equation (6.56).
end for

7. CONCLUDING REMARKS
This paper has made several contributions to the modeling of degradation of materials due to
the presence of an adverse chemical species and temperature. First, a consistent mathematical
model has been derived that has ﬁrm continuum thermodynamics underpinning. The constitutive
relations, which give rise to coupled deformation-thermal-transport equations, have been derived by
appealing to the maximization of the rate of dissipation, which is a stronger version of the second
law of thermodynamics. The proposed model is hierarchical in the sense that it recovers many
existing models as special cases. Second, the materials parameters have been calibrated with an
experimental dataset available in the literature. Third, it has been shown that the unsteady solutions
to the proposed degradation model are bounded and stable in the sense of Lyapunov even under
large deformations and large strains. Last but not the least, using several canonical problems in
degradation mechanics, we illustrated the eﬀects of chemical degradation and thermal degradation
on the response of a body that is initially hyperelastic. Some of the main features of degradation
and of the proposed model can be summarized as follows:
(C1) Degradation introduces spatial inhomogeneity. That is, a material which is originally homogeneous may cease to be homogeneous due to degradation.
(C2) The proposed mathematical model can provide the variation of important quantities like chemical potential within the body, which is essential in incorporating chemical reactions into the
modeling.
35

(C3) The extent of damage in a structural member can be both qualitatively and quantitatively
diﬀerent under strong and weak couplings between mechanical, thermal and transport processes. More importantly, weak coupling may over-predict the material degradation in some
cases while in other cases it may under-predict the degradation. It is, therefore, of paramount
importance to select the extent of coupling between the mechanical, thermal and chemical
processes.
(C4) The usual assumptions on either kinematics or stresses, which may be justiﬁed for nondegrading members, may no longer hold under degradation. For example, assumptions on
the location of neutral axis or the location of the maximum stress on the outer ﬁbers in beam
bending will not hold under degradation.
(C5) Degrading structural members may exhibit some responses that are typically associated with
viscoelasticity. In particular, we have shown that degradation can induce stress relaxation
and creep in the response of the materials even in the case of ﬁnite-sized bodies. In contrast
to a viscoelastic body (which creeps continuously upon the application of a load) the body
undergoing chemical degradation ceases to creep for practical purposes after a certain period
of time. This the moment when the transport of chemical species is close to a steady-state, if
there is no volumetric source and the boundary conditions are unchanged over time. A similar
trend holds even in the case of thermal degradation. This characteristic behavior of degrading
solids can be used to diﬀerentiate the creep associated with viscoelasticity and degradation.
Moreover, stress relaxation due to degradation depends on the geometry of the specimen,
which is also diﬀerent from the case due to viscoelasticity.
A possible future research work can be towards incorporating fatigue and fracture into the
degradation modeling. A related scientiﬁc question can be towards addressing the eﬀect of material
degradation on the crack initiation and its propagation.
ACKNOWLEDGMENTS
The authors acknowledge the support from the Department of Energy through Nuclear Energy
University Programs (NEUP). The opinions expressed in this paper are those of the authors and
do not necessarily reﬂect that of the sponsor(s).
References
ABAQUS/CAE/Standard, Version 6.14-1.
Simulia, Providence, Rhode Island, USA,
www.simulia.com, 2014.
COMSOL Multiphysics User’s Guide, Version 5.0-1. COMSOL, Inc., Burlington, Massachusetts,
USA, 2014.
ANSYS Multiphysics, Version 16.0.
ANSYS, Inc., Canonsburg, Pennsylvania, USA,
www.ansys.com, 2015.
S. B. Adler. Chemical expansivity of electrochemical ceramics. Journal of the American Ceramic
Society, 84:2117–2119, 2001.
S. M. Allam, H. M. F. Elbakry, and A. G. Rabeai. Behavior of one-way reinforced concrete slabs
subjected to ﬁre. Alexandria Engineering Journal, 52:749–761, 2013.
S. S. Antman. Nonlinear Problems of Elasticity. Springer-Verlag, New York, 1995.
A. W. Batchelor, L. N. Lam, and M. Chandrasekaran. Materials Degradation and Its Control by
Surface Engineering. Imperial College Press, London, third edition, 2003.
36

S. Bhowmick and V. B. Shenoy. Eﬀect of strain on the thermal conductivity of solids. The Journal
of Chemical Physics, 125:164513, 2006.
F. Björk, C. A. Eriksson, S. Karlsson, and F. Khabbaz. Degradation of components in ﬂooring
systems in humid and alkaline environments. Construction and Building Materials, 17:213–221,
2003.
E. Blond and N. Richet. Thermomechanical modelling of ion–conducting membrane for oxygen
separation. Journal of the European Ceramic Society, 28:793–801, 2008.
H. Bouadi and C. T. Sun. Hygrothermal eﬀects on the stress ﬁeld of laminated composites. Journal
of Reinforced Plastics and Composites, 8:40–54, 1989.
H. Bouadi and C. T. Sun. Hygrothermal eﬀects on structural stiﬀness and structural damping of
laminated composites. Journal of Materials Science, 25:499–505, 1990.
R. M. Bowen. Theory of mixtures. In A. C. Eringen, editor, Continuum Physics, volume III.
Academic Press, New York, 1976.
M. Buonsanti, G. Leonard, and F. Scoppelliti. Equilibrium state of a binary granular solids mixture.
Applied Mechanics and Materials, 52:389–392, 2011.
L. W. Cai and Y. J. Weitsman. Non–Fickian moisture diﬀusion in polymeric composites. Journal
of Composite Materials, 28:130–154, 1994.
D. W. Cho and K. Kim. The mechanisms of moisture damage in asphalt pavement by applying
chemistry aspects. KSCE Journal of Civil Engineering, 14:333–341, 2010.
B. D. Coleman and E. H. Dill. On thermodynamics and the stability of motions of materials with
memory. Archive for Rational Mechanics and Analysis, 51:1–53, 1973.
O. Coussy. Poromechanics. John Wiley & Sons, Inc., New York, 2004.
J. C. Criscione, J. D. Humphrey, A. S. Douglas, and W. C. Hunter. An invariant basis for natural
strain which yields orthogonal stress response terms in isotropic hyperelasticity. Journal of the
Mechanics and Physics of Solids, 48:2445–2465, 2000.
S. Darbha and K. R. Rajagopal. Unsteady motions of degrading or aging linearized elastic solids.
International Journal of Non-Linear Mechanics, 44:478–485, 2009.
C L. Dym. Stability Theory and Its Applications to Structural Mechanics. Dover Publications, New
York, 2002.
J. L. Ericksen. Thermoelastic stability. In Proceedings 5th US National Congress of Applied Mechanics, pages 187–193, 1966.
D. Gawin, F. Pesavento, and B. A. Schreﬂer. Modeling deterioration of cementitious materials exposed to calcium leaching in non-isothermal conditions. Computer Methods in Applied Mechanics
and Engineering, 198:3051–3083, 2009.
F. P. Glasser, J. Marchand, and E. Samson. Durability of concrete degradation phenomena involving
detrimental chemical reactions. Cement and Concrete Research, 38:226–246, 2008.
S. Grasberger and G. Meschke. Thermo-hygro-mechanical degradation of concrete: From coupled 3D
material modeling to durability–oriented multiﬁeld structural analyses. Materials and Structures,
37:244–256, 2004.
J. D. Gu, T. E. Ford, N. S. Berke, and R. Mitchell. Biodeterioration of concrete by the fungus
Fusarium. International Biodeterioration & Biodegradation, 41:101–109, 1998.
M. E. Gurtin. Thermodynamics and stability. Archives of Rational Mechanics and Analysis, 59:
53–96, 1975.
J. K. Hale and H. Kocak. Dynamics and Bifurcations. Springer-Verlag, New York, 1991.
37

C. A. Harper. Handbook of Plastics, Elastomers, & Composites. McGraw–Hill, New York, fourth
edition, 2002.
M. T. Heath. Scientiﬁc Computing–An Introductory Survey. McGraw-Hill, New York, USA, second
edition, 2005.
A. W. Herrmann. ASCE 2013 Report Card for America’s Infrastructure. In IABSE Symposium
Report, volume 99, pages 9–10. International Association for Bridge and Structural Engineering,
2013.
G. A. Holzapfel. Nonlinear Solid Mechanics. John Wiley & Sons, Inc., Chichester, 2000.
E. Jarkova, H. Pleiner, H. W. Müller, A. Fink, and H. R. Brand. Hydrodynamics of nematic
ferroﬂuids. The European Physical Journal E, 5:583–588, 2001.
Y. G. Jung, I. M. Peterson, D. K. Kim, and B. R. Lawn. Lifetime-limiting strength degradation
from contact fatigue in dental ceramics. Journal of Dental Research, 79:722–731, 2000.
L. Kachanov. Introduction to Continuum Damage Mechanics. Springer Science & Business Media,
Dordrecht, Netherlands, 1986.
M. F. Kaplan. Concrete Radiation Shielding: Nuclear Physics, Concrete Properties, Design and
Construction. John Wiley & Sons, Inc., New York, 1989.
S. Karra and K. R. Rajagopal. Degradation and healing in a generalized neo-hookean solid due to
infusion of a ﬂuid. Mechanics of Time-Dependent Materials, 16:85–104, 2012.
D. Klepach and T. I. Zohdi. Strain assisted diﬀusion: Modeling and simulation of deformationdependent diﬀusion in composite media. Composites Part B: Engineering, 56:413–423, 2014.
R. Kolberg and A. Wineman. Response of beams of non-linear viscoelastic materials exhibiting
strain-dependent stress relaxation. International Journal of Non-linear Mechanics, 32:863–883,
1997.
W. M. Lai, J. S. Hou, and V. C. Mow. A triphasic theory for the swelling and deformation behaviors
of articular cartilage. Journal of Biomechanical Engineering, 113:245–258, 1991.
J. Lemaitre and R. Desmorat. Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle
Failures. Springer Science & Business Media, Berlin, Heidelberg, Germany, 2005.
V. C. Li. On engineered cementitious composites (ECC). Journal of Advanced Concrete Technology,
1:215–230, 2006.
A. I. Lurie. Nonlinear Theory of Elasticity. North Holland Series in Applied Mathematics and
Mechanics, Elsevier Science, Netherlands, 1990.
G. A. Maugin. The Thermomechanics of Nonlinear Irreversible Behaviours: An Introduction. World
Scientiﬁc Publishing Company, New Jersey, USA, 1998.
K. B. McAfee. Stress-enhanced diﬀusion in glass I. Glass under tension and compression. Journal
of Chemical Physics, 28:218–226, 1958a.
K. B. McAfee. Stress-enhanced diﬀusion in glass II. Glass under shear. Journal of Chemical Physics,
28:226–229, 1958b.
A. N. Morozovska, E. A. Eliseev, A. K. Tagantsev, S. L. Bravina, L. Q. Chen, and S. V. Kalinin.
Thermodynamics of electromechanically coupled mixed ionic–electronic conductors: Deformation
potential, Vegard strains, and ﬂexoelectric eﬀect. Physical Review B, 83:195313, 2011.
M. K. Mudunuru and K. B. Nakshatrala. A framework for coupled deformation-diﬀusion analysis with application to degradation/healing. International Journal for Numerical Methods in
Engineering, 89:1144–1170, 2012.

38

A. Muliana, K. R. Rajagopal, and S. C. Subramanian. Degradation of an elastic composite cylinder
due to the diﬀusion of a ﬂuid. Journal of Composite Materials, 43:1225–1249, 2009.
E. R. Myers, W. M. Lai, and V. C. Mow. A continuum theory and an experiment for the ion–induced
swelling behavior of articular cartilage. Journal of Biomechanical Engineering, 106:151–158, 1984.
D. J. Naus. Primer on durability of nuclear power plant reinforced concrete structures–A review of
pertinent factors. Technical report, Oak Ridge National Laboratory (ORNL), NUREG/CR–6927,
2007.
R. W. Ogden. Nonlinear Elastic Deformations. Dover publications, New York, 1997.
L. Onsager. Reciprocal relations in irreversible processes. i. Physical Review, 37:405, 1931a.
L. Onsager. Reciprocal relations in irreversible processes. ii. Physical Review, 38:2265, 1931b.
S. T. Peng and R. F. Landel. Induced anisotropy of thermal conductivity of polymer solids under
large strains. Journal of Applied Polymer Science, 19:49–68, 1975.
R. Picard and R. Leis. Some remarks on the horizontal line method. Mathematical Methods in the
Applied Sciences, 2:471–479, 1980.
K. R. Rajagopal, A. R. Srinivasa, and A. S. Wineman. On the shear and bending of a degrading
polymer beam. International Journal of Plasticity, 23:1618–1636, 2007.
E. Rothe. Zweidimensionale parabolische randwertaufgaben als grenzfall eindimensionaler randwertaufgaben. Mathematische Annalen, 102:650–670, 1930.
M. H. Sadd. Elasticity: Theory, Applications, and Numerics. Academic Press, Oxford, third edition,
2014.
J. Plešek and A. Kruisová. Formulation, validation and numerical procedures for Hencky’s elasticity
model. Computer and Structures, 84:1141–1150, 2006.
G. C. Sih, J. G. Michopoulos, and S. C. Chou. Hygrothermoelasticity. Martinus Nijhoﬀ Publishers,
Dordrecht, The Netherlands, 1986.
R. M. Springman and J. L. Bassani. Mechano-chemical coupling in the adhesion of thin-shell
structures. Journal of the Mechanics and Physics of Solids, 57:909–931, 2009.
S. H. Strogatz. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry,
and Engineering. Westview press, Cambridge, MA, 2001.
R. N. Swamy. The Alkali-Silica Reaction in Concrete. CRC Press, New York, USA, 2002.
J. F. Ulm, O. Coussy, L. Kefei, and C. Larive. Thermo–chemo–mechanics of ASR expansion in
concrete structures. Journal of Engineering Mechanics, 126:233–242, 2000.
R. Černy and P. Rovnaníková. Transport Processes in Concrete. CRC Press, New York, USA, 2002.
D. C. Venerus, J. D. Schieber, V. Balasubramanian, K. Bush, and S. Smoukov. Anisotropic thermal
conduction in a polymer liquid subjected to shear ﬂow. Physical Review Letters, 93:098301, 2004.
G. Z. Voyiadjis and P. I. Kattan. Damage Mechanics. Taylor & Francis Group, CRC Press, Boca
Raton, Florida, USA, 2005.
S. Wang and V. C. Li. High-early-strength engineered cementitious composites. ACI Materials
Journal, 103:97–105, 2006.
Y. J. Weitsman. Coupled damage and moisture-transport in ﬁber-reinforced, polymeric composites.
International Journal of Solids and Structures, 23:1003–1025, 1987.
Y. J. Weitsman. Anomalous ﬂuid sorption in polymeric composites and its relation to ﬂuid-induced
damage. Composites Part A: Applied Science and Manufacturing, 37:617–623, 2006.
Y. J. Weitsman and Y. J. Guo. A correlation between ﬂuid-induced damage and anomalous ﬂuid
sorption in polymeric composites. Composites Science and Technology, 62:889–908, 2002.
39

S. Wiggins. Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer Science &
Business Media, New York, second edition, 2003.
K. Willam, I. Rhee, and Y. Xi. Thermal degradation of heterogeneous concrete materials. Journal
of Materials in Civil Engineering, 17:276–285, 2005.
R. Zheng, R. I. Tanner, and X. J. Fan. Injection Molding: Integration of Theory and Modeling
Methods. Springer, Berlin, Heidelberg, 2011.
H. Ziegler. An Introduction to Thermomechanics. North Holland Publishing Company, Amsterdam,
Netherlands, 1983.

Ri
pi

Ro

po = atm
Figure 5. Inflation of a degrading spherical shell: A pictorial description of degrading shell
in the reference configuration. The shell is subjected to an inner pressure pi and an outer
pressure po .

0.4

No degradation
t = 0.0
t = 0.01
t = 0.05
t = 0.10
t = 0.50

0.3

Tθθ

0.2
0.1
0.0
-0.1
-0.2
0.5

0.6

0.7

0.8

0.9

1.0

R
Figure 6. Inflation of a degrading spherical shell: This figure shows the hoop stress T θθ as
a function of R at various instants of time due to an inner pressure of pi = 0.5. Analysis
is performed under strongly coupled chemo-thermo-mechano degradation. Note that the
stress is increasing with time under degradation.

40

0.4

pi = 0.2

0.3

Tθθ

pi = 0.3
0.2

pi = 0.4

0.1

pi = 0.5

0.0
-0.1
-0.2
0.5

0.6

0.7

0.8

0.9

1.0

R
Figure 7. Inflation of a degrading spherical shell: This figure shows the hoop stress T θθ as
a function of R at t = 0.1 for various inner pressures pi . Analysis is performed for strongly
coupled chemo-thermo-mechano degradation. T θθ increases in a non-linear fashion as the
pressure loading increases, which is different from the case as time progresses.

Chemical potential

0.0

pi = 0.2, ω = 0.05
pi = 0.3, ω = 0.05
pi = 0.4, ω = 0.05

-0.5

pi = 0.5, ω = 0.05
pi = 0.2, ω = 1.00
pi = 0.3, ω = 1.00

-1.0

pi = 0.4, ω = 1.00
pi = 0.5, ω = 1.00
No degradation

-1.5
0.5

0.6

0.7

0.8

0.9

1.0

R
Figure 8. Inflation of a degrading spherical shell: This figure shows the chemical potential
as a function of the reference location R at t = 0.2 due to various inner pressures pi under
different cases. One can see that for non-degrading shell, the chemical potential is unchanged
with respect to pressure loading. However, for strong coupling, it increases with pi in a
non-linear fashion when ω is small enough. This is because for small ω, diffusion takes
the dominance in the coupling effect. When pressure loading increases, the diffusivity is
increasing due to the growing tr[E]. For large ω, the deformation is dominant in the coupling,
which is −I E term in chemical potential. Since the first invariant I E is always positive in
this problem, chemical potential is decreasing when the pressure loading increases.

41

0.70

0.75

ri

0.65
0.60

0.65

ri

0.70

Strong coupling at t = 0.1
Strong coupling at t = 0.5
Strong coupling at t = 2.0
Weak coupling at t = 0.1
Weak coupling at t = 0.5
Weak coupling at t = 2.0
No degradation

0.60

Strong coupling at t = 0.1
Strong coupling at t = 0.5
Strong coupling at t = 2.0
Weak coupling at t = 0.1
Weak coupling at t = 0.5
Weak coupling at t = 2.0
No degradation

0.55

0.55
0.50
0.2

0.3

0.4

0.50
0.2

0.5

0.3

pi

0.4

0.5

pi

(a) Thermo-dominated degradation

(b) Chemo-dominated degradation

Figure 9. Inflation of a degrading spherical shell: This figure shows the plot of ri as a
function of the inner pressure pi for strongly and weakly coupled chemo-thermo-mechano
degradation problem. Note that in weak coupling the heat conductivity and diffusivity are
both constants, while the Lamé parameters still depend on concentration and temperature.
We take µ1 = 0.3 and µ2 = 0.4 for thermo-dominated degradation. For chemo-dominated
degradation, we have µ1 = 0.7 and µ2 = 0.1. For a given pi , one can see that r i for
weak coupling is larger than strong coupling when thermal degradation dominates. This is
because I E is always positive in this problem, the thermal conductivity decreases due to
the increase in I E . However, when moisture-induced degradation dominates, r i for weak
coupling is smaller than strong coupling problem. From this figure, we can observe creep-like
behavior for all the case studies.
1.0

Extent of damage (µ/µ0)

Extent of damage (µ/µ0)

1.0

0.8

0.6

0.4

0.2
0.5

No degradation
t = 0.00
t = 0.01
t = 0.05
t = 0.10

0.6

0.7

0.8

0.9

1.0

R

0.8

0.6

0.4

0.2
0.5

No degradation
t = 0.00
t = 0.01
t = 0.05
t = 0.10

0.6

0.7

0.8

0.9

R

(a) Thermo-dominated degradation

(b) Chemo-dominated degradation

Figure 10. Inflation of a degrading spherical shell: This figure shows the extent of damage as
a function of the reference location at various instants of time due to inner pressure pi = 0.5.
Different values are chosen for µ1 and µ2 for thermo-dominant and chemo-dominant degradation. Analysis is performed for strongly coupled case. For thermo-dominated problem,
healing-like behavior is observed at early time steps. This is because at initial times, we
have variable heat sinks in the entire body. As ϑ ≤ ϑ0 , the material damage is less than
that of at time t = 0 (but still remains below that of the virgin material). However, this
heal-like behavior becomes less distinct (or even doesn’t exist) when the chemo-degradation
achieves the dominance.
42

1.0

0.8

Extent of damage (µ/µ0)

Extent of damage (µ/µ0)

0.8
0.7
0.6
0.5

pi = 0.2

0.4

pi = 0.4

pi = 0.3

0.7
0.6

pi = 0.2

0.5

pi = 0.3
pi = 0.4

0.4

pi = 0.5
0.3
0.5

0.6

pi = 0.5
0.7

0.8

0.9

0.3
0.5

1.0

0.6

0.7

R

0.8

0.9

R

(a) Weak coupling degradation

(b) Strong coupling degradation

Figure 11. Inflation of a degrading spherical shell: This figure shows the extent of damage
as a function of the reference location at t = 1 for various inner pressures ‘pi ’. Analysis
is performed for thermo-dominated degradation. As the pressure increases, for the weakly
coupled problem, the extend of damage decreases. This means that when the inflation
pressure pi increases, the body degrades more significantly. However, this is not the case for
the strongly coupled problem. In this particular case, thermo-mechano coupling dominates
and plays a vital role. As I E ≥ 0, the strain-dependent thermal conductivity decreases as the
pressure loading increases. Hence, there is less damage in the material due to the decrease
in temperature values as compared to weakly coupled chemo-thermo-mechano degradation
problem.

PSfrag replacements

y

2L

x
r

Y
2W

θ
Oref

X

rc
x

Ocurr

Reference conﬁguration

Current conﬁguration

Figure 12. Bending of a degrading beam: A pictorial description of degrading beam in
both reference and current configurations. Bending moment is applied at the two ends of
the beam just after time t = 0. Oref and Ocurr correspond to the origin in the reference and
current configurations.

43

1.0

5

5

0

-5
-1

-0.5

0
X

No degradation
t = 0.0
t = 0.1
t = 0.2
t = 0.5
t = 1.0
1
0.5

No degradation
t = 0.0
t = 0.1
t = 0.2
t = 0.5
t = 1.0

Tθθ

10

Tθθ

10

0

-5
-1

(a) Weak coupling degradation

0
X

-0.5

0.5

(b) Strong coupling degradation

Figure 13. Bending of a degrading beam: This figure shows the plot of T θθ as a function of
the reference location of the cross-section at various instants of time. The stress distribution
is not linear, which is the case for finite deformation beam bending problem. Herein, we
observe that the neutral axis shifts further to the left. Moreover, in case of weak coupling
for some instants of time the maximum stress does not occur at either tensile or compressive
sides of the beam after the onset of degradation.

6.0

Three-way strong coupling
Three-way weak coupling

M

5.0

4.0

3.0
0.0

0.2

0.4

0.6

0.8

1.0

t
Figure 14. Bending of a degrading beam: This figure shows the plot of bending moment at
various instants of time for both strong and weak coupling chemo-thermo-mechano degradation. Moment relaxation is observed for both cases, however, in weak coupling the moment
declines at a much faster rate than that of the strong coupling case. Note that the bending
moment is a constant without degradation.

44

1

Chemical potential

0.0

-0.5

Strong coupling, t = 0.0
Strong coupling, t = 0.1
Strong coupling, t = 1.0
No degradation, t = 0.0
No degradation, t = 0.1
No degradation, t = 1.0

-1.0

-1.5

-2.0
-1

0
X

-0.5

1

0.5

Figure 15. Bending of a degrading beam: This figure shows the plot of chemical potential as
a function of the reference location of the cross-section at various instants of time when there
is no degradation and for strong coupling cases. In the strong coupling scenario, although
diffusion process is dominant, one can still observe that the deformation has a significant
effect on chemical potential as compared with non-degradation case.
0.8

Extent of damage (µ/µ0)

Extent of damage (µ/µ0)

1.0
0.8
0.6
No degradation
t = 0.0
t = 0.1
t = 0.2
t = 0.5
t = 1.0

0.4
0.2
0.0
-1

-0.5

0
X

0.5

1

(a) Strong coupling degradation

t = 0.0
t = 0.1
t = 0.2
t = 0.5
t = 1.0

0.6

0.4

0.2

0.0
-1

-0.5

0
X

0.5

(b) Weak coupling degradation

Figure 16. Bending of a degrading beam: This figure shows the extent of damage as a
function of the reference location of the cross-section at various instants of time (due to the
application of bending moment). Note that analysis is performed for both strongly coupled
and weakly coupled chemo-thermo-mechano degradation. One can see that a virgin beam
which is initially homogeneous after degradation is not homogeneous anymore. In addition,
the extent of damage is monotonic for weak coupling, which is not the case for strong
coupling. Such a phenomena has implications in damage control and retrofitting of the
degrading beams.

45

1

Ri

Ro

M

R

Z
θ
X

Y

Figure 17. Torsional shear of a degrading cylinder: A pictorial description of the degrading
cylinder under torsion in the reference configuration. Ri and Ro are, respectively, the inner
and outer radii of the cylinder. X, Y , and Z are the Cartesian coordinates in the reference
configuration. The bottom of the cylinder is fixed and a twisting moment is applied at the
top of the cylinder for t ≥ 0.

46

1.0

µ1 = 0.1

0.9

µ1 = 0.2

0.7

µ1 = 0.4

0.6

µi = 0.5

M

0.8

µ1 = 0.3

0.5
0.00

0.02

0.04

0.06

0.08

0.10

t
(a) Moment under different µ1

1.2

Ri = 0.2
Ri = 0.4

0.9

Ri = 0.6

M

Ri = 0.8

0.6

0.3
0.00

0.02

0.04

0.06

0.08

0.10

t
(b) Moment under different Ri

Figure 18. Torsional shear of a degrading cylinder: This figure shows the twisting moment
at various instants of time due to a given angle of twist per unit length of the cylinder,
Ψ1 = 0.75. One can see that as µ1 increases the twisting moment required to keep Ψ1
unchanged, decreases. Similar type of behavior is observed when µ1 is kept constant and
µ2 is varied. Herein, the main observation is that moment relaxation not only depends on
material degradation but also on the geometry of the degrading body.

47

0.20

No degradation
Strong coupling
Weak couling

N

0.15

0.10

0.05
0.00

0.02

0.04

0.06

0.08

0.10

t
Figure 19. Torsional shear of a degrading cylinder: This figure shows the non-dimensionalized
normal force N due to Poynting effect at various instants of time. Analysis is performed
for a given angle of twist per unit length of the cylinder, Ψ1 = 0.75. When there is no
degradation the normal force is constant. However, due to degradation one can see that
the normal force relaxes over time. The decrease in this normal force for weak-coupling is
higher than that of the strong coupling.

48

