Document [original]

Continuum Mech. Thermodyn.

https://doi.org/10.1007/s00161-020-00939-4

ORIGINAL ARTICLE

Bilen Emek Abali ·Jan Vorel ·Roman Wan-Wendner

Thermo-mechano-chemical modeling and computation

of thermosetting polymers used in post-installed fastening

systems in concrete structures

Received: 9 June 2020 / Accepted: 13 October 2020

Abstract As thermoset polymers find frequent implementation in engineering design, their application in

structural engineering is rather limited. One key reason relies on the ongoing curing process in typical applica-

tions such as post-installed adhesive anchors, joints by structural elements or surface-mounted laminates glued

by adhesive polymers. Mechanochemistry including curing and aging under thermal as well as mechanical

loading causes a multiphysics problem to be discussed. For restricting the variety of material models based

on empirical observations, we aim at a thermodynamically sound strategy for modeling thermosets. By pro-

viding a careful analysis and clearly identifying the assumptions and simplifications, we present the general

framework for modeling and computational implementation of thermo-mechano-chemical processes by using

open-source codes.

Keywords Continuum mechanics ·Thermosetting polymers ·Thermomechanics ·Curing ·Finite element

method (FEM)

1 Introduction

Legal documents require a long service life, especially in structural engineering, no less than 50 years is usual.

As a result, experimental procedures are costly and often not feasible. Nonetheless, a few testing techniques

and models exist aiming at a performance prediction for that time horizon given the coupled phenomena of

post-curing, creep, and fatigue at periodic thermal or mechanical loading. Further complications arise from the

high internal moisture of, e.g., concrete or other building materials which are in contact with the thermosets

and that give rise to degradation phenomena, such as hydrolytic aging [1–6] playing a significant role in the

Communicated by Luca Placidi.

B. E. Abali

Chair of Continuum Mechanics and Constitutive Theory, Institute of Mechanics, Technische Universität Berlin, Einsteinufer 5,

10587 Berlin, Germany

B. E. Abali (B

)·R. Wan-Wendner

Christian Doppler Laboratory LiCRoFast, University of Natural Resources and Life Sciences, Peter-Jordan-Straße 82,

1190 Vienna, Austria

E-mail: [email protected]

J. Vorel

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

R. Wan-Wendner

Magnel–Vandepitte Laboratory, Department of Structural Engineering and Building Materials, Ghent University,

Technologiepark-Zwijnaarde 60, 9052 Ghent, Belgium

B. E. Abali et al.

design-life of the structure. Different phenomena in many length scales [7–11] necessitate incorporation of

multiscale approaches. Chemical reaction leading to hardening (hydration) needs to be considered by including

dissipative effects [12] also caused by the internal friction during deformation [13,14]. Mechanical response

change related to the loading rate [15,16] or modeling damage in concrete [17–19] has to be modeled and

solved together with the hydration in order to reach accurate prediction of the design-life. Differing bonding

mechanisms between structural parts like concrete and steel is an active research area. In the case of a mechanical

interlocking, this so-called apparent adhesion [20] based on friction contact provides only a limited bond with

respect to a chemical or dispersive adhesion. Amending adhesion is possible by gluing concrete and steel [21–

23] that increases also the time to failure [24]. Such a chemical bond is achieved by a special type of polymer,

specifically in this work, we concentrate on this adhesive material at the interface [25] between concrete and

steel anchor.

Effected by their high stiffness properties, thermosetting polymers are proposed to be used at the interface.

They are already available in the market; mostly an epoxy or vinyl ester is used. By mixing two fluid-like

components, one is a multifunctional comonomer and the other one is the hardener, the chemical reaction called

curing is ignited forming cross-links (chemical bonds) between polymers leading to an increased stiffness. The

curing process never stops until it reaches the fully cured state. This process is of paramount importance for

correctly characterizing thermosets. Until the so-called gel point, the process is fast and causes a significant

shrinkage leading to residual stresses [26], after the gel point we consider the material as a solid—to be precise

a fluid with an infinite viscosity—leading to increased stiffness with further cross-linking. It is possible to

track the curing process experimentally [27] by introducing a degree of curing, experimental determination

of mechanical response [28–30] is an important parameter to calculate the long term response [31]. Mechan-

ical characterization [32] and thermomechanical considerations [33] especially below the gel point [34,35],

or considering aging [36,37] alongside to the curing have been studied. The effect of the degree of cure on

the mechanical response is inevitable. Various models have been used by exploiting numerical multiscale

approaches [38–40] or also simple linear models. Multiphysics approaches with adequate numerical imple-

mentations are developed as well [41], not only for thermosets using small displacement assumption [42,43],

but also for large displacements [44–46], has been applied in several applications [47–49], where the curing

is established in a controlled environment (i.e. heat treatment of composite elements in the autoclave) such

that the material characteristics yield the chemically best configuration. In other words, the polymer attains a

fully cured state and the curing process stops. However, this is rarely possible in case of thermosetting poly-

mers that find application in the construction industry and especially reinforced concrete structures. Owing

to the nature of constructing a structure in an uncontrolled environment under definitely not optimal condi-

tions, the thermoset polymer will typically fail to reach a fully cured state after the installation such that the

curing phenomenon continues during the life-time of the structure. Considering the composition altering due

to the curing and its effect on the design-life of the system is of paramount importance. In order to model

the whole system response accurately, in this work, we start with a thermodynamically consistent modeling

of mechanochemistry in thermosetting polymers. We present a computational implementation with the finite

element method by using the open-source packages developed under the FEniCS project [50,51].

2 Modeling the curing process

The underlying thermosetting polymer is a resin being hardened due to the chemical reaction, called curing,

leading to cross-links between existing polymer chains. Cross-linking affects (often decreases) the volume per

mass, an increase in the overall stiffness, and at the same time altering the temperature due to an exothermal

reaction. After mixing an agent to the resin, curing starts by using energy from the environment—mostly heat

or radiation energy (UV light) is supplied to level off the temperature. Once the curing starts, from the viscous

resin state until the so-called gel point, the curing process is fast. At the gel point the viscosity increases

drastically so that we may call the material a solid and the curing process slows down; its rate asymptotically

reaches zero at the fully cured state. Degree of cure is a difficult variable to measure. We start with a definition

of this variable following [52] by using theory of mixtures in a simplified manner. In a unit volume, V, masses

of resin, curing agent, and cured solid are introduced mR,mA,andmS, respectively. Their values vary in time;

but the sum is constant, m=mR+mA+mS, throughout the curing process. We introduce mass fractions:

YR=mR

m,YA=mA

m,YS=mS

m,(1)

Thermo-mechano-chemical modeling and computation of thermosetting polymers

leading to



Yα=YR+YA+YS=1.(2)

Now we introduce the degree of cure, ω=ω(t), as “measured” in time, giving the mass fraction of the solid

YS=ω, (3)

under the assumption that the initial mass of solid is zero, YS(t=0)=0, such that we have

YR(t=0)+YA(t=0)=1.(4)

Consider that the ratio of masses, YR(t=0)/YA(t=0), remains the same during the reaction. In other words,

same molar mass is used for resin and agent. Now by using z=YR(t=0), we obtain

YR(t)=z(1−ω),

YA(t)=(1−z)(1−ω),

YS(t)=ω.

(5)

Hence, the whole formulation for calculating the masses of constituents has been subsumed to the degree

of cure (conversion degree), ω. For obtaining its value, an evolution equation is developed dependent on the

chemical reaction type. Epoxy resin reaction is based on autocatalysis mechanism [53], effected by hydroxy

groups formed during catalysis. Hence, the following evolution equation [54] is found to be useful

ω•=(k1+k2ωm)(1−ω)n,k1=A1exp −E1

RT ,k2=A2exp −E2

RT ,(6)

where amplitudes, A×, activation energies, E×, and powers, m,n, are constant parameters to be determined.

