DOI:10.14311/APP.2021.30.0131
Acta Polytechnica CTU Proceedings 30:131–134, 2021 © Czech Technical University in Prague, 2021

available online at http://ojs.cvut.cz/ojs/index.php/app

MODELLING OF TIME-DEPENDENT BEHAVIOUR OF
PARTICULATE THERMOSET POLYMERS

Jan Vozáb∗, Jan Vorel

Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7,
166 29 Prague 6, Czech Republic

∗ corresponding author: jan.vozab@fsv.cvut.cz

Abstract. A preliminary study of a numerical model describing the behaviour of polymer-based
composites is presented. The numerical model consists of three main parts. The first is the microplane
M4 model, which is the main part of the model and is used to simulate elastoplastic behaviour and
damage. The second part consists of a generalized Maxwell model, which adds the effect of linear
creep of the material to the calculation. The last part is a free volume model that extends the linear
creep to the nonlinear creep. The creep is calculated on the deviatoric part of the normal stress of each
microplane, which allows the model to capture the polymer behaviour adequately without adjusting
the free volume of the model.

Keywords: Creep, free volume, microplane, thermosets, viscoelasticity.

1. Introduction
Many different materials are used in building con-
structions. One specific type that is increasingly
used is epoxies. These polymers are often used to
form different types of composites, as universal adhe-
sives, rigid foams, or as structural adhesives [1]. Like
many polymeric materials, they are subject to non-
linear creep, which is associated with susceptibility
to changes in temperature and humidity [2–4].

To properly design components made of this mate-
rial, the appropriate numerical tool is essential. Many
models, describing broad nonlinear time-dependent
responses of these material types, have been devel-
oped [5, 6]. By combining these approaches with
the generally accepted material model for mechani-
cal loading, the overall material response can be ad-
equately described.

In this paper, we introduce the novel approach
combining a microplane material model M4 [7, 8],
which is enhanced with the free volume approach on
the microplane level to capture the time-dependent
behaviour [5].

2. Numerical model
Simulating nonlinear behaviour is challenging be-
cause no nonlinear theory has been introduced that
can describe all materials. However, many models
have been developed that describe broad nonlinear
responses [5, 6]. By combining these approaches, it is
then possible to describe the behaviour of these ma-
terials. The numerical model described in this paper
consists of three main parts. The microplane material
model M4 [7, 8] with improvements presented in [9]
is used as a fundamental part. The Maxwell chain on
the microplane level is used to capture the creep be-
haviour [10]. The additional free volume approach [5]

allows for describing nonlinear viscoelastic behaviour
under small to moderate strains.

2.1. Microplane model M4
The microplane material model M4 [7] is the fourth
version of the microplane model introduced by
Bažant [11, 12]. This constitutive model is character-
ized by the relationships between stresses and strains
projected onto a surface of a microplane, which is ori-
ented by its own normal vector. The basic idea is that
the strain vector −→εN on the microplane is a projection
of a strain tensor ε. Then the normal strain on the
microplane have the form

εN = Nij εij . (1)

Similarly, the shear components, which are char-
acterized in the direction M and L, are given by or-
thogonal vectors, −→m and

−→
l , take the form

εM = Mij εij ; εL = Lij εij . (2)

in which Mij = (minj ) + (mj ni)/2 and Lij = (linj +
lj ni)/2 [13, 14]. All three components can be seen in
Fig.1(a).

The static equivalence is computed approximately
by the principle of virtual work written for the surface
Ω of a unit hemisphere

2π
3

σij δεij =
!

Ω
(σN δεN + σLδεL + σM δεM )dΩ. (3)

This equation represents that the virtual work of
macro-stresses within a unit sphere must be equal
to the virtual work of micro-stresses regarded as the
tractions on the surface of the sphere [7]. Substituting
δεN = Nij δεij , δεM = Mij δεij and δεL = Lij δεij ,
and noting that the last variational equation must

131

http://dx.doi.org/10.14311/APP.2021.30.0131
http://ojs.cvut.cz/ojs/index.php/app


Jan Vozáb, Jan Vorel Acta Polytechnica CTU Proceedings

x
3

x
2x

1

(a)

e
2

e
1

e
3

