Smart Materials and Structures

A Novel Statistical Spring-Bead Based Network Model for
Self-Sensing Smart Polymer Materials
Jinjun Zhang1, Bonsung Koo1, Yingtao Liu2, Jin Zou1, Aditi Chattopadhyay1 and Lenore Dai1
1

School for Engineering of Matter, Transport, and Energy, Arizona State University, Tempe, AZ, USA 85287
2

School of Aerospace and Mechanical Engineering, University of Oklahoma, Norman, OK, USA 73019

This paper presents a multiscale modeling approach to simulate the self-sensing
behavior of a load sensitive smart polymer material. A statistical spring-bead based
network model is developed to bridge the molecular dynamics simulations at the nanoscale
and the finite element model at the macroscale. Parametric studies are conducted on the
developed network model to investigate the effects of the thermoset crosslinking degree on
mechanical response of the self-sensing material. The comparison between experimental
and simulation results shows that the multiscale framework is able to capture the global
mechanical response with adequate accuracy and the network model is also capable of
simulating the self-sensing phenomenon of the smart polymer. Finally, the molecular
dynamics simulation and network model based simulation are implemented to evaluate
damage initiation in the self-sensing material under monotonic loading.

Nomenclature
I

N

̅

=
=
=
=
=
=
=
=
=
=

fluorescence intensity density
original length of the ith spring
variation of the spring length
number of springs in this subzone
number of actual potential crosslinked covalent bonds
number of all potential crosslinked covalent bonds
loading strain
equivalent strain in the ith spring
average equivalent strain in the bead-centered subzone
crosslinking degree

Keywords: spring-bead model, network model, statistical distribution of crosslinking degree, self-sensing smart
material

I. Introduction

L

oad sensitive materials that exhibit specific chemical reactions upon external loading (also known as
mechanophores) have received considerable attention in recent years. [1-4] Mechanophores, which are covalently
linked into polymers, is a result of the transmission of external forces at the macroscale to the debonding of covalent
links at the molecular level. [4,5] Therefore, itprovides a novel means for designing smart materials with self-sensing
capability. Cho et al. [6,7] have developed a cyclobutane-based mechanophore comprised of carbon-carbon covalent
bonds that transform into cinnamate groups and emit visible fluorescence under external loading. Zou et al.[8] have
developed a novel self-sensing technique in which Tris-(Cinnamoyl oxymethyl)-Ethane (TCE) is embedded as the
cyclobutane-based mechanophore material within a polymer matrix. In this study, cyclobutane-containing
crosslinked polymers are incorporated into an epoxy matrix and the effect on their mechanical properties are
examined, followed by a demonstration of early damage detection through mechanically induced fluorescence. This
experimental research has shown significant promise in developing mechanophore embedded polymer materials that
1

Smart Materials and Structures
can be used to detect damage initiation in polymeric composite materials. [8] However, the discovery of new
mechanophores (i.e., mechanochemically reactive units) is currently time-consuming, expensive, low-yielding, and a
process that is not easily scaled to large quantities. High fidelity modeling techniques are necessary to understand
the physics of mechanophores. A will developed modeling approach can simulate the multiscale response of these
“smart” polymers, which will be useful to design new polymers with mechanophores as well as apply the available
materials for large scale engineering applications.
In epoxy based polymers, the crosslinked bond clusters between the resin and the hardener molecules generate a
large-scale network structure. The heterogeneous crosslinked network at the microscale has significant impact on the
local mechanical properties.[9-11] In order to design and synthesize polymers with multifunctional capabilities, it is
important to develop a fundamental understanding of the effects of microscale variability on material properties and
response, such as the Young’s modulus, heterogeneous stress/strain distribution and damage initiation mechanisms.
Currently, there are numerous approaches to the modeling of polymeric systems integrated with nanoparticles; they
include homogenization techniques, the finite element (FE) method, molecular dynamics (MD) simulation, Monte
Carlo (MC) simulation, and the Mori-Tanaka approach.[12-14] MD simulations and FE analysis are two of the most
widely used methods for understanding the material behavior of epoxy-based systems at the nano and microscales.
MD simulation is capable of effectively capturing the heterogeneous material properties of a polymer that are
influenced by its crosslinking degree.[9-11] Similarly, Fan et al.[11] performed MD simulations for characterizing
material properties of a crosslinked epoxy resin compound. Using this model, they were able to capture linear
thermal expansion coefficients and Young's modulus of the material. MD simulations have also been used to
investigate the strength and molecular mechanisms of adhesion between an inorganic substrate and a cured epoxy
resin.[15] In this model, the crosslink density and corresponding material properties of the crosslinked system were
accurately estimated. However, the high computational cost associated with large-scale MD simulations is a major
restriction and often limits its application beyond the nanoscale.
FE modeling has been a common method to model material performance, primarily because it is
computationally efficient.[16-18] Onck et al. developed a two-dimensional model to study the strain stiffening of
crosslinked filamentous polymer.[19] The FE method was used to discretize the filament with Euler-Bernoulli beam
elements that accounted for stretching and bending. Jiang et al. performed a parametric study using FE analysis to
compute the material properties and scratch behavior of a crosslinked polymer. [20] The simulation results in this
work corresponded well with the mechanical properties obtained from the experiments. Chen and Lagoudas
developed a constitutive theory for polymers in which the crosslinking density was considered as the major
influence on the mechanical behavior. [21] This model has been implemented and extended by other researchers to
study thermoelasticity, viscoelasticity, and nonequilibrium relaxation. [22-24] While the FE method provides
advantages in high simulation efficiency to model polymers, it cannot be used to model heterogeneous
microstructure because the definition of discontinuous material properties is often infeasible in FE approach.
Therefore, new techniques with improved computational efficiency and high accuracy are desired for investigating
the material behavior of polymers at the microscale.
In addition to FE and MD models, another method commonly adopted is the discrete material techniques, such
as the spring-bead model, have been developed to study polymer material. [25-27] Underhill and Doyle developed a
new method for generating coarse-grained models of polymers. [28] To represent the larger scale material behavior of
polymers, spring-bead chains have been further crosslinked to generate a network model. [29] A spring-bead based
network model that bridges nanoscale information to higher length scales, has been developed to model the
structural epoxy embedded with TCE monomer. [8,30] The material was partitioned into discrete mass beads, which
were then linked using linear springs at the microscale. By integrating multiple spring-bead models, a network
model was constructed to represent the material properties at the mesoscale. A series of MD simulations were
performed initially to define the spring stiffness in the statistical network model using a well-known MD simulation
tool called LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator), and a force field for general
organic materials called MMFF (Merck molecular force field). [31,32] Compared with the MD model, only the
equivalent mechanical responses that significantly influence the damage initiation was captured in the spring-bead
based network model. Additional material properties, such as thermal and chemical properties, were omitted in order
to improve simulation efficiency. The developed multiscale methodology was demonstrated as being
computationally efficient and able to provide a possible means to bridge various length scales (from 10 nm in MD
simulation to 10 mm in FE model) without significant loss of accuracy. In this paper, the statistical network model
is further investigated to understand the effects of the model parameters on material response and damage in the
TCE embedded polymer. The model is also used to develop a relationship between mechanical strain and
experimentally observed fluorescence, with damage growth in the smart polymer under multiaxial loading.