Universal gas constant, R, is known. Temperature, T, is the key variable altering kinetic rates, k×, and thus

curing rate, ω•, significantly. We refer to [55, Table 1] for parameter values of different materials. Obviously,

ω=ω•dtis used for determining the current degree of cure in each material point. This model fails to

incorporate vitrification such that different mechanisms below and above glass transition temperature, Tg,are

not modeled accurately. Especially curing at low temperatures (below 100◦C) leads to incomplete conversion

that is of paramount interest in fastening systems, where the hardening occurs in environmental conditions. A

possible approach as in [56] characterizes glassy and rubbery states occurring simultaneously. The idea was

introduced in [57] based on phenomenological observation combined with the free volume reduction during

cure leading to decrease in mobility. Mobility plays a dominant role in the rubbery state as molecules collide

and form a network via cross-linking. Beyond a cross-linking density, material vitrifies and in this glassy

state the chemical rate is more dominant. By following [58], we use only primary and secondary amines via

autocatalysis and impurity catalysis (denoted by cin index), collective kinetic rates, K1,K1c, involve diffusion

controlled mechanism, Kdiff, as well as chemically steered mechanism, K1,chem,K1c,chem, as suggested in [59]

as follows:

ω•=K1(1−ω)2ω+K1c

K1,1

K1,chem

Kdiff

K1c

K1c,chem

Kdiff

K1,chem =A3exp −E3

RT ,K1c,chem =A4exp −E4

RT ,Kdiff =kTgexp C1(T−Tg)

C2+|T−Tg|.

(7)

The diffusion rate is modeled by the Williams–Landel–Ferry (WLF) equation based on [60,61]inthe

rubbery regime and greater than chemical rate in many orders in magnitude. In the glassy state, diffusion rate

is nearly zero regarding chemical rate. This interplay delivers a realistic prediction of conversion degree, ω,

at low temperatures, where a so-called partial freezing inhibits to attain ω=1. This model incorporates glass

transition temperature, which is often determined by a fit function [62] based on calorimetric measurements.

For consistency, we follow [58] and use the relation from [63]

Tg=exp (1−ω) ln(Tg,0)+Cωln(Tg,∞)

(1−ω) +Cω,(8)

with C=c∞/c0as the ratio of the heat capacity change at the glass transition temperature of the

fully cured network with ω=1, Tg,∞per the heat capacity change at the glass transition temperature of the

B. E. Abali et al.

Table 1 Parameters for the evolution equation of the conversion degree, ω,giveninEq.(6) (upper part) and in Eq.(7) (lower

part) compiled from [58]

Parameter Variable Value Unit

Amplitude A11.7×1091/s

Amplitude A2280 ×1031/s

Universal gas constant R8.314 J/(mol K) ˆ= G/K

Activation (energy) E1/R12 ×103K

Activation (energy) E2/R7×103K

Power constant m0.6–

Power constant n1.2–

Amplitude A3103.41/s

Activation (energy) E341.5kJ/molˆ= kG

Amplitude A4103.61/s

Activation (energy) E451.8kJ/molˆ= kG

Amplitude kTg3.8×10−51/s

WLF parameter C140 –

WLF parameter C252 K

Ratio C0.57 –

Monomer transition temperature Tg,0-10 ◦C

Fully-cured transition temperature Tg,∞165 ◦C

monomer with ω=0, Tg,0. The necessary coefficients for the evolution equations (6), (7) are determined by

inverse analysis to data obtained by differential scanning calorimetry, we compile them from the literature as

presented in Table 1. By knowing this mass fraction and the current mass density of the mixture, we can extract

masses of each constituent by using Eq.(5) with zbeing the mass fraction of resin initially.

3 Thermodynamical formulation

Under the simplification that deformation is small such that no distinction is necessary between configurations,

we use ()•as the partial time derivative and ,ifor the space derivative with respect to Xi, where space coordinates,

X, denote material particles composing a material system. Hence, the mass balance is satisfied and the mass

density, ρ, is a given function in space—for homogeneous materials we set ρ=const in space, the presented

formalism holds true for heterogeneous materials as well. We need to solve balances of momentum and internal

energy,

ρu••

i−σji,j−ρgi=0,

ρu•+qi,i−ρr=σjiε•

(9)

respectively, where uin m is the displacement to be calculated; stress, σin Pa ˆ=N/m2, heat flux, qin W/m2,

and specific internal energy (energy per mass), uin J/kg, need to be defined; gin N/kg is the gravitational

specific force, and rin W/kg is the given specific supply term. For the sake of simplicity, we use a linear strain

measure:

εij =1

2(ui,j+uj,i). (10)

However, we emphasize that the formulation is also suitable for nonlinear measures. Internal energy needs

to be defined by exploiting thermodynamics. There are several methodologies; we refer to [64]forsuch

examples based on thermodynamics of irreversible processes and non-equilibrium thermodynamics. Briefly,

we introduce the so-called specific Helmholtz free energy:

f=u−Tη, (11)

with the specific entropy, η, to be defined indicating the reversible part of heat flux and temperature, T,tobe

calculated. Free energy is modeled by assuming that it depends on temperature, T,strain,ε, and degree of

cure, ω,

f=fT,ε,ω

.(12)

Thermo-mechano-chemical modeling and computation of thermosetting polymers

We stress that we choose to exclude a dependency on the rate of these variables. Now by inserting the free

energy into the balance of internal energy, dividing by T, and after rewriting, we obtain

ρf•+ρηT•+ρTη•+qi,i−ρr=σjiε•

ij ,

T∂f

∂T+ηT•+ρ

∂f

∂ωω•+ρη•+qi,i

T−ρr

T=1

Tσji −ρ∂f

∂εij ε•

ij ,(13)

Based on the nonequilibrium thermodynamics, we decompose

η=ηeq +ηneq ,σ=σeq +σneq ,(14)

into equilibrium and nonequilibrium parts and then use the well-known relations for equilibrium parts

ηeq =−∂f

∂T,σ

ji =ρ∂f

∂εij

.(15)

For simplicity we assume that rate of temperature and rate of strain fail to alter mechanical response, leading

to ηneq =0, σneq =0. Equation (15)1is often introduced as the definition of the entropy and Eq.(15)2is

Cattaneo’s theorem in its general form for small displacements, we refer to [65] for its derivation from a

Lagrange function. We stress that these relations show that (equilibrium) entropy and (equilibrium) stress

depend on the same arguments as the free energy, η=ηT,ε,ω

and σ=σT,ε,ω

, which is often called

the principle of equipresence [66, §293.η]. By using these relations, we arrive at the balance of entropy:

ρη•+qi

T,i−ρr

T=−qi

T2T,i−ρ

∂f

∂ωω•,(16)

with the right-hand side called the entropy production that is zero for reversible, and it is positive for irreversible

processes. By using this assertion—the second law of thermodynamics—we obtain Fourier’s law and the

restriction

qi=−κT,i,−∂f

∂ωω•≥0,(17)

assuring that the entropy production is (zero or) positive as long as κ≥0. The former is well known in ther-

modynamics. The latter is equivalent to the postulated Karush–Kuhn–Tucker relation originally introduced

in plasticity, for a general discussion we refer to [67]. We emphasize that the chosen evolution equation for ω

in Eq.(6)or(7) is positive such that ∂f/∂ω has to be negative. Heat is inserted into the system in a standard

calorimetric measurement, δQ=cdT, with specific capacity, c. Now by assuming the relation Tdη=δQ,

as proposed in [68] and used in all thermodynamical formulations, we obtain

T=∂η

∂T=− ∂2f

∂T∂T.(18)

All these relations will be fulfilled by defining the free energy adequately. After inserting Eq.(17)1into Eq.(13)2

and using η=ηT,ε,qas well as Eq.(15), we acquire

ρc

TT•+∂η

∂εij

ε•