e
M

e
L

e
N

e
T

e

(b)

Figure 1. Microplane model: (a) distribution of
integration points (microplane normals) - system of
21 microplanes per hemisphere; (b) microplane strain
components

include any variation δεij , the result is the following
basic equilibrium equation [12]:

σij =
3

2π

!

Ω
sij dΩ ≈ 6

Nm"

µ=1
wµs

(µ)
ij , (4)

with
sij = σN Nij + σM Mij + σLLij . (5)

This integral is an approximation by an optimal
Gaussian integration formula for a spherical surface,
as you can see in Fig.1(b). This surface is represent-
ing a weighted sum over the microplanes. In the finite
element algorithms, this integral has to be computed
at each integration point of each finite element in each
time step [7].

The most generalized constitutive relations on the
microplane level have the forms

σN (t) = F tτ =0 [εN (τ ), εM (τ ), εL(τ )] , (6)
σL(t) = Gtτ =0 [εN (τ ), εM (τ ), εL(τ )] , (7)

σM (t) = Htτ =0 [εN (τ ), εM (τ ), εL(τ )] , (8)

where F , L and H are functionals of the microplane
strains in time t.

Volumetric/deviatoric decomposition of the normal
strain component [15] needs to be employed to cap-
ture the Poisson ratio of polymers adequately

εD = εN − εv ; εv = εkk/3, (9)

where εv stands for the volumetric strain and D de-
notes the deviatoric component.

2.2. Maxwell chain
The formulation of the microplane model M4 briefly
described in the previous section is extended to in-
corporate time-dependent behaviour of epoxies. Ap-
proach similar to [10] is utilized. Creep is developed
at every integration point, i.e. on each microplane of
every integration point. On the level of microplanes,
creep is implemented as linearly viscoelastic.

The behavior of each Maxwell element (µ), where
a spring and a dashpot are connected in series, is
expressed as

ε̇µ =
σ̇µ
Eµ

+
σµ
ηµ

, (10)

where Eµ is the elastic modulus and ηµ is the viscosity
of the element. In the proposed model, the viscosity
is assumed only for the deviatoric part of the nor-
mal strain component on each microplane separately.
Thus, according to the approximation by the gener-
alized Maxwell model, the equation for calculation
of deviatoric stress predictor in the incremental form
reads

σveD = σ
i
D + E

′′
D (εD − ε

i
D ) − ∆σ

′′
D =

= σiD + E
′′
D (∆εD − ∆ε

′′
D ), (11)

where i means the initial value in the time step ∆t
and ve stands for the viscoelastic stress increment,
which is the generalization of the elastic stress incre-
ment. Variables labeled by ′′ are calculated according
to the equations

E′′D =
N"

µ=1

1 − e−∆t/τµ
∆t/τµ

ED,µ, (12)

∆σ′′D =
N"

µ=1
(1 − e−∆t/τµ )σiD,µ, (13)

∆ε′′D = ∆σ
′′
D /E

′′
D , (14)

where N is the number of Maxwell elements and τµ =
ηµ/ED,µ relaxation time of the µ-th element. If a
constant Poisson’s ratio ν is assumed

ED,µ =
5Eµ

(2 + 3γ)(1 + ν)
, (15)

where γ is a parameter that may be chosen or can be
optimized so as to match the given test data best (this
value is set to 1 in this study) [16]. Note that remain-
ing constitutive relations remain unchanged and the
moduli are updated based on E′′ = E′′D (2+3γ)(1+ν).

2.3. Free volume
The last extension of our model is a free volume ap-
proach to extend the viscoelastic behaviour to a non-
linear scale, which is more convenient for the mate-
rial. The description of this approach can be found

132



vol. 30/2021 Modelling of particulate thermoset polymers

Figure 2. Reological model of the Maxwell chain.

τµ Eµ [MPa]

1 10 2.6174
2 102 107.5458
3 103 455.0239
4 104 933.0051
5 105 562.4520
6 106 273.8973
7 107 198.4917
8 108 135.9774
9 109 41.7351
10 1010 4.0175
11 1011 4.7209
12 1012 -0.0819
13 ∞ 39.9691

Table 1. Used values of Maxwell chain. Values con-
verted from [5, Tables 1 and 2].