2

Smart Materials and Structures

In a crosslinked network model, material performance is affected by parameters, such as crosslinking degree,
chain flexibility and chain lengths. [33-35] Among all variabilities, crosslinking degree, which represents the degree of
microscopic crosslinked covalent bonds, is considered as one of the major effects on the mechanical properties of
polymers. [33, 36-38] Urbaczewski-Espuche et al. developed a crosslinked network model and studied the influence of
the crosslinking degree on the mechanical properties of epoxy. [39] The results showed that the mechanical properties
including Young’s modulus, ultrasonic modulus, thermomechanical relaxation, and plastic behaviors are all affected
by crosslinking degree. Svaneborg et al. investigated the effect of microscopic disorder on the mechanical response
of network models with the same crosslinking degree. [40] They concluded that the heterogeneities in the randomly
generated crosslinked network models cause significant differences in the localization of monomers. An important
objective of this work was to understand the microscopic performance of smart polymeric material through a
parametric study of the network model. A series of simulations were performed to compare the effects of the
crosslinking degree on the mechanical response, such as strain distribution. The knowledge obtained was important
for enhancing material performance, and for methods required for controlling crosslinking degree during the
manufacturing process. In addition, this study led to a preliminary investigation for evaluating damage accumulation
and growth under monotonic loading.
In this paper a multiscale modeling approach is developed to simulate the mechanophore reactions of an epoxy
based smart polymer system bridging between the nanoscale and micro/macroscales. In Section II, the developed
spring-bead network model is introduced in details. The parametric study of the spring-bead model is discussed in
Section III. The multiscale network model is presented in Section IV. The comparison of the experimental validation
of the simulated damage is shown in Section V.

II. Development of spring-bead based network model
Zhang and Chattopadhyay recently developed a novel multiscale model that can characterize the material
behavior of the crosslinked polymer. [30] The schematic of the multiscale approach adopted in this work is shown in
figure 1. A brief description of the model is presented here as a comprehensive overview. At the microscale, a
spring-bead model was used to capture the mechanical response of a crosslinked cluster of polymer chains. MD
simulations were performed to provide local material properties represented by springs. At the mesoscale, the
network model consisting of a series of springs and beads was developed to represent larger scale crosslinked
polymer chains and for capturing the mechanical responses. Mechanical equivalence optimization was then
implemented to identify the relationship between the crosslinking degree and its corresponding representative spring
stiffness in the network model. This network model was capable of bridging the high accuracy MD model at the
microscale and the computationally efficient FE model at the macroscale, serving as a simplified mechanical
equivalent of MD models with significant reduction in computational effort.

Figure 1. Modeling of smart polymer materials at the microscale, mesoscale, and macroscale.
The multiscale approach developed in this paper was as follows: First, the material behavior of the
heterogeneous crosslinked microstructure was determined using MD simulation. Di-Glycidyl Ether of Bisphenol F
(DGEBF) and Di-Ethylene Tri-Amine (DETA), respectively, were used for the molecular structures of the epoxy
resin and the hardener that form the polymer. In order to capture the variability in the crosslinking degree and the
mechanical properties of the polymer system, a statistical approach was used, whereby a representative volume
3

Smart Materials and Structures

element (RVE) model was constructed with molecules of the resin, hardener, and smart material, which were
randomly distributed in the initial configuration (see figure 2). The weight percentages of the components in this
smart polymer were 70.9% for DEGBF, 19.1% for DETA, and 10.0% for TCE. In the crosslinked polymer material,
the crosslinking degree, which is defined as the percentage of actual crosslinked covalent bonds over all potential
crosslinked bonds, significantly influences the mechanical properties at the microscale. It is important to accurately
capture the crosslinking degree using MD simulations. The crosslinking degree can be defined in Eq. (1):
(1)
where is the crosslinking degree and
and
are the number of the actual and all potential crosslinked
covalent bonds in the same zone at the microscale, respectively. In this approach, a total of 500 MD simulations
were performed to represent the crosslinking degree of the target system, each with the same molecules (42,350
atoms in each RVE) but randomly generated initial configuration (different initial positions). The size of RVE
models was 14nm × 14nm × 14nm. After performing energy minimization based on the conjugate gradient
algorithm to avoid undesired repulsive forces between atoms, iso-baric/iso-thermal (NPT) ensemble simulations was
conducted to equilibrate each RVE at standard room condition (300K and 1atm). Periodic boundary conditions
(PBC) were applied to the boundaries of the RVE. The equilibrated RVE was cured using a bond formation
command in LAMMPS with a cut-off radius of 4.0 Angstroms; this cut-off radius was selected because it is more
than twice the length of a carbon-nitrogen covalent bond. Based on the simulation results there was a quasi-normal
distribution for the crosslinking degrees, of which the range was between 40% and 65% and the average value was
52.72% (see figure 3). MD simulations were further performed on the crosslinked polymeric systems to obtain the
mechanical responses (tensile, compressive, and shear moduli) associated with variable crosslinking degrees. The
moduli can be described as functions of the crosslinking degree. These results were used to construct the springbead based network model.