ij +∂η

∂ωω•−κT,ii

T−ρr

T=−ρ

∂f

∂ωω•,

−ρ∂2f

∂T∂TT•+∂2f

∂εij∂Tε•

ij +∂2f

∂ω∂Tω•−κT,ii

T−ρr

T=−ρ

∂f

∂ωω•.

(19)

This governing equation for the temperature is a general formulation under the assumption that the free energy

depends “only” on temperature, strain, and degree of cure. By defining the free energy adequately, the governing

equation will be closed and then calculated numerically.

B. E. Abali et al.

4 Thermo-mechano-chemistry of thermosets

We aim at calculating the deformation by satisfying

ρu••

i−ρ∂f

∂εij ,j

−ρgi=0,(20)

and the temperature by fulfilling

−ρ∂2f

∂T∂TT•+∂2f

∂εij∂Tε•

ij +∂2f

∂ω∂Tω•−κT,ii

T−ρr

T=−ρ

∂f

∂ωω•,(21)

for a hardening thermosetting polymer, whose state is characterized by the degree of cure, ω, modeled by the

evolution equation (6)or(7). In order to close these governing equations, we need to model the system by

defining a specific free energy, f. The whole formulation is reduced to define a correct model for the free energy.

In the following, we give a possible model and present its relation to the existing models in the literature.

We emphasize that the free energy may depend on temperature, strain, and degree of cure. The real

dependence is possibly determined by experiments; there are no thermodynamical restriction about its form.

But we know, based on the numerical stability and also partly on our understanding of energy as well as

intuition, the energy needs to be positive definite. In full accordance with Eqs.(18), (17)2, we propose to use

the following free energy expression,

f=−cTln T

Tref −Tref+1

2ρεijCijklεkl −(T−Tref)

ρεijCijklαkl −ω

ρεijCijklβkl −ωHref ,(22)

where all material parameters, c,C,α,β,Href, may depend on the same arguments as energy, namely tem-

perature, strain, and degree of cure. This model is the simplest case for a thermoelastic material; the form is

quadratic providing stability. The reference temperature, Tref, is often set as the room temperature under two

assumptions: the material possesses no thermal stresses and the coefficient of thermal expansion (CTE), α,

is measured regarding this temperature. Curing-related shrinkage is given by the material parameter, β. Heat

release because of the exothermic reaction is given by Href, the total heat (per mass) generated by a fully

cured specimen. If the material response is such that c,C,α,βare constant in εand T, in other words we

neglect hyperelasticity as well as temperature-related changes in material coefficients, then we acquire the

(generalized) Hooke’s law with a linear shrinkage and well-known Duhamel- - Neumann extension

σij =ρ∂f

∂εij

=Cijklεkl −(T−Tref)αkl −ωβkl.(23)

Analogously, the specific entropy and its derivative with respect to the strain read

η=−∂f

∂T=cln T

Tref +1

ρεijCijklαkl ,

∂η

∂T=− ∂2f

∂T∂T=c

T,∂η

∂εij

=− ∂2f

∂εij∂T=1

ρCijklαkl .

(24)

In the general case, we expect to have material parameters depending on the cross-linking. A rather obvious

observation is made for stiffness, the material hardens as the cross-linking density increases. Accurate modeling

of the relationship between material parameters and degree of cure needs to be achieved by experiments [69]by

means of detecting ∂f/∂ω in the aforementioned general formulation. A possible way of introducing the degree

of cure is simply to introduce two phases, solid (hardened) and fluid (non-hardened), creating the stiffness in

the following simple rule of mixture: C=(1−ω)Cf+ωCs,α=(1−ω)αf+ωαs,whereCs,αsand Cf,

αfmean the maximum stiffness, maximum CTE reached after full cross-linking; and, initial stiffness, initial

CTE of the non-hardened fluid alike material, respectively. Even the specific heat capacity, c, may depend on

the degree of cure, ω,wedirectto[70] for such a study. The simple relation between material parameters and

degree of cure—motivated by the mixture theory—has a formal benefit of constructing thermodynamically

sound material properties [71, Sect. 4]. Since the parameter for shrinkage, β, is inherently related to the curing

Thermo-mechano-chemical modeling and computation of thermosetting polymers

process, the explicit (linear) dependence is established in the free energy definition instead of introducing

phase depending coefficients. In this manner, we have

∂f

∂ω =1

2ρεijCs

ijkl −Cf

ijklεkl −(T−Tref)

ρεijCs

ijkl −Cf

ijklαkl

−(T−Tref)

ρεijCijklαs

kl −αf

kl−ω

ρεijCs

ijkl −Cf

ijklβkl −1

ρεijCijklβkl −Href ,

∂η

∂ω =− ∂2f

∂ω∂T=− ∂2f

∂T∂ω =1

ρεijCs

ijkl −Cf

ijklαkl +1

ρεijCijklαs

kl −αf

kl,

(25)

such that we reach from Eq.(21)

ρc

TT•+Cijklαklε•

ij +εijCs

ijkl −Cf

ijklαklω•+εijCijklαs

kl −αf

klω•−κT,ii

T−ρr

=− 1

2TεijCs

ijkl −Cf

ijklεklω•+(T−Tref)

TεijCs

ijkl −Cf

ijklαklω•

+(T−Tref)

TεijCijklαs

kl −αf

klω•+1

Tεij2ωCs

ijkl +(1−2ω)Cf

ijklβklω•+ρHref

Tω•,

(26)

leading to

ρc

TT•−κT,ii

T−ρr

T=ρHref

Tω•−ωCs

ijkl +(1−ω)Cf

ijklωαs

kl +(1−ω)αf

klε•

−1

Tεij1

2Cs

ijkl −Cf

ijklεkl +TrefCs

ijkl −Cf

ijklωαs

kl +(1−ω)αf

kl

+TrefωCs

ijkl +(1−ω)Cf

ijklαs

kl −αf

kl−2ωCs

ijkl +(1−2ω)Cf

ijklβklω•.

(27)

Until this stage, we have used the following assumptions:

•empirical fit as in Eq.(6) for the evolution equation of degree of cure,

•small displacements,

•material parameters’ linear dependence on the degree of cure, C=(1−ω)Cf+ωCs,α=(1−ω)αf+ωαs,

•linear material equations provided by the quadratic free energy definition modeled as in Eq. (22).

Especially for modeling fastening systems in concrete, all assumptions are easily justified. Furthermore, for

the sake of a better comprehension, we introduce more simplifications as follows:

1. In the material state before cross-linking, stiffness of the fluid-like or gel-type material vanishes, Cf=0,

and partly as a consequence, a reversible thermal expansion is negligible, αf=0. Then Eq.(27)forthe

temperature is obtained as follows:

ρcT•−κT,ii −ρr=ρHrefω•−Tω2Cs

ijklαs

klε•

ij −εij1

2Cs

ijklεkl +2ωTrefCs

ijklαs

kl −2ωCs

ijklβklω•

=ρHrefω•−Tω2Cs

ijklαs

klε•

ij −1

2εijCs

ijklεklω•+(ω2)•εijCs

ijklβkl −Trefαs

kl.

(28)

2. A further simplification, especially beyond the gel point, relies on the fact that the shrinkage and thermal

expansion are small. By neglecting these effects, we obtain

ρcT•−κT,ii −ρr=ρHref −1

2εijCs

ijklεklω•.(29)

Especially by neglecting deformation, ε=0, we end up with the formulation used in the literature often

for thermal analysis during hardening.

B. E. Abali et al.

3. Further assumptions emerge in calorimetric measurements, without a supply term, r=0. We write the

global form of the latter, after using Gauss- - Ostrogradskiy theorem and Fourier’s law,

B

ρcT•dV−∂B

qinidA=BρHref −1

2εijCs

ijklεklω•dV.(30)

The surface normal, ni, is outward the body such that that the second term is the heat released from the