in [5]. In the employed approach, the dilation is in-
troduced into the time scale, which results in a non-
linear system of equations. Note that for the case
of infinitesimal dilatation, the nonlinear free volume
theory reverts to linear viscoelasticity. In general,
temperature (T ), moisture content (c) and mechani-
cal dilatation (θ) influence the time scale of the ma-
terial in a similar manner such that

a = a(T, c, θ). (16)

The shift factor a is connected to the free volume as
suggested in [17, 18]

log10 a =
B

2.303

#
1
f

−
1
f0

$
, (17)

where f is fractional free volume, f0 represent the
fractional free volume at a reference state. The cou-
pling of the free volume to material parameters that
affect its macroscopic volume can be expressed using
a linear dependence

f = f0 + αv ∆T + βv ∆c + δθ, (18)

where αv is the volumetric thermal expansion of the
free volume, βv represent the volumetric expansion
due to a change in moisture content and δ relates the
change if the free volume due to mechanical volume
changes of the polymer. Combining Eq. (17) with
Eq. (18) we get

log a =
−B

2.303f0

#
αv ∆T + βv ∆c + δθ

f0 + αv ∆T + βv ∆c + δθ

$
. (19)

This equation includes the influence of all essential
properties related to the development of the volume
changes. The influence of the shift factor on the stress
response is through convolution, or superposition, in-
tegrals relating the time dependent stress response to
the strain history [5] by

σD =
! t

0
E′′D (t

′ − τ ′)
∂εD
∂τ

dτ, (20)

where the kernel takes the form

t′ − τ ′ =
! t

τ

dξ

a[T (ξ), c(ξ), θξ]
(21)

to account for the temperature, moisture and dilata-
tion histories. Note that in [5], the modified free vol-
ume approach is utilized to simulate shear-dominated
loading scenarios with little or no dilatation. Such
an approach is not needed in the current formulation
since the creep behaviour is defined on microplanes.

Parameter Value
B 0.5
f0 0.1
δ 1

αv 2.857 × 10−3 °C−1
T0 55 °C

Table 2. Used material parameters. Values taken
from [5].

3. Results
In this phase of the development process, we used a
single-point calculation to verify the numerical model
using [19]. The values for the Maxwell chain from Ta-
ble 1 and the material parameters from Table 2 were
used in the simulation. The element was loaded by
pure shear deformation. Fig. 3 shows the force vs.
displacement diagram. The blue line represents the
elastic loading of the element, and it represents the
upper bound for the presented results. The red line
shows the elastic loading of the material, including
the effect of the Maxwell chain. As expected, creep
affects the magnitude of the reaction, and therefore
the curve has a smaller slope. The last green line rep-
resents the investigated microplane model. Initially,
it has the same inclination as the elastic model with
the Maxwell chain, but as the deformation increases,
the response deviates. Note that for this study, the
microplane model behaviour is assumed to be purely
elastic.

133



Jan Vozáb, Jan Vorel Acta Polytechnica CTU Proceedings

Figure 3. Results from finite element analysis: Cube
with an edge 1 cm, loading pure shear.

4. Conclusions
This paper presented a new model that can be used to
describe the nonlinear creep of thermosets, which is
composed of three parts: (i) the microplane model
M4, which calculates elastoplasticity and damage;
(ii) the Maxwell chain model, which adds the effect
of linear creep to the calculation; (iii) the free vol-
ume model, which extends creep to a nonlinear scale.
In Fig.3, we can see that according to preliminary
results, the model behaves as expected because the
effect of creep and softening is noticeable under shear
load. At the same time, it is necessary to test the
model on a larger scale at different load cases to en-
sure its accuracy. The Arcan experiment presented
in [5] would be appropriate. Furthermore, it is nec-
essary to test the model under a combination of dif-
ferent temperatures and humidity for verification of
the free volume model parameters. If the accuracy
of the model is still not sufficient, it could be caused
by curing of the material which could be additionally
connected to our numerical model.

Acknowledgements
The financial support provided by the GAČR grant
No. 19-15666S and by the Czech Technical University in
Prague within SGS project with the application registered
under the No. SGS20/155/OHK1/3T/11 is gratefully ac-
knowledged.