Figure 2. RVE model of smart polymer; comprises epoxy resin, hardener, and smart material (grey, cyan, blue, red
sphere represent hydrogen, carbon, oxygen, and nitrogen atoms, respectively).
The next step was to develop a spring-bead based network model to characterize the material behavior of a
crosslinked cluster of polymer chains with inputs from MD simulations. In the spring-bead model (see figure 1)
each of the two beads served as half of the bond cluster mass; the linear spring was used to represent the elastic
mechanical response (tensile and compressive) of the bond cluster. Variable spring stiffness (ratio of force to
relative deformation) was used to represent bond clusters with different crosslinking degrees. To represent larger
scale material behavior, a network model consisting of a series of spring-bead models was developed (see figure 1).
A set of network models with different combination methods was generated. The “8 neighbor based” network model
was then found to be optimal having structural stability, good representation of isotropic property, and high
computational efficiency. [30] Therefore, in this paper, all parametric studies and analyses were performed using the 8
neighbor based network model.
Mechanical equivalence optimization was implemented to integrate the mechanical responses from the MD
simulations and the network model. A multiobjective optimization technique was used to establish a relation
between the local crosslinking degree and the associated representative spring stiffness in the network model.

4

Smart Materials and Structures

Further details of the multiscale modeling approach of the smart polymer material can be found in the authors’
previous work. [30]

Figure 3. Statistical distribution of crossing-linking degrees based on MD simulation.

Figure 4. Polynomial fitting curve of spring stiffness based on optimization results.

III. Parametric study of spring-bead based network model
A. Homogeneous and heterogeneous distribution of crosslinking degree
In this section, a parametric study of the network model is conducted to understand the effect of the most
critical microscopic variable, the crosslinking degree, which influences the distribution of stress, strain, strain
energy, etc. The network model was constructed by discrete beads and crosslinked springs instead of continuous
material. Therefore, the traditional strain concept was not applicable. A new definition of strain is required prior to
the parametric study and analyses of the network model. In traditional continuous materials strain is defined as a
geometrical measure of deformation representing the relative displacement between particles in a material body. In
the spring-bead based network model the “equivalent strain” in one spring is defined as the deformation in terms of
original and relative length of this spring, as shown in Eq. (2):

5

Smart Materials and Structures

‖

(2)

‖

where is the original length of the ith spring,
is the variation of the spring length, and is the equivalent
strain in the spring. In the 8 neighbor based network model, each bead is linked with eight neighbor beads by springs
(see figure 1 and figure 5). The black solid lines and yellow circles in figure 5 represent the springs and beads,
respectively. To evaluate material deformation in a zone, the average equivalent strain must be defined in terms of
relative springs within the same zone. The average equivalent strain in the bead-centered area (the square zone with
grey color in figure 5) is defined as the averaged value of all strain in the eight neighbor springs, which are then lin
to the centered bead. For the 8 neighbor based network model, the average equivalent strain in the grey area is
defined in Eq. (3):
(3)

̅

∑

where is the equivalent strain in the ith spring, ̅ is the average equivalent strain in the bead-centered subzone, and
N is the number of springs in this subzone. N = 8 for the 8 neighbor based network model.

Figure 5. Beads and springs in the 8 neighbor based network model.
The influence of heterogeneous distribution of the crosslinking degree on the mechanical response is studied. A
set of four network models (size: 1mm × 1mm) were generated. There are 100 (10 × 10) subzones in each of the
four network models. As shown in figure 6(a), the schematic network model consists of two subzones (A and B)
representing two different crosslinking degrees. The degree of heterogeneity, which is defined as the difference in
crosslinking degree between neighbor subzones, for all four network models is present in Table 1. It can be seen
that although the average crosslinking degree value remains the same value (53%), the degree of heterogeneity
increases gradually from model 1 to model 4. Note that the crosslinking degree value is directly related to the
material’s properties, which are represented by spring stiffness. The spring stiffness in subzones A and B, for all
four models, is shown in figure 6(b). The increasing difference in spring stiffness between neighbor subzones also
indicates that the material properties change from homogeneous to heterogeneous.
Table 1. Crosslinking degree and degree of heterogeneity in the four network models
crosslinking degree in A

crosslinking degree in B

heterogeneous degree

Model 1

53%

53%

0

Model 2

56%

50%

6%

Model 3

59%

47%

12%

Model 4

62%

44%

18%

6

Smart Materials and Structures

a)
b)
Figure 6. a) Specific network model to study effect of degree of heterogeneity and b) spring stiffness in subzone A
and B of the network model 1, 2, 3 and 4.
Simulations were conducted using uniaxial tensile load applied to all four network models. The loading is strain
controlled and the maximum loading strain is 5%. The equivalent strain is used to evaluate the mechanical responses
in the network models. The simulation results are shown in Table 2 and figure 7. In all four models, the average
equivalent strain is nearly the same. It can be observed that the difference between the maximum and the minimum
equivalent strain in the network models is also influenced by the degree of heterogeneity (see Table 2). The
difference increases with degree of heterogeneity (see figure 7).