body. In an experiment, the heat supplied to the body is measured,

Q=−∂B

qinidA.(31)

In the case of an experimental setup, where the specimen is almost free on all sides, we may assume that

no energy is stored. Furthermore, as the specimen size is small, we assume constant temperature within

the body as well as constant rate of curing degree. Heat flow, Q, and temperature, T, are controlled and

measured in the experiment for the specimen of mass, m=BρdV. For a known specific heat capacity—

even if the heat capacity or Hdepends on the temperature, we stress that the temperature is constant within

the body, B—we calculate

cT•+Q

m=Hrefω•.(32)

The mass, m, and mass density, ρ, are constant, and we obtain Href and parameters in the evolution equation

for ω•. Two special cases arise in the measurement. The first case is the isothermal calorimetric measurement,

T•=0, where the specific heat, Q/m, is used for determining Href and parameters for the evolution equation,

ω•. For the former, a full curing state from ω=0toω=1 is accomplished and Href is calculated. Then,

for the latter, a partial curing (even from the same data) is exploited to determine parameters in the chosen

evolution equation. The second case is to perform a dynamic test where T•is set constant in the equipment

in order to characterize the specific heat capacity, c.

Various models in the literature are covered in the methodology presented herein. We continue with the general

governing equations to demonstrate their predictive capability by computing academic examples.

5 Generating the weak form for the finite element method

We consider a general case and set the objective to compute the primitive variables, ω,u,T, in space, x,and

time, t. By starting with the initial conditions of a partly cured body given by ωref as follows:

ω=ωref ,u=0,T=Tref ∀x∈B,t=0,(33)

within the continuum body B, we search for the transient solution of {u,T}, by solving Eqs. (20), (21)and

determine ωby updating

ω=ωref +ω•dt,(34)

with the evolution equation for ω•. We use the quadratic free energy as in Eq.(22), leading to linear material

equations; yet emphasize that the implementation is designed in such a way that more advanced models are

possible to be employed as well. Therefore, we leave the free energy as fin the following. All continuous

functions are approximated by using their discrete representations in space and time. In order to discretize in

time, we use constant time steps and compute in a time series

t={0,t,2t,3t,...}(35)

with (·)0,(·)00 denoting the calculated values of one, two time steps before, respectively. Time derivative is

approximated by the finite difference method (Euler backwards scheme)

∂(·)

∂t=(·)−(·)0

t.(36)

Thermo-mechano-chemical modeling and computation of thermosetting polymers

The governing equations become

ρui−2u0

i+u00

t2−ρ∂f

∂εji ,j

−ρgi=0

−ρ∂2f

∂T∂T

T−T0

t−ρ∂2f

∂εij∂T

εij −ε0

t−ρ∂2f

∂ω∂Tω•−κT,ii

T−ρr

T=−ρ

∂f

∂ωω•

ω:= ω+ω•t.

(37)

We remark the linearized strain measure in order to find out the value one time step before,

εij =1

2(ui,j+uj,i), ε

ij =1

2(u0

i,j+u0

j,i). (38)

For space discretization, we use the finite element method with a triangulation of the continuum space, B,and

exchanging continuous primitive functions with their counterparts in a Hilbertian Sobolev space [72]by

using nth polynomial degree form functions for each discrete element,

V={u,T}∈[Hn(B)]4:{u,T}∂B=given.(39)

For the sake of notational easiness, we simply skip to distinguish between continuous functions and their

discrete representations. We use the standard procedure to construct the integral forms by building residuals

from governing equations and multiplying (weighting) them by corresponding test functions, δu,δT,where

they are expressed in the same space as the primitive variables, called the Galerkin procedure,

Bρui−2u0

i+u00

t2−ρ∂f

∂εji ,j

−ρgiδuidV=0

B−ρ∂2f

∂T∂T

T−T0

t−ρ∂2f

∂εij∂T

εij −ε0

t−ρ∂2f

∂ω∂Tω•−κT,ii

T−ρr

T+ρ

∂f

∂ωω•δTdV=0.

(40)

We emphasize that stress already possesses space derivative of primitive variables, because of the additional

divergence on it, we need to choose at least n=2 in the space discretization. This continuity condition is

weakened by integrating by parts for all second derivative terms,

Bρui−2u0

i+u00

t2δui+ρ∂f

∂εji

δui,j−ρgiδuidV−∂B

tiδuidA=0

B−ρ∂2f

∂T∂T

T−T0

t

δT−ρ∂2f

∂εij∂T

εij −ε0

t

δT−ρ∂2f

∂ω∂Tω•δT+κT,iδT

T,i

−ρr

δT+ρ

∂f

∂ωω•δTdV+∂B

ˆq

δT

TdA=0.

(41)

wherewehaveinsertedˆ

ti=njσji =nj∂f/∂εij and ˆq=niqi=−niκT,ion boundaries. For the displacement

with a given traction, ˆ

tin Pa, two sets of boundaries appear: Dirichlet boundaries, ∂BD,andNeumann

boundaries, ∂BN. We restrict the formulation such that these two sets fail to intersect, ∂BD∩∂BN= {},so

the whole boundary is the union of them, ∂B=∂BD∪∂BN.AttheDirichlet boundary, the displacement

itself is given, hence, no test functions are necessary, δu=0∀x∈∂BD. For the rest, traction is given,

t=given ∀x∈∂BN. Analogously, the temperature may be modeled by given temperature or heat flux;

however, we use a more natural approach by using a mixed boundary condition called Robin boundaries,

herein applied to the whole boundary,

ˆq=h(T−Ta)∀x∈∂B,(42)

B. E. Abali et al.

with the convection parameter, h, and the ambient (environment) temperature we use as Ta=Tref. Two weak

forms are in different units, the former is in the unit of energy and the latter is in the unit of power. After

multiplying by the constant time step t,

Fu=Bρui−2u0

i+u00

t2δui+ρ∂f

∂εji

δui,j−ρgiδuidV−∂B

tiδuidA,

FT=B−ρ∂2f

∂T∂T(T−T0)δT−ρ∂2f

∂εij∂T(εij −ε0

ij)δT−tρ∂2f

∂ω∂Tω•δT+tκT,iδT

T,i

−tρr

δT+tρ

∂f

∂ωω•δTdV+∂B

th(T−Tref)

δT

TdA.

(43)

we can sum them up

Form =Fu+FT(44)

and use for computing thermomechanics of a curing thermoset.

6 Computational implementation and numerical examples

In order to examine the strength and accuracy in predicting the behavior of realistic systems, we implement the

weak form in Eq.(44) by using open source packages developed within the FEniCS project. The formulation

is implemented in Python, wrapped into C++, compiled, and solved by mpirun in parallel in a Linux computer

(Ubuntu 18.04). This scalable implementation is suitable for larger systems as well. Herein we demonstrate sim-

ple examples to comprehend the multiphysics in thermo-mechano-chemical systems. Especially two important

features need to be emphasized leading to the implementation presented herein.

First, the weak form is nonlinear such that we use the conventional Newton- - Raphson linearization

approach for solving a linear problem and update the solution by starting from the last converged solution.

Even if the system is governed by stiff differential equations, using small time steps leads to a convergence

in each iteration. The linearization is handled by exploiting symbolic differentiation that allows to implement

also highly nonlinear forms—for example a hyperelastic material model can be easily implemented with the

approach developed herein, we refer to [64] as well as to the supply code of [65] for examples of such an

implementation.

Second, the specific free energy in Eq.(22) has been implemented also with the aid of symbolic differ-

entiation. Technically, the whole material formulation is reduced to the free energy definition such that the

implementation is very general and can be easily adapted for other systems simply by changing the free

energy. We stress that the material modeling is limited by the assumed free energy formulation; however,

implementation is novel in the sense that any more advanced formulations are possible to utilize. We make the

