Modeling hard-soft block copolymers as a liquid crystalline polymer
M. Manav
∗
Mechanical and Process Engineering, ETH Zurich, Switzerland
M. Ponga
†
Mechanical Engineering, University of British Columbia, Vancouver, Canada
M. Ortiz
‡
California Institute of Technology, Pasadena, USA
(Dated: May 16, 2023)
We report a new computational approach to model hard-soft block copolymers like polyurea as a
liquid crystalline polymer to understand their microstructural evolution due to mechanical loading.
The resulting microstructure closely resembles the microstructure observed in polyurea. The stress-
strain relations in uniaxial compression and tension loading obtained from the model are also in close
quantitative agreement with the experimental data for polyurea. We use the model to elucidate the
evolution of the hard and the soft domains during loading, which is consistent with the experimental
measurements characterizing microstructural evolution in polyurea.
I. INTRODUCTION
Thermoplastic elastomers (TPEs) such as polyurea [1]
exhibit excellent mechanical strength and toughness as
well as resistance to shock and impact loading [2–
5]. These properties have enabled the application of
polyurea in bulletproof vests and helmets [6, 7], blast
mitigation [8–10], among others, and have led to a
surge in scientific interest in modeling it and under-
standing the microstructural origins of its mechanical re-
sponse [5, 11, 12]. These remarkable mechanical prop-
erties of polyurea are often attributed to its hierarchical
microstructure. A polyurea molecule consists of alter-
nate blocks of flexible aliphatic segments and stiff seg-
ments made of aromatic moieties. In polyureas with long
flexible segments, the two blocks phase-segregate form-
ing ribbon-shaped hard domains (
≈
5
−
10 nm wide with
≈
7 nm separation in polyurea P1000) made of stiff seg-
ments dispersed in a soft matrix of flexible blocks [13–18].
The hard domains are glassy due to the hydrogen bond
network and
π
−
π
stacking of the neighboring aromatic
moieties. They act as physical crosslinks as well as stiff-
eners. However, the hierarchical arrangement also makes
it challenging to understand the microstructural origin of
its properties.
A number of experimental studies have been conducted
to characterize the evolution of polyurea microstructure
due to mechanical loading. Wide angle X-ray diffraction
(WAXD) studies have shown that the mean spacing be-
tween hard segments does not change during quasi-static
uniaxial tensile loading, though it leads to their align-
ment along tensile direction [14, 17]. The alignment is
strain-rate dependent and observed to be absent at high
∗
mmanav@ethz.ch
†
mponga@mech.ubc.ca
‡
ortiz@caltech.edu
strain-rates [13]. Small angle X-ray scattering (SAXS)
measurements reveal an axial increase and a transversal
decrease in the interdomain spacing of hard domains with
increasing strain [14]. They also indicate a breakdown in
the hard domain network structure at true strains above
∼
0
.
4. In situ AFM tapping-phase study of polyurea also
has revealed fragmentation of hard domains upon rapid
loading and various nano and mesoscale phase transitions
upon relaxation under fixed tensile strains [19]. A study
of the microstructure of polyurea under hydrostatic com-
pression in a diamond anvil cell using WAXD and SAXS
revealed the presence of interchain hydrogen-bonded net-
work even at pressures well above 1 GPa though crys-
tallinity is not observed [20]. The study also reported
phase-mixing with increasing pressure and pointed to
the contrasting behavior of polyurea under hydrostatic
compression and in tension where hard domains become
oriented. Investigation of the atomic-scale morphology
of the domains using multi-angle energy-dispersive X-
ray diffraction at pressures up to 6 GPa shows reduced
mean spacing between nearby hard segments in contrast
to minimal deformation within a hard segment and a de-
crease in interdomain order with increasing pressure [21].
Molecular dynamics (MD) simulation has also been
employed to investigate the mechanical behavior of
polyurea. An all-atom (AA) model of polyurea lever-
aging its phase-segregated microstructure and perform-
ing molecular dynamics simulation of each domain in-
dependently before combining their properties like in a
composite yielded close agreement with their experimen-
tally measured quasistatic response [22] and shock re-
sponse [23]. Molecular scale deformation mechanisms
in polyurea during mechanical loading are investigated
using AA MD simulation in [24]. Shock response of
polyurea has also been investigated using AA MD simu-
lation [25–27] and inelastic deformation, hydrogen bond
breaking and reforming, and permanent densification
have been identified as contributors to shock energy dis-
sipation. These models can typically access length and
arXiv:2305.07673v1 [cond-mat.soft] 11 May 2023
2
time scales only up to the order of 10 nm and 10 ns,
respectively, which is insufficient to capture the evolu-
tion of polyurea microstructure. To elucidate the effect
of the phase-segregated microstructure, coarse-grained
models [28, 29] and generic bead-spring polymer mod-
els with control over the strength of pair-interaction po-
tential to simulate hard and soft segments [30–32] have
also been reported. They suggest that the multiblock
architecture of polyurea leading to the formation of in-
terconnected hard domains results in improved mechan-
ical properties. A hybrid all-atom coarse-grained model
was employed to study the toughness of polyurea, and
self-reinforcement of oriented, hard segments and stress-
dependent release of soft segments from hard domains
are identified as two key mechanisms [12].
We report here a generic coarse-grained model of
polyurea with the aim of elucidating the evolution of their
microstructure in response to mechanical loading. In a
departure from the earlier models, we leverage the rigid-
ity and the shape anisotropy of the hard segments in
a polyurea molecule to model them as ellipsoidal par-
ticles like in a liquid crystalline polymer (Figure 1).
Since liquid crystals exhibit a complex phase diagram
in temperature-pressure-density space [33], this modeling
approach also facilitates an investigation of the influence
of their phase transitions on the mechanical response of
the modeled polymer. We show that the model produces
ribbon-like hard domains resembling that in polyurea and
the stress-strain relation obtained from the model is in
close agreement with the experimental data for polyurea
P1000 (henceforth referred to as polyurea) for which the
number of repeating units in soft segments is
n
= 14 (see
Figure 1(a)). We then study its microstructural evolu-
tion during uniaxial loading.
The paper is organized as follows. Section II details
the model and methodology. Section III presents and
discusses results from the model, including stress-strain
relation and microstructural evolution with mechanical
loading. The paper ends with conclusions and future
directions in Section IV.
II. MOLECULAR DYNAMICS SIMULATION
AND ANALYSIS
We have developed a generic coarse-grained model of
a molecule as shown in Figure 1 and performed molecu-
lar dynamics simulation using LAMMPS [34]. Lennard-
Jones (LJ) units have been used throughout this paper,
where
m, σ
and
are mass, length, and energy scales,
respectively. The time scale is denoted by the charac-
teristic time
τ
. We define a coarse-grained molecule by
connecting ellipsoidal particles and points particles (Fig-
ure 1). To create a bonded connection between a point
particle and
ends
of an ellipsoidal particle, we rigidly at-
tached
dummy
point particles of very small mass to the
ends of an ellipsoidal particle. Then, point particles were
connected to these dummy particles using bonded inter-
action. Each molecule consists of eight ellipsoids con-
nected with a chain of eight point particles in between
two consecutive ellipsoids. The number of point particles
is chosen to ensure that a phase-segregated microstruc-
ture forms when a melt of the polymer is cooled [32]. The
diameter of an ellipsoidal particle is 3
σ
,
σ
, and
σ
along its
three principal axes. The masses of an ellipsoidal parti-
cle, a point particle, and a dummy particle are 1
.
5708
m
,
0
.
5236
m
, and 0
.
0005
m
, respectively.
The total potential energy consists of pair energy
(
E
pair
) and bond energy (
E
bond
).
E
=
E
pair
+
E
bond
,
(1)
E
pair
=
i,j
∈
point particles
∑
j>i
U
LJ
ij
+
k
∨
l
∈
ellipsoids
∑
l>k
U
GB
kl
.
(2)
E
pair
accounts for the interactions between a pair of point
particles, a pair of ellipsoidal particles, and an ellipsoidal
and a point particle.
E
bond
is the energy associated with
bonded interaction.
In our model, nonbonded point particles interact with
a Lennard-Jones (LJ) potential, given as follows:
U
LJ
ij
= 4
(
(
σ
r
ij
)
12
−
(
σ
r
ij
)
6
)
−
U
shift
,
(3)
where
, and
σ
are characteristic energy and length. The
potential is cut-off at a distance
r
c
= 2
.
5
σ
and shifted by
U
shift
to ensure its continuity. In the simulations, we set
= 1 and
σ
= 1.
Interaction of an ellipsoidal particle with other par-
ticles is governed by Gay-Berne (GB) potential
U
GB
kl
,
which is a product of three terms, following [35–37]:
U
GB
kl
=
U
r
·
η
·
χ,
(4)
where
U
r
is a shifted LJ interaction, and distance-
independent interaction anisotropy is characterized by
η
and
χ
.
η
and
χ
control the interaction strength based
on particle shapes and relative potential well depths, re-
spectively. Shape for a particle
i
is specified by the shape
matrix
S
i
= diag(
a
i
, b
i
, c
i
), where
a
i
, b
i
and
c
i
are ellip-
soidal semiaxes. Orientation of particle
i
is given by
A
i
which represents the rotational transformation from the
simulation cell frame to the body frame. Relative energy
matrix is given by
E
i
= diag(
ai
,
bi
,
ci
), where
ai
,
bi
and
ci
are relative well depths. Then,
U
r
is given by:
U
r
= 4
k
(
ρ
12
−
ρ
6
)
,ρ
=
σ
h
+
γσ
,
(5)
where
and
σ
are the same as in LJ potential,
k
is used
to tune the depth of the potential well,
γ
is the shift
parameter and
h
is the distance of closest approach of
interacting particles and is approximated as:
h
=
r
−
(
1
2
ˆ
r
T
G
−
1
ˆ
r
)
−
1
/
2
,
(6)
3
FIG. 1. (a) A repeating unit of a polyurea molecule.
n
= 14 for polyurea P1000 used for comparison in this work. (b) A
repeating unit of a coarse-grained molecule, and (c) a coarse-grained molecule used in the simulation. Stiff blocks with aromatic
moieties in polyurea have been replaced with ellipsoidal particles in the coarse-grained molecule and flexible blocks with point
particles. To connect the point particles to the ends of the ellipsoidal particles, a dummy point particle with negligible mass is
first rigidly attached to the ends of ellipsoids. Adjacent point particles, point particles, and dummy particles are connected by
bonds.
where
r
=
r
l
−
r
k
is the vector between centers of ellip-
soids
k
and
l
with magnitude
r
=
√
r
·
r
and unit vector
ˆ
r
=
r
/r
. Also,
G
=
A
T
k
S
2
k
A
k
+
A
T
l
S
2
l
A
l
.
η
is given as follows:
η
=
(
2
s
k
s
l
det(
G
)
)
υ/
2
,
(7)
where
s
i
=
β
i
(
a
2
i
+
b
i
c
i
)(
b
i
c
i
)
1
/
2
for particle
i
, and
χ
is
given by the following expression:
χ
= (2
ˆ
r
T
B
−
1
ˆ
r
)
μ
,
(8)
where
B
=
A
T
k
E
−
1
/μ
k
A
k
+
A
T
l
E
−
1
/μ
l
A
l
.
υ
and
μ
are
empirical exponents used to tune the potential. The for-
mulation in [35], which is implemented in LAMMPS, as-
sumes that the longest ellipsoidal semiaxis is the third
semiaxis
c
i
as also pointed out in [36].
β
i
is introduced
here as a correction factor if that is not the case.
In our simulation, parameter values for ellipsoids are
chosen to be the classical values [38] for which phase-
diagram has been studied extensively [33].
We take
a
i
= 1
.
5
σ
and
b
i
=
c
i
= 0
.
5
σ
. Also,
ai
= 0
.
2 and
bi
=
ci
= 1, where the values are for particles interact-
ing along end-to-end, side-to-side, and face-to-face semi-
axes [39], respectively. We also take
γ
= 1,
υ
= 1, and
μ
= 2. The selected cut-off for ellipsoidal particles is 4
σ
.
Furthermore, since
a
i
is taken as the longest semiaxis for
ellipsoid
i
, the correction factor is
β
i
= 4
√
3
/
10 for the
prescribed values of semiaxes.
Bonded interaction is governed by the Finitely exten-
sible nonlinear elastic (FENE) potential as given below.
U
b
=
−
1
2
KR
2
0
log
(
1
−
(
r
R
0
)
2
)
+
4
(
(
σ
r
)
12
−
(
σ
r
)
6
)
+
,
(9)
In the FENE potential
r
represents the bond length, and
R
0
= 1
.
5
σ
is the maximum bond length. At initial ele-
vated temperatures during simulation to prepare a well
mixed melt,
K
= 20
T
∗
/σ
2
was employed to control
maximum stretch in bonds, where
T
∗
is reduced tem-
perature and is given by
T
∗
=
T/
(
k
B
)
. Once the melt
was cooled down to
T
∗
= 1
.
5, we set
K
= 30
/σ
2
. We
employed time-reversible integrators presented in [40] to
evolve the translational and rotational motion of the rigid
bodies consisting of an ellipsoid and dummy particles,
and velocity-Verlet scheme for the translational motion
of point particles. Ellipsoidal particles have only five de-
grees of freedom (DOF) since they are not free to spin
around their longest principal axis. So, we block this
DOF in our simulation and remove it from the calculation
of temperature. We used the Nos ́e-Hoover thermostat for
4
point particles. Chain thermostat and barostat based on
generalized Nos ́e-Hoover approach [41] were used to con-
trol the temperature of the rigid bodies and the pressure
in the simulation cell. Note that while barostating is
based on the rigid bodies, the positions of all the parti-
cles in the simulation cell are updated. A timescale of
0
.
5
τ
was used for the thermostat and 5
τ
was used for
the barostat. Throughout the simulation, a time step of
0
.
005
τ
was employed.
A. Representative volume element (RVE)
To obtain glassy hard domains, we developed RVEs for
k
= 2
,
3 and 5, ensuring stronger interaction between
ellipsoids compared to the interaction between point par-
ticles and following [31, 32]. In our simulation, we placed
1000 molecules in a simulation box in a random config-
uration at a high temperature to obtain a well mixed
melt. The initial melt was equilibrated at a high tem-
perature and a small pressure, before being cooled down
to
T
∗
= 1
.
5 and releasing the pressure. Simulation cells
become unstable at high temperatures in an NPT simu-
lation due to the limited cut-off range of the pair po-
tentials. The pressure was applied to keep the simu-
lation cells stable. The raised temperature and pres-
sure for
k
= 2 were
T
∗
= 2 and
p
∗
= 0
.
1 (
p
∗
is re-
duced pressure,
p
∗
=
p/
(
σ
3
)
). They were set to be
T
∗
= 2
.
5 and
p
∗
= 0
.
3 for
k
= 3, and
T
∗
= 3 and
p
∗
= 0
.
6 for
k
= 5. Then the cells were equilibrated
at
T
∗
= 1
.
5 for 10
6
time steps. All the equilibrated
melt at
T
∗
= 1
.
5 were cooled to a reduced temperature
of
T
∗
= 0
.
55. For
k
= 5, ellipsoidal particles phase
separate even at
T
∗
= 1
.
5. But for the lower values
of
k
, we observed phase-segregated microstructure with
ribbon-shaped hard domains at
T
∗
= 0
.
55. The simu-
lation cell at
T
∗
= 0
.
55 was further equilibrated for at
least 5
×
10
6
time steps before using it as an RVE for the
production run. Since the glass transition temperature
of a polymer melt of LJ particles with
= 1 and
σ
= 1
is
≈
0
.
45 [31, 42], it is in a melt state at
T
∗
= 0
.
55.
However, the ellipsoidal particles are in solid state, as
indicated by their translational diffusion and rotational
autocorrelation.
B. Mean squared displacement and autocorrelation
We calculate mean squared displacement (
MSD
) and
rotational autocorrelation (
C
R
) for the ellipsoidal parti-
cles, which are given as follows:
MSD
(∆
t
) =
〈|
r
i
(
t
0
+ ∆
t
)
−
r
i
(
t
0
)
|
2
〉
,
(10)
C
R
(∆
t
) =
〈
u
i
(
t
0
+ ∆
t
)
u
i
(
t
0
)
〉
,
(11)
where ∆
t
is the time interval,
r
i
(
t
) and
u
i
(
t
) are the
position and orientation vectors of the
i
th
ellipsoid at
time
t
, and
〈·〉
is phase average.
C. Microstructure characterization
To characterize the hard domain, its radial distribution
function (RDF) and orientational order parameter are
calculated. The nematic tensor is defined as:
Q
=
1
2
〈
3
u
i
u
T
i
−
I
〉
,
(12)
where
I
is identity matrix. The largest eigenvalue of
Q
is
the order parameter
S
and the corresponding eigenvector
is the orientation of the director alignment. In uniaxial
compression, directors align in the plane perpendicular to
the direction of compression. Then the eigenvalue with
a negative value has the largest magnitude and we take
that as the order parameter.
To characterize the deformation of the soft domain,
we investigate the mean and distribution of the end-to-
end distance
R
ee
of the soft segments connecting two
ellipsoidal particles in a molecule.
D. Uniaxial mechanical loading
In uniaxial tension and compression tests, the sample is
stretched (compressed) at a fixed engineering strain rate
along
x
-direction. In
y
and
z
-directions, pressure is main-
tained at zero. We have performed tests at three strain
rates: 5
×
10
−
3
τ
−
1
, 5
×
10
−
4
τ
−
1
and 2
.
5
×
10
−
5
τ
−
1
.
III. RESULT AND DISCUSSION
In this section, we present and discuss the results from
the simulation, including initial configuration, mechan-
ical response, and microstructure evolution due to the
applied compressive and tensile uniaxial loading.
A. Initial configuration
Figure 2(a) shows a well-mixed and homogeneous melt
at
T
∗
= 1
.
5 for
k
= 3. It forms phase-segregated do-
mains at
T
∗
= 0
.
55 as shown in Figure 2(b). Also,
the phase-segregated hard domains are ribbon-shaped
as exhibited by Figure 2(c) in close agreement with
the phase-segregated domains in polyurea. Furthermore,
even though the ribbon-shaped domains appear to be dis-
connected in Figure 2(b), cluster analysis confirms that
they are connected and form a scaffold as seen in Fig-
ure 2(c).
To elucidate the glassy nature of the hard domain,
mean squared displacement and rotational autocorrela-
tion of the ellipsoidal particles are shown in Figure 3.
The small magnitude of
MSD
and very slow decay of
C
R
over the simulation time, which is of the order of the
inverse of the slowest applied strain rate in uniaxial tests,
allude to their being in a glassy state. Also, a stronger
5
(b)
(a)
(c)
FIG. 2. (a) Homogeneous polymer melt at
T
∗
= 1
.
5, (b) phase-separated RVE at
T
∗
= 0
.
55, and (c) same as (b) but with
translucent soft domain exhibiting ribbon-shaped hard domain.
0 2000 4000 6000 8000 10000
12000 14000
0
2
4
6
8
10
12
0
2000 4000 6000 8000 10000
12000 14000
0.9
0.92
0.94
0.96
0.98
1
FIG. 3. Evolution of the mean squared displacement (top),
and the rotational autocorrelation (bottom) of the ellipsoidal
particles with time in equilibrium condition at
T
∗
= 0
.
55.
interaction between the ellipsoidal particles (higher
k
)
leads to a slower decay of
C
R
.
B. Mechanical response
We applied uniaxial strain in tension and compression
at three different strain rates to obtain the mechanical
response of the RVE. Figure 4 compares the resulting
stress-strain response for
k
= 2
,
3 and 5 with the ex-
perimental data for polyurea [2, 13, 43]. Stress has been
0
0.2
0.4
0.6
0.8
0
0.2
0.4
0.6
0.8
1
0
0.2 0.4 0.6 0.8
1
0
0.2
0.4
0.6
0.8
1
FIG. 4. Uniaxial stress-strain curves in compression (top) and
tension (bottom) at engineering strain rate ̇
e
= 5
×
10
−
5
τ
−
1
in simulation is compared with the experimental data for
polyurea. Experimental data in compression and tension are
taken from [43] and [2, 13], respectively.
normalized by dividing it with the stress at the maximum
simulated true strain (
ε
T
=
−
0
.
98 in compression and
ε
T
= 1
.
2 in tension). Since stress
∼
energy/length
3
, the
normalization helps remove the effect of energy scale (
)
and length scale (
σ
) and bares the effect of microstruc-
ture and time scale. Since strain rates is simulation is
expected to be higher than the experimental data, the
simulated response in compression conforms to the ef-
fect of strain rate. The response in tension however does
not conform to the expected strain rate effect. Never-