Figure 7. The difference in neighbor crosslinking degree influences the difference between maximum and minimum
equivalent strain in the network model.
Thus, the parametric study involving heterogeneous crosslinking degree distribution revealed that a larger
degree of heterogeneity leads to a higher difference of equivalent strain. It must be noted that in brittle materials,
strain gradient plays an important role in singular or non-singular stress concentrations and local damage. [41] It has
been found that a higher local strain gradient can accelerate damage initiation. [42,43] From the simulations and
analyses of results, the increasing degree of heterogeneity can increase the strain gradient, which leads to a higher
probability of damage evolution. These observations indicate that a reduction in the degree of heterogeneity can help
to reduce the damage initiation.

7

Smart Materials and Structures

Table 2. Equivalent strain in network models to study degree of heterogeneity effect
Model 1

Model 2

Model 3

Model 4

Average equivalent strain

0.0339

0.0339

0.0339

0.0338

Maximum equivalent strain

0.0358

0.0357

0.0358

0.0361

Minimum equivalent strain

0.0326

0.0320

0.0319

0.0319

Difference between maximum and
minimum equivalent strain

8.93%

10.36%

10.89%

11.63%

B. Homogeneous and statistical distribution of crosslinking degree
At the microscale of the self-sensing smart polymer, the crosslinking degrees are statistically distributed. This is
also the limitation of the traditional FE model, which is not capable of representing the crosslinked microstructure,
thereby serving as motivation for the network model presented here. To represent the statistical distribution of
crosslinking degrees at the microscale, a statistical network model is developed, as shown in figure 8(a). To validate
the microstructural representation capability of the statistical network model, a homogeneous network model is
constructed and compared with the statistical network model. All crosslinking degrees in the homogeneous network
model are kept the same, as shown in figure 8(b). For both the statistical and the homogeneous network models, the
same size and geometry as the model used in section III.A (10 × 10 subzones in 1mm × 1mm) is maintained. As
shown in figure 8(a) and figure 8(b), the numbers in the variable subzones represent the corresponding local
crosslinking degrees (unit: %). The weight percentage of the variable crosslinking degrees in the statistical network
model is governed by the probability density function (PDF) obtained from the MD simulations (see figure 3). More
[
details about the PDF of the crosslinking degrees can be found in the authors’ previous work. 30] For the
homogeneous network model, uniform crosslinking degree, which is equal to the averaged value of all crosslinking
degrees in the statistical network model, (52.72%) is used for all subzones.

a)
b)
Figure 8. Crosslinking degree distribution in a) statistical network model and b) homogeneous network model
(unit: %).
Uniaxial tensile test is performed in both statistical and homogeneous network models. The loading is strain
control and the maximum loading strain is 5%. The simulation results are shown in figure 9(a) and figure 9(b),
which provide the equivalent strain distribution in the statistical and homogeneous network models, respectively.
Results of the average, maximum and minimum equivalent strain in the statistical and homogeneous network
models, are presented in Table 1.
As shown in Table 3, the average values of equivalent strain in the statistical network model and the
homogeneous network model are close. This implies that the overall mechanical response of the statistical network
model corresponds well with the homogeneous model under monotonic loading. The difference between the
8

Smart Materials and Structures

maximum and the minimum strain for statistical network model (10.36%) is greater than the corresponding value for
homogeneous network model (8.40%). Based on the simulation results, a similar conclusion is obtained, as in
section III.A. The statistical distribution of crosslinking degrees results in a larger distinguishable disorder in strain
distribution indicating a high strain gradient. Consequently, compared with the homogeneous microstructure the
statistical microstructure can lead to a higher damage initiation.

a)
b)
Figure 9. Equivalent strain distribution in a) statistical network model and b) homogeneous network model.
Table 3. Equivalent strain in statistical and homogeneous network models
Statistical

Homogeneous

Average equivalent strain

0.0338

0.0339

Maximum equivalent strain

0.0357

0.0357

Minimum equivalent strain

0.0320

0.0327

Difference between maximum and
minimum equivalent strain

10.36%

8.40%

C. Effect of spring length in network model
As discussed earlier, the two steps in developing a network model are: i) designing the configuration variables
including arrangement, spring size, and spring stiffness; ii) bridging MD simulation and network model using
mechanical equivalence optimization. In the first step, the definition of spring stiffness is related to the spring length.
In the second step, the mechanical equivalence provides the network model with the most proximate mechanical
response of the MD simulation. Technically, the spring length has a negligible effect on the mechanical response of
the network model due to the mechanical equivalence optimization conducted in the second step of the development
of the network model. To understand the spring length effect on the material behavior representation, two network
models with different spring lengths were generated. The sizes of the two network models were 8 μm × 8 μm and 1
mm × 1 mm respectively, and the corresponding spring lengths in the two network models were 0.8 μm and 0.1 mm
respectively. Based on the simulation results shown in Table 4, the material characteristics including tensile modulus,
compressive modulus and shear modulus represented by the two network models are close (error is less than 1%),
validating the hypothesis that there is negligible effect of spring length on the mechanical response of the network
model. However, shorter spring lengths can capture mechanical responses at shorter length scales, and thereby
enhance simulation accuracy. Additionally, corresponding computational efficiency is reduced because of the
increased number of springs. Spring length, therefore, plays a similar role as that of mesh size in an FE model.
Table 4. Mechanical properties in network models of different spring lengths.
Network model of short
springs (0.8 μm)
Network model of long

Tensile Modulus

Compressive Modulus

Shear Modulus

1.778 GPa

1.752 GPa

0.502 GPa

1.769 GPa

1.744 GPa

0.503 GPa

9

Smart Materials and Structures

springs (0.1 mm)
Error of mechanical
properties