code publicly available in [73] for encouraging its use under GNU Public licensing [74] for an efficient and

transparent scientific exchange.

6.1 Curing models comparison

Two evolution equations (6), (7) have been introduced and implemented. Both equations consist of parameters

depending on the current temperature such that the weak form FTin Eq. (43) is nonlinear. Material parameters

belong to typical epoxy type of materials with coefficients compiled in Tables 1,3. We compute temperature

anddegreeofcureforasimplegeometryof50mm×10mm×10mm, made of homogeneous material. First,

we solve the temperature distribution by using model A in Eq. (6), and secondly, we obtain the temperature

by applying model B in Eq.(7). Three-dimensional simulation results are recorded, and both results at the

same location (geometric midpoint of the domain) are compared qualitatively. We emphasize that the material

parameters in Table 1are realistic; but from different sources such that models A and B fail to be compared

quantitatively.

We begin with the well-known effect of the ambient temperature on the curing phenomenon. By setting the

convection parameter high, h=1kW/(m2K), we enforce nearly a constant surface temperature equivalent to

the ambient temperature, Tamb. Applying model A and model B leads to results in Fig.1. As expected, higher

temperatures increase the conversion rate such that the possible upper limit is attained in less than 90 minutes

Thermo-mechano-chemical modeling and computation of thermosetting polymers

Fig. 1 Degree of cure, ω, over time, t, under different ambient temperatures computed by the model A (left) as in Eq. (6)and

model B (right) as in Eq. (7)

Fig. 2 Heat flow in a temperature ramp with two evolution equations: left model A in Eq.(6) and right model B in Eq.(7)

at 90◦C than in 20 hours at 50◦C. These simulations model an isotherm calorimetry measurement. Model A is

capable of quantifying the process only in the case of complete curing. We point out the case of an incomplete

curing taking place at low temperatures. The model B given by the evolution equation (7) is indeed capable of

simulating this freezing behavior effected by the vitrification. Such a vitrification, 10–20 ◦C above the curing

temperature, is typically observed in isotherm calorimetry measurements. This glass transition temperature,

Tg, is also detected and measured with the aid of this observation. During cross-linking, the exotherm reaction

produces heat which is taken out of the system in order to hold the temperature constant. For simplicity,

we use a constant value for the specific heat capacity in simulations. In reality, the change of specific heat

capacity between glassy and rubbery states causes an overshoot in this heat energy, which is used to determine

the numerical value of Tgin a non-fully cured structure. Material showing a distinct Tgjustifies the use of

model B. A typical example are the epoxy systems used in post-installed fastening applications subject to low

(ambient) temperature curing.

For a better understanding of the role of the evolution equation on the released heat, we simulate a tem-

perature ramp test, where structures with different degrees of cure are simulated as being in a temperature

chamber where the ambient temperature increases by 10K per minute. Figure 2demonstrates the specific heat

flow, Q/m, obtained by Eq.(32) from the temperature and degree of cure computation. Obviously, model A is

symmetric leading to the sigmoid relation, whereas model B shows a kink at the glass transition temperature

such that the symmetric character is lost.

In general, the computational implementation results in reliable simulations with realistic parameters. Such

a tool is of importance to study and comprehend the behavior of different curing models. We have selected

two models without assessment of the quantitative accuracy to a specific material.

B. E. Abali et al.

Table 2 Used temperature alternation on one day approximated to the data for Brussels (BE) on October 1, 2019

Time Temperature in ◦C

12 am 16

2am 14

9am 17

12 pm 15

5pm 20

12 am 15

Fig. 3 For a comparative analysis, curing of polymer adhesive material as used in Brussels (Belgium) on October 1, 2019. Left:

Using the temperature change over the day by using 5 picked values as given in Table 2. Right: Mean temperature value applied

on the boundary and shown in the core of the sample

6.2 Curing under alternating thermal loading

For a realistic simulation, we use model B and consider that the epoxy with parameters in Table 1has been

used for adhering components in an outdoor facility. Complete curing of the material is desired for an ade-

quate stiffness leading to high adhering. We obtain hourly temperature, as an example on October 1, 2019

in Brussels,1and use with a rough approximation as given in Table 2. By using a convection parameter

h=100W/(m2K), we model a mild windy ambient for heat exchange over the surface. Again for a simple

geometry of 50mm×10mm×10 mm, we compute transiently the temperature distribution in tree-dimensional

continuum coupled to the degree of cure given by the model B. The simulation starts directly after casting,

ω=0, computes for 24 hours temperature and degree of cure. The temperature values in the core of the

structure are recorded, which are nearly the same as the surface temperature since the structure is small. We

emphasize that ambient temperature and structure’s temperature are not identical. Convection over the surface

as well as curing produced heat are responsible for this difference.

Important observations are attained from results demonstrated in Fig. 3. After one day of curing in envi-

ronmental temperatures, only ω=0.35 has been reached in Fig. 3(left). With high mobility, rate of conversion

is accelerating in the beginning, especially after the early morning when the temperature starts rising. As a

coincidence, on this day the temperature drops before noon such that the material presumably vitrifies and

afternoon temperature increase is only partly amending the conversion rate. If the material has a gel point of

ω=0.5, then 24 hours of curing in a mild autumn day reaches not even beyond the gel point. Obviously, the

expected stiffness is not delivered in this condition, at least for this material. Therefore, different commercial

products use modified compositions, allowing an increased reaction kinetics in these temperatures. By utilizing

the implementation herein, we enable investigating the degree of cure of a material with known parameters.

Moreover, we stress that the temperature history is of importance as well. When we use the mean value for

the ambient temperature, instead of the alternating temperature, we observe that a significant difference of

approximately 10% is captured in Fig.3(right). Even a linear relation between stiffness and the degree of cure

is assumed, this 10% deviation is of importance to consider in design. This analysis is for a small sample, where

the temperature conductivity is high enough to generate a homogeneous temperature distribution. However, the

exothermal reaction is still altering the temperature in the middle of this sample as visible in Fig.3(right). But

1Taken from https://www.wunderground.com/dashboard/pws/IBRUSS7/graph/2019-10-1/2019-10-1/daily.

Thermo-mechano-chemical modeling and computation of thermosetting polymers

this change is insignificant. Herein, we demonstrate how valuable such a study becomes especially whenever

applied to more realistic geometries.

6.3 Mechanical loading superposed by hardening

Under a mechanical loading, the amount of deformation depends on the stiffness that changes during hardening.

Two possible implementations can be developed. One approach is to use a stiffness, Cs, for the fluid phase

leading to a soft matter even before the gel point. Then we start with ω=0 and with the aid of this initial

stiffness, the numerical implementation converges. Another approach is to set the fluid stiffness to zero,

Cs=0, αs=0, and start with an initial curing, ωref, greater than zero such that the simulation begins

with an initial stiffness again leading to a convergence. We choose the second approach, for example the gel

point, ωref =0.5, where the material response is more solid type than fluid. We use the latter approach in

the following simulations. Materials response is chosen to be elastic and the initial geometry is used as the

reference configuration to which the continuum body recovers after unloading. This modeling is based on

the fact that the hardening beyond ωref fails to affect the reference configuration. We leave an experimental

verification of this assumption to further research.

A simple, thick beam of 50mm×10mm×10mm is modeled. Material is assumed to be isotropic such that

the solid properties are given by Lame’s constants, λ,μ, coefficient of thermal expansion, α, and shrinkage

parameter, β,

ijkl =λδijδkl +μδikδjl +μδilδjk ,α

ij =αδij ,β

ij =βδij .(45)

The Lame parameters are related to the (measurable) engineering constants, E,ν, as follows:

λ=Eν

(1+ν)(1−2ν) ,μ=E

2(1+ν) .(46)

We use fictitious but realistic material coefficients as compiled in Table 3. Nearly all coefficients can be

determined by experiments on the corresponding material. But the thermal convection parameter modeling the