References
[1] M. A. Meyers, K. K. Chawla. Mechanical Behavior of

Materials. Cambridge University Press, 2nd edn., 2008.
[2] J.-P. Pascault, H. Sautereau, J. Verdu, R. J. Williams.

Thermosetting polymers, vol. 64. CRC Press, 2002.
[3] R. F. Landel, L. E. Nielsen. Mechanical properties of

polymers and composites. CRC press, 1993.
[4] P. Silva, T. Valente, M. Azenha, et al. Viscoelastic

response of an epoxy adhesive for construction since its

early ages: Experiments and modelling. Composites
Part B: Engineering 116:266–277, 2017.

[5] C. F. Popelar, K. M. Liechti. Multiaxial Nonlinear
Viscoelastic Characterization and Modeling of a
Structural Adhesive. Journal of Engineering Materials
and Technology 119(3):205–210, 1997.
https://asmedigitalcollection.asme.org/
materialstechnology/article-pdf/119/3/205/
5681124/205_1.pdf doi:10.1115/1.2812245.

[6] G. Cusatis, D. Pelessone, A. Mencarelli. Lattice
discrete particle model (ldpm) for failure behavior of
concrete. i: Theory. Cement and Concrete Composites
33(9):881–890, 2011.

[7] Z. P. Bažant, F. C. Caner, I. Carol, et al. Microplane
model m4 for concrete. i: Formulation with
work-conjugate deviatoric stress. Journal of
Engineering Mechanics 126(9):944–953, 2000.

[8] F. C. Caner, Z. P. Bažant. Microplane model m4 for
concrete. ii: Algorithm and calibration. Journal of
Engineering Mechanics 126(9):954–961, 2000.

[9] G. Di Luzio. A symmetric over-nonlocal microplane
model m4 for fracture in concrete. International journal
of solids and structures 44(13):4418–4441, 2007.

[10] Z. P. Bažant, F. C. Caner, M. D. Adley, S. A. Akers.
Fracturing rate effect and creep in microplane model
for dynamics. Journal of engineering mechanics
126(9):962–970, 2000.

[11] Z. Bažant. Microplane model for strain-controlled
inelastic behaviour, chapter 3. Mechanics of
engineering materials, chapter 3 1984.

[12] Z. P. Bažant, P. G. Gambarova. Crack shear in
concrete: Crack band microplane model. Journal of
Structural Engineering 110(9):2015–2035, 1984.

[13] Z. P. Bažant, P. C. Prat. Microplane model for
brittle-plastic material: I. theory. Journal of
Engineering Mechanics 114(10):1672–1688, 1988.

[14] Z. P. Bažant, P. C. Prat. Microplane model for
brittle-plastic material: Ii. verification. Journal of
Engineering Mechanics 114(10):1689–1702, 1988.

[15] I. Carol, Z. Bažant. Damage and plasticity in
microplane theory. International Journal of Solids and
Structures 34(29):3807–3835, 1997.

[16] Z. Bažant, Y. Xiang, P. Prat. Microplane model for
concrete. I: Stress-strain boundaries and finite strain.
Journal of Engineering Mechanics 122(3):245–254,
1996.

[17] A. K. Doolittle. Studies in newtonian flow. ii. the
dependence of the viscosity of liquids on freeâĂŘspace.
Journal of Applied Physics 22(12):1471–1475, 1951.
https://doi.org/10.1063/1.1699894
doi:10.1063/1.1699894.

[18] A. K. Doolittle, D. B. Doolittle. Studies in newtonian
flow. v. further verification of the freeâĂŘspace viscosity
equation. Journal of Applied Physics 28(8):901–905,
1957. https://doi.org/10.1063/1.1722884
doi:10.1063/1.1722884.

[19] ES3. Mars - modeling and analysis of the response
structures.

134

https://asmedigitalcollection.asme.org/materialstechnology/article-pdf/119/3/205/5681124/205_1.pdf
https://doi.org/10.1115/1.2812245
https://doi.org/10.1063/1.1699894
https://doi.org/10.1063/1.1699894
https://doi.org/10.1063/1.1722884
https://doi.org/10.1063/1.1722884