0.51%

0.46%

0.20%

IV. Multiscale integration based on a statistical network model
A. Reliability of statistical network model
A set of three meso statistical network models were constructed to validate the simulation accuracy (see figure
10(b), figure 10(c) and figure 10(d)). For each the three network models, the same size and geometry as the model
used in section III.A (10 × 10 subzones in 1mm × 1mm) is maintained. There are 121 beads, 220 short springs, and
200 long springs in the 8 neighbors based network model. The lengths of the short and long springs in the network
models are 0.1 mm and 0.1414 mm, respectively. The numbers shown in the square subzones in figure 10(b), figure
10(c), and figure 10(d) represent the local crosslinking degrees.

a)

b)

c)

d)

Figure 10. a) Probability distribution of crosslinking degree in the statistical network model, b) network model 1,
c) network model 2, and d) network model 3.
Simulations were conducted using uniaxial strain controlled tensile load; the maximum loading strain was 10%.
The average equivalent strain in the three network models was calculated. It was observed that the average
equivalent strain in the three statistical network models was very close and the error was less than 1% (see figure
11). The results establish the validity and robustness of the representative models. It should be noted that the spring10

Smart Materials and Structures

bead based network model improves the computational efficiency significantly. For a similar uniaxial tensile loading
simulation 20 hours of runtime is required using MD modeling whereas a network model completes the simulation
in 4 minutes. The simulation efficiency is improved by 300 times using the network model.

Figure 11. Average equivalent strain in three statistical network models vs. loading strain.
B. Local strain and crosslinking degree in the network model
Simulations were conducted with the statistical network model to study the relation between local strain and
crosslinking degree using biaxial strain controlled tensile load; the maximum loading strain was 10% in both
directions. Results are presented for the statistical model 1, shown in figure 12(a). The distribution of equivalent
strain is shown in figure 12(b). There is an obvious relation between local strain and crosslinking degree: high strain
concentration (the red band in figure 12(b)) occurs mostly where there is low crosslinking degree (the red square
zone in figure 12(a)) especially where the crosslinking degree is below 52.72% (average crosslinking degree in the
model). These results are physically meaningful; lower material stiffness leads to a higher strain under the same
external loads. The locations with high strain concentration also correlate to potential local defects. The damage
evaluation based on MD & network model simulations will be discussed in section V.

a)
b)
Figure 12. Simulation of statistical network model under biaxial tensile load: a) crosslinking degree distribution in
network model and b) distribution of equivalent strain in network model.
C. Comparison between network model and FE model
The spring-bead based network model was developed in order to represent the smart polymer at the mesoscale.
To capture the material behavior from mesoscale to macroscale, a coupled network and an FE model are required.
The meso network model is explicitly modeled and embedded at the hot spot within a FE mesh of the structure so

11

Smart Materials and Structures

that different length scales can be bridged. A notched dog bone specimen was designed for simulation under
uniaxial tensile loads. The dimensions of the notched dog bone specimen are shown in figure 13. The network
model and FE model were used concurrently to represent the material properties at the notch tip and the remaining
section, respectively (see figure 14(a)). An FE model without embedded springs and beads was developed (see
figure 14(b)). To validate the network model based multiscale modeling approach two measures were considered:
hot spot location and far field distributions. To ensure the network model is representative of the bulk material and
mechanical responses both measures must remain consistent with or without a network model implemented in the
FE model.

Figure 13. Dimensional drawing of dog bone specimen for simulation.

a)
b)
Figure 14. Strain distribution in the hot spot using a) network model and b) FE model.
Comparison of the strain distributions obtained using the two models are shown in figure 14. In both, the
network model (figure 14(a)) and FE model (figure 14(b)), the high strain concentrations are observed near the
notch tip. However, the strain distribution in the network model shows larger variability due to the variability in the
crosslinking degree. The local strain distribution obtained using the network model is expected to provide a more
accurate estimate of the local defects at the microscale. figure 15 shows the comparison of strain distributions
outside the hot spots in the simulations, with and without a network model. A good correlation can be observed
between the network model and the FE model, validating the accuracy of the network model with the FE model
outside of the hot spot (notch tip).

12

Smart Materials and Structures

a)
b)
Figure 15. Strain distribution within notched specimen a) with network model and b) without network
model.

V. Damage estimation using statistical network model
A. Evaluation of self-sensing intensity of a smart polymer
In this research TCEs were used as the stress-sensing material to detect the defects in the smart polymer. The
epoxy resin is brittle with a propensity to crack if exposed to sufficient deformation, leading to irreversible damage
in the form of micro cracking. In the high strain concentration locations of the material the force can be transferred
through the bulk material to individual polymer chains, and then to the cleavable bonds on the mechanophore units.
The local deformation in these regions exceeds the critical limitation of the polymer chains and results in the
[
permanent deformation under rupture of crosslinked bonds. 44] This early damage precursor can be detected through
the employment of TCE-containing polymers. Note that the fluorescence response of the mechanophore results from
the rupture of the covalent bonds in cyclobutane. It is also related to the onset of permanent strain (deformation) of
the bonds at the microscale. The fluorescence intensity is also directly dependent on the density of the breaking
bonds. The novel spring-bead based network model representing this TCE embedded polymer is also required to
simulate this self-sensing feature. Therefore, there is a correlation between the local strain and the color changing
phenomena. This relation between fluorescence intensity in the experiments and the mechanical response in the
simulation using the network model was studied.
Uniaxial compressive tests were performed to obtain the fluorescence intensities of the smart material under
different loads. A Test Resources Electromechanical System (with the maximum loading capability 4.4 kN) was
used to conduct the compressive tests; stress-strain curves were also obtained. The tests were performed with
displacement control in a longitudinal direction at a loading rate of 1 mm/min. The maximum loading strain was
10% to ensure that the material deformation remained in the elastic range. To monitor the fluorescence from defects
in the polymer, the cubic samples of polymer blends were compressed to different strains and observed under the
fluorescence microscope (Nikon Elipse TE300 inverted video microscope) by exposure to 330-380 nm UV light. An
image processing program, ImageJ, was used to process the fluorescent micrographs and capture the fluorescence
intensities. The micrographs were taken at the same location of the specimen surface, with the size of 427 μm × 325
μm. The integrated fluorescence intensity was calculated through the summation of all pixel values in the
micrograph. A quasi-linear relation was observed between the fluorescence intensity density and the loading strain
(see figure 16(a)). Figure 16(b) and figure16(c) show the microscopic images of the specimen surface and the
corresponding fluorescence responses when the loading strain is 4% and 6%, respectively. The relation between
fluorescence intensity density and the loading strain can be described using a linear fit:
(4)