rate of heat exchange between the surface of the material and environment depends on the material, surface

roughness, and ambient state [75]. For example, moving air is known to extract more heat from the material

than still air. Analogously, humid air possesses a higher heat transfer coefficient than dry air. We consider

a uniaxial tensile test, the beam is clamped on one end and pulled on the other end with a controlled force,

first, relatively quickly until it reaches its maximum and then, secondly, it is held at this force throughout the

whole simulation of 10 minutes. In the case of an elastic material, the response is instantaneous, the beam

gets stretched and remains at the same length under constant force. This case occurs in Fig.4(top), where the

ambient temperature is less than Tgsuch that the material remains at the same degree of cure, ωref.Inthe

case of an ambient temperature higher than the glass transition temperature, with the additional assumption

that Young’s modulus remains the same, response to the applied force is the same; however, the material

undergoes further curing. Hence, the material hardens, stiffness increases, and the beam gets shorter under the

same force held constantly, see Fig.4(bottom). Since we have modeled stiffness to degree of cure linearly, the

length change is linear as well. In reality, Young’s modulus depends on the temperature, especially around the

glass transition. The resulting decrease of the modulus is significant. Herein we neglect the viscous character

and model this behavior by using inverse tangents function and constant moduli for rubbery and glassy states,

Erand Eg, respectively, as follows:

E=Eg−Er

21+2

πarctan(Tg−T)+Er.(47)

For the simulation, we use again fictitious but realistic parameters, Eg=2GPa and Er=0.5GPa, leading

to a stark decrease at the Tg, for example as in Fig.4. We emphasize that the glass transition temperature,

Tg, increases with evolving degree of cure given in Eq.(8). Therefore, we repeat the last analysis and apply

a given boundary condition, hold this displacement, increase the ambient temperature, measure the necessary

force. In the case of constant Young’s modulus, as we model (linear) elastic material, the response remains

the same as long as no postcuring occurs. Since Young’s modulus decreases at the time crossing the current

glass transition temperature—we stress that the simulation starts off with ωref—stiffness drops because of

decreasing Young’s modulus. Furthermore, curing begins and material hardens as well. Such a response is

B. E. Abali et al.

Table 3 Coefficients and material parameters used in the simulation with SI units

Parameter Variable Value Unit

Mass density ρ2400 kg/m3

Thermal conductivity κ2.3 J/(smK)

Specific heat capacity c1100 J/(kgK)

Coefficient of thermal expansion α20 ×10−61/K

Coefficient of shrinkage β2×10−5–

Thermal convection parameter h2–20 J/(sm2K)

Heat release Href 300 J/g

Young’s modulus E2000 MPa

Poisson’s ratio ν0.2–

Fig. 4 Uniaxial test, (engineering) stress, σ11,degreeofcure,ω, engineering strain, ε11, and glass transition temperature, Tg,for

a given ambient temperature. Top: Tamb =30 ◦C. Bottom: Tamb =60 ◦C

Fig. 5 Young’s modulus change during the transition from glassy to the rubbery state, modeled by an arctan function for the

case of Tg=50◦C

Thermo-mechano-chemical modeling and computation of thermosetting polymers

Fig. 6 Stress relaxation under constant stretching in three regimes: first, the modulus decreases by crossing to the rubbery state,

second, postcuring hardens the material and stiffness increases, third, modulus increases by crossing to the glassy state

Table 4 Variation of convection parameter, h, and the computed time necessary to attain the state with 95% of cured solid material

hin W/(m2K) 2.5 5 10 20

tin min for ω=0.95 51.7 58.6 61.7 63.1

experimentally observed, we refer to [76] for a discussion in this sense. Herein, we observe analogous response

in Fig.6.

The curing rate depends on the local temperature that converges to the ambient temperature steered by the

convection parameter, h. Even for the same ambient temperature Tref =300 K, a variation of the convection

parameter shows that the increased heat exchange, larger hin the simulation, increases the time necessary for

attaining the same degree of cure. This phenomenon is demonstrated in Table4quantitatively, where the time

is determined, at which, in the middle of the same structure, ω=0.95 has been reached. The justification of

this phenomenon is as follows: Increased heat exchange leads to the increasing energy loss to the environment

in form of heat energy that is generated as a consequence of the exothermal curing reaction. Obviously, smaller

heat transport parameter represents a better isolation from the environment, practically sealing the structure.

The generated heat during cross-linking is captured within the material, as seen in the temperature evolution

at the middle point of the beam demonstrated in Fig.7. In other words, the heat energy, which is released free

at the moment of curing, increases locally the temperature and accelerates the curing process. The evolution

equation for curing depends on the temperature.

7Conclusion

A general framework has been developed for thermomechanics of thermosetting polymers under curing gov-

erned by an evolution equation. The numerical implementation of such a thermo-mechano-chemical process

has been established by using open-source packages. Simple geometries are utilized for promoting a better

understanding of underlying multiphysics as well as capability of the computational implementation. The

following assumptions have been undertaken:

•Small strain formulation is used; there are several applications especially in composite materials, where

this assumption is accurate.

•Dissipative parts in the thermal and mechanical terms are neglected. Depending on the chosen material, a

viscoelastic term might be necessary.

•The evolution equation for curing depends solely on the temperature. The temperature dependency is

intuitive; however, it is not known, if a mechanical deformation would increase the rate of degree of cure.

•The reference configuration is chosen at the gel point, which has been the initial condition as well. Although

this assumption is very tempting, we leave an experimental validation of this assumption to further research.

Thermodynamically sound modeling of thermo-mechano-chemical processes allows us to accurately model the

coupling effects between variables like temperature, deformation, and degree of cure. An important contribution

is that we have subsumed all modeling by using the free energy such that any possible restrictions of the

B. E. Abali et al.

Fig. 7 Local temperature evolution for different cases by varying heat transfer coefficient for the same ambient temperature

simulation is recovered by revisiting the free energy definition. The computational implementation has been

utilized by using an automatized generation of relations with the aid of symbolic differentiation of the assumed

free energy. Therefore, any more advanced formulation is easily implemented into the framework developed

herein. For demonstrating the strength of the numerical implementation as well as comprehending the material

model’s capabilities, three different simulations have been utilized:

1. Two evolution equations, model A and B, are compared by means of the released heat in order to present

the kink in the evolution denoting the glass transition temperature.

2. A realistic scenario has been used for demonstrating the vitrification in the case of an uncontrolled ambient

temperature like the daily alternating temperature in Brussels.

3. Superposed hardening and mechanical loading is simulated during the cross-linking related glass transition

change leading to a softening. This experimentally observed interplay between hardening and softening

because of thermal and chemical contribution is possible by the proposed multiphysics simulation by adding

glass transition dependence to the Young’s modulus. Additionally, the role of the heat convection across

the boundaries has been discussed.

We stress that the code is publicly available in [73] for encouraging its use under GNU Public licensing [74].

Acknowledgements This research was funded by the Christian Doppler Forschungsgesellschaft, Grant LICROFAST to Roman

Wan-Wendner. The financial support by the Austrian Federal Ministry for Digital and Economic Affairs and the National Foun-

dation for Research, Technology, and Development is gratefully acknowledged, as is the partial financial support by the Austrian

Science Fund (FWF) under Grant P 31203-N32. Jan Vorel acknowledges the financial support provided by the Czech Technical

University in Prague within SGS project with the application registered under the No. SGS20/155/OHK1/3T/11.

Funding Open Access funding enabled and organized by Projekt DEAL.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use,

sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original

author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other

third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit

line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted

Thermo-mechano-chemical modeling and computation of thermosetting polymers

by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To

view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References

1. Hong, H.: Assessment of reliability of aging reinforced concrete structures. J. Struct. Eng. 126(12), 1458–1465 (2000)