13

Smart Materials and Structures

where I represents the observed fluorescence intensity density, and

is the loading strain.

b)

a)
c)
Figure 16. Fluorescence observation of self-sensing smart polymer. a) Integrated intensity density vs. loading
strain in fluorescence observation, b) Microscopic images and detected fluorescence of specimen at 4% strain,
and c) Microscopic images and detected fluorescence of specimen at 6% strain.
To capture the self-sensing performance, simulation was performed using the statistical network model 1, as
shown in figure 10(b). Uniaxial compressive loads of the same magnitude as that of experiments were applied to the
network model. The equivalent strain in the network model was calculated, which shows a quasi-linear relation
between the average equivalent strain and the loading strain (figure 17(a)). The distribution of equivalent strain in
the network model under different loading strain can be seen in figure 17(b)-(e). The simulation result shows that the
strain distribution is influenced by the local crosslinking degree. Also, both local and average equivalent strain
increases linearly with loading. A linear fit is used to describe the relation between average equivalent strain and the
global loading strain:
(5)

̅

where is the loading strain, and ̅ is the average equivalent strain in the statistical network model. As shown in
Eqs. 4 and 5, both the fluorescence intensity and the strain in the network model exhibit a quasi-linear relationship
with respect to the loading strain (see figure 16(a) and figure 17(a)). A function describing the relation between
average strain ̅ in the network model and fluorescence intensity density I from experimental observation can be
derived by combining Eqs. (4) and (5):
169.09 ̅ +123.17

(6)

Thus, the average equivalent strain from the network model based simulation can be used to describe the observed
fluorescence intensity of the TCE-embedded polymer.
These results indicate that the average equivalent strain obtained from the network model can be used as an
index to describe the self-sensing fluorescence intensity. This approach of evaluating fluorescence intensity
numerically using the high efficiency network model can reduce the experimental cost significantly.

14

Smart Materials and Structures

b)

c)

a)
d)
e)
Figure 17. Network model based simulation results: a) Average equivalent strain in network model; and equivalent
strain distribution at different global strain at b) 3%, c) 5%, d) 8%, and e) 10%.
B. Strain based damage evaluation to predict failure
In the MD simulation uniaxial tensile tests were implemented on a series of RVE models with different
crosslinking degrees. The RVE size, molecular constituents and MD simulation conditions have been discussed in
section II. The simulation results showed that the stress (load) in the RVE model increased in a fluctuating curve
with strain (see figure 18). One important conclusion from the MD results is that the first significant kink of the
stress-strain curves takes place where the strain value is around 8% for a range of crosslinking degrees. This
phenomenon indicates that there is a notable microstructural defect resulting from the rupture of polymer chains at
the kink point. The relationship between the mechanical response at the kink point and the damage initiation has also
been observed in experiments and simulations conducted by other researchers. Kingsbury et al. obtained the
threshold kink point of the stress-strain curve of mechanophore crosslinked polymer in their experiments. 4 Their
results also corresponded with the onset of the mechanochemical activation in the fluorescence observation.
Silberstein et al. found a non-equilibrium state resulting in local force fluctuations in glassy polymers.5 The MD
simulation results showed good agreement between the onset of strain softening (first kink) and the onset of
mechanophore activation.

Figure 18. Tensile mechanical responses of smart polymer at variable crosslinking degrees based on the MD
simulation.

15

Smart Materials and Structures

Based on the damage analysis in the MD simulation, there are micro cracks when the local strain is 8%,
associated with the observed initial kink point of the stress-strain curve. Next, simulations were conducted using the
network model to observe the strain evolution. Statistical network model 1 (see figure 10(b)) was used under a
uniaxial tensile load. The equivalent strain distribution from the simulation results is shown in figure 19. It is
observed that the local equivalent strain increases with loading. The loading strain implemented on the network
model is 11.2%, when the maximum local equivalent strain exceeds 8%, indicating permanent tensile strain for this
smart material.

a)

b)

c)

d)

Figure 19. Equivalent strain distribution in statistical network model 1 at different loading strain: a) 5%, b) 8%, c)
10%, and d) 11%.
Experiments were conducted to validate the damage evaluation phenomenon of the smart polymer. A set of
three uniaxial tensile tests were conducted. The Test Resources Electromechanical System (with the maximum
loading capability 4.4 kN) was also used to conduct the tensile tests. From the experimental results the material
showed elastic behavior until failure (see figure 20). The average permanent tensile strain was 10.6% based on
experimental observation. The error of the permanent tensile strain between experiment and simulation was 5.36%.
This statistical network model, therefore, was capable of evaluating the damage initiation in the polymer material.

Figure 20. Experimental results of tensile tests.

VI. Conclusion
A multiscale approach is developed here to model the mechanical properties of a novel mechanophore
crosslinked polymer that is capable of detecting damage precursor. A spring-bead based network model is developed
to represent the self-sensing smart polymeric material. Three different length scales are considered. At the
microscale, the spring-bead model is constructed based on the MD simulation results to represent crosslinked
polymer chains. At the mesoscale, a network model is developed to capture the statistical material properties. At the
macroscale, a coupled network and the FE model is constructed to represent the mechanical response in the
structure. A series of parametric studies are performed to understand the microstructural behavior of the smart
material. The results show that the local crosslinking degree has a significant effect on the mechanical properties of
the material. A high crosslinking degree value leads to low local strain, which demonstrates stronger material
16

Smart Materials and Structures

performance. A larger difference between neighbor crosslinking degrees can increase the strain gradient, which then
leads to higher local defects. Results from the simulations and experiments were combined to establish a correlation
between the mechanophore fluorescence and local strain. A linear relation between average strain in the network
model and fluorescence intensity density in the experimental observation is also obtained. The computationally
efficient statistical network model is seen as capable of describing the self-sensing phenomenon using the
relationship between local equivalent strain in the network model and the observed fluorescence intensity. The
statistical network model is also used to estimate the damage initiation by analyzing the MD and network simulation
results. Good agreement was observed between simulation and experiment. Although the thermal and chemical
properties of the mechanophore are not addressed in the network model to improve simulation efficiency, this study
shows that the statistical network model can successfully capture the dynamic behavior of the self-sensing smart
polymer material.

Acknowledgements
This work is partially supported by the Air Force Office of Scientific Research (Grant number: FA9550-12-1-0331),
and the program manager is Dr. David Stargel. This work is also partially supported by the Office of Naval Research
(Grant number: N00014-14-1-0068), and the program manager is Dr. William Nickerson.

References
1

Thostenson, E. T., and Chou T. W. “Carbon nanotube networks: sensing of distributed strain and damage for
life prediction and self healing,” Advanced Materials, Vol.18, No.21, 2006, pp. 2837-2841.
2
Wu, D. Y., Meure S., and Solomon D. “Self-healing polymeric materials: a review of recent developments,”
Progress in Polymer Science, Vol. 33, No.5, 2008, pp. 479-522.
3
Black, A. L., Lenhardt, J. M., and Craig, S. L., “From molecular mechanochemistry to stress-responsive
materials.” Journal of Materials Chemistry, Vol. 21, No. 6, 2011, pp. 1655-1663.
4
Kingsbury, C. M., et al. “Shear activation of mechanophore-crosslinked polymers,” Journal of Materials
Chemistry, Vol. 21, No.23, 2011, pp. 8381-8388.
5
Silberstein, M. N., et al. “Modeling mechanophore activation within a crosslinked glassy matrix,” Journal of
Applied Physics, Vol. 114, No.2, 2013, pp. 023504.
6
Cho, S. Y., Kim J. G., and Chung C. M., “A fluorescent crack sensor based on cyclobutane-containing
crosslinked polymers of tricinnamates,” Sensors and Actuators B: Chemical, Vol. 134, No. 2, 2008, pp. 822-825.
7
Cho, S. Y., Kim J. G., and Chung C. M., “Photochemical Crack Healing in Cinnamate based Polymers,”
Journal of Nanoscience and Nanotechnology, Vol. 10, No. 10, 2010, pp. 6972-6976.
8
Zou, J., et al. “Early damage detection in epoxy matrix using cyclobutane-based polymers,” Smart Materials
and Structures, Vol. 23, No.9, 2014, pp. 095038.
9
Flory P. J., and Rehner Jr. J., “Statistical Mechanics of Cross‐Linked Polymer Networks II. Swelling,” The
Journal of Chemical Physics, Vol. 11, No.11, 1943, pp. 521-526.
10
Krumova M., et al. “Effect of crosslinking on the mechanical and thermal properties of poly (vinyl alcohol),”
Polymer, Vol. 41, No.26, 2000, pp. 9265-9272.
11
Fan H. B., and Yuen M. M., “Material properties of the crosslinked epoxy resin compound predicted by
molecular dynamics simulation,” Polymer, Vol. 48, No.7, 2007, pp. 2174-2178.
12
Ghosh, S., Lee, K., Moorthy, S., “Multiple scale analysis of heterogeneous elastic structures using
homogenization theory and voronoi cell finite element method,” Int. Journal of Solids and Structures, Vol. 32, No.
1, 1995, pp. 27–62.
13
Borodin O., et al, “Multiscale modeling of viscoelastic properties of polymer nanocomposites,” Journal of
Polymer Science Part B: Polymer Physics, Vol. 43, No. 8, 2005, pp. 1005-1013.
14
Fermeglia M., and Pricl S., “Multiscale modeling for polymer systems of industrial interest,” Progress in
organic coatings, Vol. 58, No.2, 2007, pp. 187-199.
15
Yarovsky I., and Evans E., “Computer simulation of structure and properties of crosslinked polymers:
application to epoxy resins,” Polymer, Vol. 43, No.3, 2002, pp. 963-969.
16
Neerukatti, R. K., Liu K.C., Kovvali N. and Chattopadhyay A., “Fatigue life prediction using hybrid prognosis
for structural health monitoring,” Journal of Aerospace Information Systems, Vol. 11, No. 4, 2014, pp. 211-232.
17
Peng, T., et al., “A novel Bayesian imaging method for probabilistic delamination detection of composite
materials,” Smart Materials and Structures, Vol. 22, No.12, 2013, pp. 125019.

17

Smart Materials and Structures
18