2. Wang, C., Li, Q., Ellingwood, B.R.: Time-dependent reliability of ageing structures: an approximate approach. Struct.

Infrastruct. Eng. 12(12), 1566–1572 (2016)

3. Wan, L., Wendner, R., Liang, B., Cusatis, G.: Analysis of the behavior of ultra high performance concrete at early age. Cem.

Concr. Compos. 74, 120–135 (2016)

4. Czernuschka, L., Nincevic, K., Boumakis, I., Wan-Wendner, L., Wan-Wendner, R.: Aging behavior of normal and high

strength concretes. In: Computational Modelling of Concrete Structures: Proceedings of the Conference on Computational

Modelling of Concrete and Concrete Structures (EURO-C 2018), February 26-March 1, 2018, Bad Hofgastein, Austria, p.

197, CRC Press (2018)

5. González, A., Barrias, A., Teixeira, R., Martinez, D., Heitner, B., Antonopoulou, S., Zou, G., Gonzalez Merino, A., Sourav,

S.N.A., Casas, J.R., et al.: Truss, a european innovative training network dealing with the challenges of an aging infrastructure

network. In: The 7th Asia-Pacific Workshop on Structural Health Monitoring (APWSHM 2018), Hong Kong, China, 12–15

November 2018 (2018)

6. Kim, S.-G., Park, Y.-S., Lee, Y.-H.: Rate-type age-dependent constitutive formulation of concrete loaded at an early age.

Materials 12(3), 514 (2019)

7. Ulm, F.-J., Constantinides, G., Heukamp, F.: Is concrete a poromechanics materials? A multiscale investigation of poroelastic

properties. Mater. Struct. 37(1), 43–58 (2004)

8. Aigner, E., Lackner, R., Pichler, C.: Multiscale prediction of viscoelastic properties of asphalt concrete. J. Mater. Civil Eng.

21(12), 771–780 (2009)

9. Unger, J.F., Eckardt, S.: Multiscale modeling of concrete. Arch. Comput. Methods Eng. 18(3), 341 (2011)

10. Contrafatto, L., Cuomo, M., Greco, L.: Meso-scale simulation of concrete multiaxial behaviour. Eur. J. Environ. Civil Eng.

21(7–8), 896–911 (2017)

11. Giorgio, I., Scerrato, D.: Multi-scale concrete model with rate-dependent internal friction. Eur. J. Environ. Civil Eng. 21(7–8),

821–839 (2017)

12. Scerrato, D., Giorgio, I., Della Corte, A., Madeo, A., Dowling, N., Darve, F.: Towards the design of an enriched concrete

with enhanced dissipation performances. Cem. Concr. Res. 84, 48–61 (2016)

13. Scerrato, D., Giorgio, I., Della Corte, A., Madeo, A., Limam, A.: A micro-structural model for dissipation phenomena in the

concrete. Int. J. Numer. Anal. Methods Geomech. 39(18), 2037–2052 (2015)

14. Alam, S.Y., Loukili, A.: Transition from energy dissipation to crack openings during continuum-discontinuum fracture of

concrete. Int. J. Fract. 206(1), 49–66 (2017)

15. Chiaia, B., Kumpyak, O., Placidi, L., Maksimov, V.: Experimental analysis and modeling of two-way reinforced concrete

slabs over different kinds of yielding supports under short-term dynamic loading. Eng. Struct. 96, 88–99 (2015)

16. dell’Isola, F., Bragov, A.M., Igumnov, L.A., Abali, B.E., Lomunov, A.K., Lamzin, D.A., Konstantinov, A.Y.: Mechanical

Response Change in Fine Grain Concrete Under High Strain and Stress Rates, New Achievements in Continuum Mechanics

and Thermodynamics, pp. 71–80. Springer, Berlin (2019)

17. Van Mier, J.G.: Concrete Fracture: A Multiscale Approach. CRC Press, Cambridge (2012)

18. Contrafatto, L., Cuomo, M., Gazzo, S.: A concrete homogenisation technique at meso-scale level accounting for damaging

behaviour of cement paste and aggregates. Comput. Struct. 173, 1–18 (2016)

19. Wan-Wendner, L., Wan-Wendner, R., Cusatis, G.: Age-dependent size effect and fracture characteristics of ultra-high per-

formance concrete. Cem. Concr. Compos. 85, 67–82 (2018)

20. Sauer, R.A.: A survey of computational models for adhesion. J. Adhes. 92(2), 81–120 (2016)

21. Boumakis, I., Marcon, M., Nincevic, K., Czernuschka, L.-M., Wan-Wendner, R.: Concrete creep effect on bond stress in

adhesive fastening systems. In: 3rd International Symposium on Connections between Steel and Concrete, pp. 396–406

(2017)

22. Podroužek, J., Vorel, J., Wan-Wendner, R.: Design for lifecycle robustness of fastening systems. Beton-und Stahlbetonbau

113, 62–66 (2018)

23. Wan-Wendner, R., Podroužek, J.: Robust power law extrapolation for adhesive anchors under sustained load. ACI Struct. J.

116(1), 71–81 (2019)

24. Ninˇcevi´c, K., Boumakis, I., Meissl, S., Wan-Wendner, R.: Consistent time-to-failure tests and analyses of adhesive anchor

systems. Appl. Sci. 10(4), 1527 (2020)

25. Marcon, M., Vorel, J., Ninˇcevi´c, K., Wan-Wendner, R.: Modeling adhesive anchors in a discrete element framework. Materials

10(8), 917 (2017)

26. Kiasat, M.: Curing shrinkage and residual stresses in viscoelastic thermosetting resins and composites, Ph.D. thesis, Technical

University of Delft (2000)

27. Yagimli, B., Lion, A.: Experimental investigations and material modelling of curing processes under small deformations.

ZAMM-J. Appl. Math. Mech./Zeitschrift für Angewandte Mathematik und Mechanik 91(5), 342–359 (2011)

28. Liao, Z., Hossain, M., Yao, X., Mehnert, M., Steinmann, P.: On thermo-viscoelastic experimental characterisation and

numerical modelling of VHB polymer. Int. J. Non-Linear Mech. 118, 103263 (2019)

29. Singer, G., Sinn, G., Lichtenegger, H.C., Veigel, S., Zecchini, M., Wan-Wendner, R.: Evaluation of in-situ shrinkage and

expansion properties of polymer composite materials for adhesive anchor systems by a novel approach based on digital image

correlation. Polymer Test. 79, 106035 (2019)

30. Bradler, P.R., Fischer, J., Wan-Wendner, R., Lang, R.W.: Shear test equipment for testing various polymeric materials by

using standardized multipurpose specimens with minor adaptions. Polymer Test. 75, 93–98 (2019)

B. E. Abali et al.

31. Fischer, J., Bradler, P.R., Schmidtbauer, D., Lang, R.W., Wan-Wendner, R.: Long-term creep behavior of resin-based polymers

in the construction industry. Mater. Today Commun. 18, 60–65 (2019)

32. Hossain, M., Steinmann, P.: Modelling and simulation of the curing process of polymers by a modified formulation of the

arruda-boyce model. Arch. Mech. 63(5–6), 621–633 (2011)

33. André, M., Wriggers, P.: Thermo-mechanical behaviour of rubber materials during vulcanization. Int. J. Solids Struct.

42(16–17), 4758–4778 (2005)

34. Mahnken, R.: Thermodynamic consistent modeling of polymer curing coupled to visco-elasticity at large strains. Int. J.

Solids Struct. 50(13), 2003–2021 (2013)

35. Dal, H., Zopf, C., Kaliske, M.: Micro-sphere based viscoplastic constitutive model for uncured green rubber. Int. J. Solid.

Struct. 132, 201–217 (2018)

36. Johlitz, M., Lion, A.: Chemo-thermomechanical ageing of elastomers based on multiphase continuum mechanics. Continuum

Mech. Thermodyn. 25(5), 605–624 (2013)

37. Nguyen, T.-T., Waldmann, D., Bui, T.Q.: Phase field simulation of early-age fracture in cement-based materials. Int. J. Solids

Struct. 191, 157 (2019)

38. Klinge, S., Bartels, A., Steinmann, P.: The multiscale approach to the curing of polymers incorporating viscous and shrinkage

effects. Int. J. Solids Struct. 49(26), 3883–3900 (2012)

39. Klinge, S., Hackl, K.: Application of the multiscale fem to the modeling of nonlinear composites with a random microstructure.

Int. J. Multiscale Comput. Eng. 10(3), 36–391 (2012)

40. Otero, F., Oller, S., Martinez, X.: Multiscale computational homogenization: review and proposal of a new enhanced-first-

order method. Arch. Comput. Methods Eng. 25(2), 479–505 (2018)

41. Zohdi, T., Wriggers, P.: Phenomenological modeling and numerical simulation of the environmental degradation of multi-

phase engineering materials. Arch. Appl. Mech. 70(1–3), 47–64 (2000)

42. Hossain, M., Possart, G., Steinmann, P.: A small-strain model to simulate the curing of thermosets. Comput. Mech. 43(6),

769–779 (2009)

43. Hossain, M., Steinmann, P.: Degree of cure-dependent modelling for polymer curing processes at small-strain. Part I:

consistent reformulation. Comput. Mech. 53(4), 777–787 (2014)

44. Hossain, M., Possart, G., Steinmann, P.: A finite strain framework for the simulation of polymer curing. Part I: elasticity.

Comput. Mech. 44(5), 621–630 (2009)

45. Hossain, M., Possart, G., Steinmann, P.: A finite strain framework for the simulation of polymer curing. Part II. Viscoelasticity

and shrinkage. Comput. Mech. 46(3), 363–375 (2010)

46. Sain, T., Loeffel, K., Chester, S.: A thermo-chemo-mechanically coupled constitutive model for curing of glassy polymers.

J. Mecha. Phys. Solids 116, 267–289 (2018)

47. Landgraf, R., Rudolph, M., Scherzer, R., Ihlemann, J.: Modelling and simulation of adhesive curing processes in bonded

piezo metal composites. Comput. Mech. 54(2), 547–565 (2014)

48. Hu, S., Chen, Y.: A fem coupling model for properties prediction during the curing of an epoxy adhesive for a novel assembly

of radio telescope panel, Advances in Optical and Mechanical Technologies for Telescopes and Instrumentation, vol. 9151,

p. 915132. International Society for Optics and Photonics (2014)

49. Leistner, C., Hartmann, S., Abliz, D., Ziegmann, G.: Modeling and simulation of the curing process of epoxy resins using

finite elements. Continuum Mech. Thermodyn. 32(2), 327–350 (2020)

50. Alnaes, M.S., Logg, A., Mardal, K.A., Skavhaug, O., Langtangen, H.P.: Unified framework for finite element assembly. Int.

J. Comput. Sci. Eng. 4(4), 231–244 (2009)

51. Logg, A., Mardal, K.A., Wells, G.: Automated Solution of Differential Equations by the Finite Element Method: The FEniCS

Book. Springer, Berlin (2012)

52. Lion, A., Höfer, P.: On the phenomenological representation of curing phenomena in continuum mechanics. Arch. Mech.

59(1), 59–89 (2007)

53. Shechter, L., Wynstra, J., Kurkjy, R.P.: Glycidyl ether reactions with amines. Ind. Eng. Chem. 48(1), 94–97 (1956)

54. Kamal, M.R.: Thermoset characterization for moldability analysis. Polym. Eng. Sci. 14(3), 231–239 (1974)

55. Heinrich, C., Aldridge, M., Wineman, A.S., Kieffer, J., Waas, A.M., Shahwan, K.W.: Generation of heat and stress during

the cure of polymers used in fiber composites. Int. J. Eng. Sci. 53, 85–111 (2012)

56. Flammersheim, H.-J., Opfermann, J.R.: Investigation of epoxide curing reactions by differential scanning calorimetry-formal

kinetic evaluation. Macromol. Mater. Eng. 286(3), 143–150 (2001)

57. Chern, C.-S., Poehlein, G.W.: A kinetic model for curing reactions of epoxides with amines. Polym. Eng. Sci. 27(11),

788–795 (1987)

58. Wise, C.W., Cook, W.D., Goodwin, A.A.: Chemico-diffusion kinetics of model epoxy-amine resins. Polymer 38(13), 3251–

3261 (1997)

59. Rabinowitch, E.: Collision, co-ordination, diffusion and reaction velocity in condensed systems. Trans. Faraday Soc. 33,

1225–1233 (1937)

60. Doolittle, A.K.: Studies in Newtonian flow. II. The dependence of the viscosity of liquids on free-space. J. Appl. Phys. 22(12),

1471–1475 (1951)

61. Williams, M.L., Landel, R.F., Ferry, J.D.: The temperature dependence of relaxation mechanisms in amorphous polymers

and other glass-forming liquids. J. Am. Chem. Soc. 77(14), 3701–3707 (1955)

62. Hesekamp, D., Broecker, H.C., Pahl, M.H.: Chemo-rheology of cross-linking polymers. Chem. Eng. Technol. Ind. Chem.

Plant Equip. Process Eng. Biotechnol. 21(2), 149–153 (1998)

63. Venditti, R., Gillham, J.: A relationship between the glass transition temperature (tg) and fractional conversion for thermoset-

ting systems. J. Appl. Polym. Sci. 64(1), 3–14 (1997)

64. Abali, B.E.: Computational Reality, Solving Nonlinear and Coupled Problems in Continuum Mechanics, vol. 55 of Advanced

Structured Materials. Springer Nature, Singapore (2017)

65. Abali, B.E., Müller, W.H., dell’Isola, F.: Theory and computation of higher gradient elasticity theories based on action

principles. Arch. Appl. Mech. 87(9), 1495–1510 (2017)

Thermo-mechano-chemical modeling and computation of thermosetting polymers

66. Truesdell, C., Toupin, R.A.: Principles of classical mechanics and field theory, Handbuch der Physik Vol. III/1 (Ed. by Flügge,

S.) (1960)

67. Placidi, L., Barchiesi, E.: Energy approach to brittle fracture in strain-gradient modelling. Proc. R. Soc. A Math. Phys. Eng.

Sci. 474(2210), 20170878 (2018)

68. Caratheodory, C.: Untersuchungen über die Grundlagen der Thermodynamik. Mathematische Annalen 67(3), 355–386 (1909)

69. O’Brien, D.J., Mather, P.T., White, S.R.: Viscoelastic properties of an epoxy resin during cure. J. Compos. Mater. 35(10),

883–904 (2001)

70. Chern, B.-C., Moon, T.J., Howell, J.R., Tan, W.: New experimental data for enthalpy of reaction and temperature- and

degree-of-cure-dependent specific heat and thermal conductivity of the Hercules 3501–6 epoxy system. J. Compos. Mater.

36(17), 2061–2072 (2002)

71. Landgraf, R.: Modellierung und Simulation der Aushärtung polymerer Werkstoffe, Ph.D. thesis, Technical University of

Chemnitz (2016)

72. Zohdi, T.I.: Finite Element Primer for Beginners. Springer, Berlin (2018)

73. Abali, B.E.: Supply code (2019). http://bilenemek.abali.org

74. GNU Public: Gnu general public license (2007). http://www.gnu.org/copyleft/gpl.html

75. Still, M., Venzke, H., Durst, F., Melling, A.: Influence of humidity on the convective heat transfer from small cylinders. Exp.

Fluids 24(2), 141–150 (1998)

76. Oleinik, E.F.: Glassy polymers as matrices for advanced composites. Polymer J. 19(1), 105–117 (1987)

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional

affiliations.