Peng, T., et al., “Integrated experimental and numerical investigation for fatigue damage diagnosis in
composite plates,” Structural Health Monitoring, 2014, 1475921714532992.
19
Onck, P. R., et al., “Alternative explanation of stiffening in crosslinked semiflexible networks,” Physical
review letters, Vol. 95, No.17, 2005, pp. 178102.
20
Jiang, H., et al., “Finite element method parametric study on scratch behavior of polymers,” Journal of
Polymer Science Part B: Polymer Physics, Vol. 45, No.12, 2007, pp. 1435-1447.
21
Chen, Y., and Dimitris C. L., “A constitutive theory for shape memory polymers. Part I: Large
deformations,” Journal of the Mechanics and Physics of Solids, Vol. 56, No.5, 2008, pp. 1752-1765.
22
Westbrook, K. K., et al., “A 3D finite deformation constitutive model for amorphous shape memory polymers:
a multi-branch modeling approach for nonequilibrium relaxation processes,” Mechanics of Materials, Vol. 43, No.
12, 2011, pp. 853-869.
23
Volk, B. L., Lagoudas, D. C., and Chen, Y. “Analysis of the finite deformation response of shape memory
polymers: II. 1D calibration and numerical implementation of a finite deformation, thermoelastic model,” Smart
Materials and Structures, Vol. 19, No. 7, 2010, pp. 075006.
24
Diani, J., et al. "Predicting thermal shape memory of crosslinked polymer networks from linear
viscoelasticity." International Journal of Solids and Structures 49.5 (2012): 793-799.
25
Buxton, G. A., and Balazs, A. C., “Lattice spring model of filled polymers and nanocomposites,” Journal of
chemical physics, Vol. 117, No. 16, 2002, pp. 7649-7658.
26
Iwata K., “Irreversible Statistical Mechanics of Polymer Chains. I. Fokker–Planck Diffusion Equation,”
Journal of Chemical Physics, Vol. 54, No. 1, 2003, pp. 12-26.
27
Schwarz U. S., Erdmann T., and Bischofs I. B., “Focal adhesions as mechanosensors: the two-spring model,”
Biosystems, Vol. 83, No. 2, 2006, pp. 225-232.
28
Underhill P. T., and Doyle P. S., “Development of bead-spring polymer models using the constant extension
ensemble,” Journal of Rheology, Vol. 49, No.5, 2005, pp. 963-987.
29
Indei, T., Schieber, J. D. and Takimoto, J., “Effects of fluctuations of crosslinking points on viscoelastic
properties of associating polymer networks,” Rheologica acta, Vol. 51, No.11-12, 2012, pp. 1021-1039.
30
Zhang, J., Koo, B., Subramanian, N., Liu, Y. and Chattopadhyay, A., “Statistical Multiscale Modeling of Smart
Polymer Materials using a Spring-Bead Based Network Model,” 55TH AIAA/ASME/ASCE/AHS/SC STRUCTURES,
STRUCTURAL DYNAMICS, AND MATERIALS CONFERENCE, 2014, January13-17, National Harbor, Maryland
31
Halgren, T. A., “Merck molecular force field. I. Basis, form, scope, parameterization, and performance of
MMFF94,” Journal of computational chemistry, Vol. 17, No.5‐6, 1996, pp. 490-519.
32
Plimpton, S., Crozier, P., and Thompson, A., “LAMMPS-large-scale atomic/molecular massively parallel
simulator,” Vol. 18, 2007, Sandia National Laboratories.
33
Espuche, E., et al., “Influence of crosslink density and chain flexibility on mechanical properties of model
epoxy networks,” Macromolecular symposia, Vol. 93., No. 1., 1995, Hüthig & Wepf Verlag.
34
Charulatha, V., and Rajaram, A., “Influence of different crosslinking treatments on the physical properties of
collagen membranes,” Biomaterials, Vol. 24, No.5, 2003, pp. 759-767.
35
Park, H. B., et al., “Effect of crosslinked chain length in sulfonated polyimide membranes on water sorption,
proton conduction, and methanol permeation properties,” Journal of membrane science, Vol. 285, No.1, 2006, pp.
432-443.
36
Zosel, A., and Ley, G., “Influence of crosslinking on structure, mechanical properties, and strength of latex
films,” Macromolecules, Vol. 26, No. 9, 1993, pp. 2222-2227.
37
Halary, J. L., “Structure–property relationships in epoxy-amine networks of well-controlled architecture,” High
performance polymers, Vol. 12, No.1, 2000, pp. 141-153.
38
Berger, J., et al., “Structure and interactions in covalently and ionically crosslinked chitosan hydrogels for
biomedical applications,” European Journal of Pharmaceutics and Biopharmaceutics, Vol. 57, No. 1, 2004, pp. 1934.
39
Urbaczewski‐Espuche, E., et al., “Influence of chain flexibility and crosslink density on mechanical properties
of epoxy/amine networks,” Polymer engineering & science, Vol. 31, No. 22, 1991, pp. 1572-1580.
40
Svaneborg, C., Grest, G. S., and Everaers, R., “Disorder effects on the strain response of model polymer
networks,” Polymer, Vol. 46, No. 12, 2005, pp. 4283-4295.
41
Li, J., et al., “A micromechanics-based strain gradient damage model for fracture prediction of brittle
materials–Part II: Damage modeling and numerical simulations,” International Journal of Solids and Structures,
Vol. 48, No.24, 2011, pp. 3346-3358.
42
Geers, M. G. D., et al., “Strain-based transient-gradient damage model for failure analyses,” Computer
Methods in Applied Mechanics and Engineering, Vol. 160, No.1, 1998, pp. 133-153.
18

Smart Materials and Structures
43

Francfort, Gilles A., and J-J. Marigo. "Revisiting brittle fracture as an energy minimization problem." Journal
of the Mechanics and Physics of Solids 46.8 (1998): 1319-1342.
44
Buehler, Markus J., Sinan Keten, and Theodor Ackbarow. "Theoretical and computational hierarchical
nanomechanics of protein materials: Deformation and fracture." Progress in Materials Science 53.8 (2008): 11011241.

19