HACON - A Program For Temperature and Stress Simulation
HACON - A Program For Temperature and Stress Simulation
HACON - A Program For Temperature and Stress Simulation
Structural
Mechanics
Detta r en tom sida!
Structural Mechanics
HACON
A PROGRAM FOR SIMULATION
OF TEMPERATURE AND STRESS
IN HARDENING CONCRETE
i
ii
Contents
1 INTRODUCTION 1
1.1 General remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Characteristics of computer program . . . . . . . . . . . . . . . . . . 1
1.3 Summary of the report contents . . . . . . . . . . . . . . . . . . . . . 2
1.4 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 PROPERTIES OF HARDENING CONCRETE 5
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Hydration of concrete . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Compressive strength . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Elastic strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Thermal strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.6 Stress induced thermal strain . . . . . . . . . . . . . . . . . . . . . . 10
2.7 Autogeneous shrinkage . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.8 Creep strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.9 Fracturing strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.10 Stress-strain relation . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.11 Matrix formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3 DESCRIPTION OF THEORY FOR COMPUTATION OF TEM-
PERATURE 25
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Finite element formulation . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3 Innite element formulation . . . . . . . . . . . . . . . . . . . . . . . 30
3.4 Numerical solution procedure . . . . . . . . . . . . . . . . . . . . . . 32
3.5 Cooling pipes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4 DESCRIPTION OF THEORY FOR COMPUTATION OF DIS-
PLACEMENTS AND STRESSES 35
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Finite element formulation . . . . . . . . . . . . . . . . . . . . . . . . 35
4.3 Numerical solution procedure . . . . . . . . . . . . . . . . . . . . . . 39
4.4 Denition of equivalent length . . . . . . . . . . . . . . . . . . . . . . 41
iii
5 COMPUTER PROGRAM 43
5.1 General remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.2 Generation of input data . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.3 Presentation of output data . . . . . . . . . . . . . . . . . . . . . . . 46
A NOTATIONS 49
B DETERMINATION OF HEAT DEVELOPMENT PARAMETERS 53
B.1 Heat development in curing box . . . . . . . . . . . . . . . . . . . . . 53
B.2 Determination of maturity . . . . . . . . . . . . . . . . . . . . . . . . 55
B.3 Determination of degree of hydration . . . . . . . . . . . . . . . . . . 55
B.4 Computer code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
iv
Chapter 1
INTRODUCTION
1
Heat exchange with an unbounded region of e.g. rock is considered, using
nite elements.
The nite element model can be extended during the analysis to consider
casting stages performed at dierent points of time.
Variations in the environmental temperature and the removal of formwork and
insulation can be considered.
The temperature or generated heat can be prescribed at specic points to take
into consideration the eect of cooling pipes or heating cables.
In the stress analysis the development of material properties during hardening
is considered.
The constitutive equations take into account creep, stress induced thermal
strain, autogeneous shrinkage and crack developmet.
Input data to the program is dened interactively and mesh generation is in-
cluded.
Output data from the program is stored in a le, which is interpreted by a spe-
cial purpose post-processor. Distribution of temperature, degree of hydration and
equivalent maturity time are presented by colour or iso-lines. The history of the
temperature at a specied point and the mean, maximum and minimum temper-
ature, in a specied region are shown in diagrams. Displacements are illustrated
by showing a deformed nite element mesh. Magnitude and direction of principal
stresses are shown by arrows. Distribution of maximum principal stress, strength
and stress-strength ratio are presented by colour. The history of stress, strength
and stress-strength ratio is shown in diagrams.
1.4 Notations
The notations are explained in the text where they rst appear. They are also listed
in Appendix A.
References to literature are quoted in the text by numbers in square brackets,
[ ]. The references are given in alphabetical order in Appendix C.
Equations are, in general, written in component form and unless otherwise indi-
cated, the summation convention is applied, see e.g. Malvern [24].
2
A component of a vector is indicated by using indices, e.g. qi. To make a
distinction between the vector itself and its components, vectors are denoted by
boldface letters, e.g. q. The time derivative is denoted by a dot over the variable
considered, e.g. T_ .
3
4
Chapter 2
PROPERTIES OF HARDENING
CONCRETE
2.1 Introduction
Theoretical simulation of displacements and stresses in hardening concrete structures
require a proper material description. The material model, which will be described
below, aims to be general enough to re
ect the phenomena of interest, without
being more complicated than necessary. In
uence of hardening on the material
properties is taken into account. In the present work, the concrete is assumed not
to be aected by drying. Therefore, in
uence of drying on the material properties
is not considered. Since the tensile strength of concrete is much lower than the
compressive strength, a common situation is that the tensile strength is exceeded,
while the maximum compressive stress is well below the compressive strength. For
convenience, compression failure is therefore not included in the present model.
Heat development and development of material properties are dependent on the
hydration process. Therefore a description of this process is given below.
5
introduce the quantity equivalent maturity time te, which facilitates the comparison
of hydration processes during dierent thermal conditions. The equivalent maturity
time te is dened by the integral
Zt
f (T )
te = d (2.2)
f (Tr )
0
where f (T ) is a maturity function and Tr is a reference temperature, chosen as
Tr = 293K (20 C ). Several maturity functions have been proposed, see Byfors [5]
for a review. Freiesleben Hansen and Pedersen [14] proposed a maturity function
based on the Arrhenius equation for thermal activation, i.e.
f (T ) = ke T (2.3)
0
Use of Eq. (2.4) with the assumption of temperature dependency of , has
shown a good correlation between the equivalent maturity time and the compressive
concrete strength, see Freiesleben Hansen and Pedersen [14] and Byfors [5]. In the
present work Eq. (2.4) is used, with the quantity given by
0
T T
= 0 r a (2.5)
T Ta
as proposed by Jonasson [20]. In Eq. (2.5), 0 and 0 are material parameters
obtained experimentally, and Ta is a constant chosen as Ta = 263K ( 10C ).
The quantity of generated heat per mass unit of cement Wc is related to the
degree of hydration by Eq. (2.1). To obtain an expression for how the generated
heat is related to the equivalent maturity time te, a relation between and te is
required. In the present work the relation proposed by Jonasson [20] is adopted, i.e.
= e (ln(1+ t ))
te
1
1
1
(2.6)
where 1 , t1 and 1 are material parameters obtained experimentally.
In Appendix B a description of a method for determination of the material
parameters in Eqs. (2.5) and (2.6) from measurements using an insulated curing
box is given.
In a concrete structure the generated heat per volume unit is of interest. There-
fore, the generated heat per mass unit of cement is multiplied by the cement content
C , i.e. the generated heat per volume unit of concrete W is given by
W = CWc (2.7)
6
Substitution of Eq. (2.1) into Eq. (2.7) yields
W = CWc0 (2.8)
In the nite element equations, the time derivative of W will be used. Dieren-
tiation of Eq. (2.8) yields
dW d dte
Q= = CWc0 (2.9)
dt dte dt
where, according to Eqs. (2.6) and (2.4)
d
1 1
dte
= t1+1t ln 1 + tte (2.10)
1 e 1
e( Tr T )
dte
dt
= 1 1
(2.11)
7
2.4 Elastic strain
During hardening, the elasticity properties change. It is, however, reasonable to
assume that the deformation of concrete which is subjected to load is maintained
during hardening. This is because chemical reactions cause development of the inter-
nal material structure, and the new connection produced can be expected to be free
from stress and not carry the previously applied load. However, these connections
yield an increased material stiness. Unloading of the material will therefore show
an irreversible strain caused by hardening of loaded concrete. Assuming isotropy,
the elastic strain rate "_e is thus related to the stress rate _ by the generalized Hooke's
law
"_eij = Cijkm
e _
km (2.15)
where C e is the isotropic compliance tensor, given by
e
Cijkm = eij km + e (ik jm + im jk ) (2.16)
Here ij is Kronecker's delta and e and ke are elastic parameters, related to the
elastic modulus E and Poisson's ratio by
e = (2.17)
E
e =
1+ (2.18)
2E
The growth of the elastic modulus E during hardening may be described in the
same manner as the development of the compressive strength fc has been described
by Byfors [5], i.e.
E = E E0 (2.19)
where b E
a1E te 1
tr
= (b b ) (2.20)
1+ a1E te 1E 2E
a2E tr
and E0 is the elastic modulus at time te = tr . The parameters a1E , b1E , a2E and
b2E depend on the cement type and the concrete composition. Eq. (2.19), with
parameters chosen as a1E = 10(8:0 2:0W =C ) , b1E = 4:0; a2E = 1:0, and b2E = 0:1,
0
is compared with experimental data by Byfors [5] in Fig. 2.1. In accordance with
the experimental data, the value of E0 has been chosen as 38:0GP a, 38:0GP a and
31:0GP a for the water-cement ratios 0:40, 0:58, and 1:00, respectively. To get an
estimation of the value of E0 in Eq. (2.19) one may use the relation
s
fc0
E0 = E1 (2.21)
fr
where E1 = 6:0GP a and fr is a reference value, chosen as fr = 1:0MP a.
8
Figure 2.1: Comparison between Eq. (2.19) and experimental data according to
Byfors [5].
According to experimental data by Byfors [5] the value of Poisson's ratio de-
creases rapidly at a very early age, and then increases. This behaviour may be
described by the relation
te
= 1 e ( tr ) + 2 1 e ( tr ) (2.22)
1
te 1
where 1 and 2 are the initial and nal values of , respectively, and 1 and 2
(1 > 2 ) are parameters which express the in
uence of hardening. In Fig. 2.2
Eq. (2.22), with 1 = 0:5, 2 = 0:25, 1 = 150, and 2 = 40, is illustrated and
compared with experimental data according to Byfors [5].
9
Figure 2.2: Comparison between Eq. (2.22) and experimental data according to
Byfors [5].
2.6 Stress induced thermal strain
When the temperature of a concrete specimen exposed to external load is raised,
the observed deformation diers from the sum of the deformation of a loaded speci-
men at constant temperature and the thermal expansion of a non-loaded specimen.
The excess deformation is often called transitional thermal creep, and is sometimes
referred to as stress induced thermal stain. For temperatures below 1000C, the phe-
nomenon has been observed by e.g. Hansen and Eriksson [15], Illston and Sanders
[18], [19] and Parrott [28]. Khoury et al. [21] have made an excellent review of
experimental observations in this area. Important characteristics are that the strain
develops rapidly after increase of temperature to a level not previously attained and
that it is irrecoverable. The magnitude of the strain is proportional to the stress and
to the temperature change and may be regarded as unaected by maturity. The rate
of the stress induced strain "_T , as proposed by Thelandersson [34], can be assumed
to be given by the relation
T T_
"_Tij = T Cijkm (2.24)
km
where C T is a tensor which, assuming isotropy, can be expressed as
T = + ( + )
Cijkm (2.25)
T ij km T ik jm im jk
In Eq. (2.24) is the stress, T_ is the time derivative of temperature, and T is a
parameter with the value 1 when the temperature increases to a level not previously
attained, and 0 in other cases. In Eq. (2.25), T and T are material parameters,
which in analogy with Eqs. (2.19) and (2.20) can be expressed as
T = T (2.26)
T E
10
T =
1 + T (2.27)
2ET
According to the data given in Refs. [18], [28] and [21], ET varies between 0:5
and 4:0T P aK and a typical value is 2:0T P aK . The value of T cannot be obtained
from the data in these references, but data for concrete heated above 100C , given by
Thelandersson [34], yields T = 0:29, i.e. about the same as the value of Poisson's
ratio for hardened concrete. In absence of detailed experimental data, T may
be assumed to have a value in the range 0:25 0:30.
in which E is the elastic modulus and 0 and t are parameters, whose values can be
obtained from diagrams given in Ref. [4]. The relations can, with good agreement,
be described by the expressions
w
0 = 01 0 02 (2.31)
C
11
5
x 10
1
0
autogeneous shrinkage
5
0 200 400 600 800 1000
maturity (h)
where Eq. (2.32) is known as a Dirichlet series. The parameters may be chosen
as 01 = 3:2, 02 = 0:3; N = 5, 1 = 104s, 2 = 105s, 3 = 106s, 4 = 107s,
5 = 108 s, 1 = 0:06, 2 = 0:03, 3 = 0:20; 4 = 0:52 and 5 = 0:15. The parameter
t considers the in
uence of the concrete age at loading. This in
uence may be
described by the relation
0t1
te
0
t = tr
(2.33)
This expression, with 0t1 = 0:27, is illustrated in Fig. 2.4. In the gure data from
several experiments as compiled by Parrot [29] is also shown (from Byfors [5]). A
comparison with the curve given in Ref. [4] is also provided.
Adopting the principle of superposition, which means that the creep strain due
to a varying stress is equal to the sum of the creep strains due to stress increments
applied at dierent times t0, see e.g. Bazant [3], the creep strain may be written
Zt
d
"c = C 0 dt0 (2.34)
dt
0
Assuming isotropy, the following generalization of Eq. (2.34), see e.g. Bazant
12
Figure 2.4: Comparison between theory and experimental data.
[3], can be used
N Zt h
X t t0 i dkm 0
"cij = cn
Cijkm 1 e n
dt0
dt (2.35)
n=1 0
where the creep compliance tensor Ccncan be expressed in a way similar to the
elastic compliance tensor in Eq. (2.18), i.e.
cn = + ( + )
Cijkm (2.36)
cn ij km cn ik jm im jk
in which cn and cn are creep compliance parameters which, analogously to Eqs.
(2.21) and (2.22) can be expressed as
cn
cn = (2.37)
Ecn
cn =
1 + cn (2.38)
2Ecn
where Ecn is given by the relations
Ecn = 0 E0t n ; n = 1
Ecn = 0E0t0n ; n 2 (2.39)
and cn are quantities analogous to Poisson's ratio . To avoid great in
uence at a
late stage, from stress variations at an early stage, Ecn has for n 2 been assumed
13
to be proportional to E0 instead of E . In absence of detailed experimental results
one may assume that cn = , see e.g. Anderson [1].
The time derivative of the creep strain "cij , according to Eq. (2.35) is given by
"_ =
c
XN
1 e tn
n (2.40)
ij ij
n=1 n
where the quantities
ij express the history of the material and are given by
Zt
t0 dkm 0
ijn = cn e n
Cijkm
dt0
dt (2.41)
0
It may be observed that when the creep strain rate at time t + t is to be
computed, the integral of Eq. (2.41) has to be evaluated only from time t to time
t + t, if the value of the integral at time t is known. This means that the stress
history before time t need not be available when the creep strain rate at time t +t
is to be computed.
f r
14
Figure 2.5: Comparison between Eq. (2.42) and experimental data according to
Byfors [5].
where fr is a reference value chosen as fr = 1:0MP a. Eq. (2.42) is illustrated in
Fig. 2.5, with the parameters chosen as a1t = 10 (6: 0 2 : 00W
)
C , b1t = 3:0, a2t = 1:0,
b2t = 0:14, and ft1 = 0:3MP a. Experimental data representing the tensile strength
obtained by Byfors [5] in splitting tests is also shown in the gure.
The crack plane of the rst crack is assumed to be normal to the direction
of the maximum principal stress. A possible second crack is assumed to develop
perpendicular to the rst crack, if the normal stress in that direction also reaches
the tensile strength. Likewise a third crack may develop perpendicular to the existing
two cracks.
For convenience, a local coordinate system is introduced, where the x1-axis is
normal to the plane of the rst crack. If a second crack has occurred, the x2-axis is
normal to the plane of that crack. The unit vectors ii in the local coordinate system
are related to the unit vectors ik in the global coordinate system by the relation
ii = aik ik (2.45)
where
aik = cos (xi , xk ) (2.46)
in which xi and xk denote local and global coordinates respectively. It may be
noted that, in the present description of fracture, barred variables relate to the local
coordinate system.
The stress normal to the plane of crack No. is assumed to decrease when
the crack width w increases. It should be noted that when or is used as tensor
index, in this description, a specic component is assumed and that the summation
convention is not applied for a repeated index or .
15
Figure 2.6: Comparison between Eq. (2.48) and experimental data by Petersson
[31]
The energy GF necessary to produce one unit area of crack is given by the integral
Zwc
GF = dw (2.47)
0
where wc is the crack width when the normal stress has dropped to zero. According
to experimental results by Petersson [31], the fracture energy GF of hardening con-
crete develops in the same way as the elastic modulus E. Thus, the fracture energy
GF may be described as
GF = E GF 0 (2.48)
where GF 0 is the fracture energy at time te = tr . To get an estimation of the value
of GF 0 one may use the relation
GF 0 = GF 1 E0 (2.49)
where GF 1 = 2:5 mm gives good agreement with experimental data according to
Petersson [31], see Fig. 2.6.
The stress is for the sake of simplicity assumed to be a linear function of w,
i.e.
= ft + N w (2.50)
16
Figure 2.7: Relation of stress and crack width w
where ft is the uniaxial tensile strength and N is a proportionality factor (N < 0).
The relation between stress and crack width w is illustrated in Fig. 2.7.
Eq. (2.50) and the fact that the normal stress is zero when the crack width is
wc yield
f
wc = t (2.51)
N
Substitution of Eq. (2.50) into Eq. (2.47), evaluation of the integral and use of
Eq. (2.51) results in
ft2
N= (2.52)
2GF
A fracturing strain tensor "f is introduced which represents the mean fracturing
strain in a region which includes the discrete crack. The fracturing strain "f normal
to the plane of crack No. is dened by
w
"f =
(2.53)
L
where L is an equivalent length associated with crack No. which is related to
the size of the nite element where the crack develops.
Substitution of Eq. (2.53) into Eq. (2.50) and dierentiation with respect to
time yields
"_f = J _ (2.54)
where 1
J = (2.55)
NL
Eq. (2.55) is valid on condition that
w_ > 0 (2.56)
17
If this condition is not satised, i.e. unloading is taking place, the crack that has
developed is assumed to be irrecoverable and the fracturing strain rate is assumed
to be zero. This means that Eq. (2.54) is still valid, but now with J given by
J = 0 (2.57)
When a crack is fully developed, i.e. when w > wc, the stress is assumed to
be zero. If Eq. (2.54) is applied in this situation the value of J is innitely large.
As mentioned above, the fracture energy is given by Eq. (2.47), which means that
the remaining fracture energy of a developing crack is changed by the rate w_
when the rate of change in the crack width is w_ . In hardening concrete it may be
assumed that new connections develop in undamaged regions of fracturing concrete.
The remaining fracture energy G is therefore assumed to increase by a fraction
G =GF of the increase of GF , as described by Eq. (2.49). These assumptions lead
to the following description of the rate of change of the remaining fracture energy
G_ :
G_ = GGF G_ F w_ ; w_ > 0 (2.58)
If a crack, fully or partly developed, is exposed to compressive stress, the hardening
process may yield new connections even in damaged regions, since contact can be
expected between particles on dierent sides of the crack plane. The remaining
fracture energy G may therefore be assumed to increase by the same rate as the
fracture energy GF , i.e.
G_ = G_ F ; 0; w_ = 0 (2.59)
For tensile stresses below the current tensile strength f linear interpolation
between Eq. (2.58) and Eq. (2.59) may be applied, i.e. the rate of change of
fracture energy is given by
_G = 1 + f GGF 1 G_ F ; 0 f ; w_ = 0 (2.60)
where the tensile strength f is proportional to the square root of the remaining
portion of the fracture energy, i.e.
r
G
f = f (2.61)
GF t
Development of a crack will reduce the ability to transfer shear stresses across
the crack plane. In analysis of fracture it is important, in addition to modelling the
behaviour normal to a crack plane, to obtain a reasonable expression for the shear
behaviour. A common way of handling the reduction of shear stiness is simply
to multiply the elastic shear modulus by a so-called shear retention factor, which
has a constant value in the interval from 0 to 1. It is, however, not satisfactory to
assume that the reduction of the shear stiness is independent of the crack width.
18
It is reasonable to assume that the shear stiness is reduced gradually as the crack
width increases. Since the fracture is localized to a thin zone, the reduction of shear
stiness is related to this zone and the shear stiness in the region outside the zone
is unaected by cracking. When a smeared crack approach is used, an eective shear
modulus, representative of a region which includes the crack, is applied. The eective
shear modulus depends on the shear stiness of the fracture zone and also on the
shear stiness of the unaected region outside the zone, and is therefore dependent
on the size of the region considered. It is thus not acceptable to use a shear retention
factor which is not related to the size of the region of which it is assumed to be
representative. In experiments by Paulay and Loeber [30] relationships between
shear stress and shear displacement at constant crack widths have been obtained.
According to their results the shear displacement is more or less proportional to
the shear stress and to the crack width. Based on these observations the shear
displacement ws in crack No. in the direction x is by Ottosen and Dahlblom
[27] assumed to be given by
w
w s = ( 6= )
(2.62)
Gs
where w is the crack width, Gs is a material parameter, so-called slip modulus, and
is the shear stress in the crack plane. To obtain a rate formulation Eq. (2.62)
is dierentiated with respect to time. This yields an expression where the rate of
shear displacement is dependent on the rate of the normal stress in addition to the
rate of shear stress. This yields, however, a nonsymmetric system of equations.
Since it is advantageous to have a symmetric system of equations, the rate of shear
displacement and the rate of normal stress may be assumed to be uncoupled. This
is obtained by replacing Eq. (2.62) by the relation
w_
_
s = w (2.63)
G s
The value of Gs may be expected to increase during development of the elastic
modulus E . Thus, the slip modulus may be described as
Gs = E Gs0 (2.64)
where Gs0 is the slip modulus at time te = tr : To get an estimation of Gs0 one may
use the expression
Gs0 = Gs1 E0 (2.65)
where Gs1 = 1:0 10 4 may be chosen on the basis of Paulay and Loeber's data [30].
As in the case of fracture normal to the crack plane, a fracturing strain component
is dened for shear displacement. The fracturing shear strain " is dened by
f
w w
s s
" =
1
2 L + L ( 6= ) (2.66)
f
19
where L is the equivalent length corresponding to crack No. as introduced above
in Eq. (2.53). Dierentiation of Eq. (2.66) with respect to time and substitution of
Eqs. (2.53) and (2.63) into the result, and considering that _ = _ , yields
"_f = A _ ( 6= ) (2.67)
where
"f + "f
A = ( 6= ) (2.68)
2Gs
Using Eqs. (2.55) and (2.67) a complete relation between fracturing strain rate
and stress rate can be obtained
"_fpq = Cpqrs
f _ rs (2.69)
where
C1111
f
= J1 (2.70)
C2222
f
= J2 (2.71)
C3333
f
= J3 (2.72)
C1212
f
= C1221
f
= C2112
f
= C2121
f
= A12 2 (2.73)
C2323
f
= C2332
f
= C3223
f
= C3232
f
= A23 2 (2.74)
C3131
f
= C3113
f
= C1331
f
= C1313
f
= A31 2 (2.75)
All other components Cpqrs
f are zero.
The stress and fracturing strain rate tensors expressed in the local coordinate
system can be transformed to the global coordinate system by the usual relations.
_ rs = ark asm _ km (2.76)
"_fij = api aqj "_fpq (2.77)
where aik is given by Eq. (2.46), Eqs. (2.69), (2.76) and (2.77) can be combined to
form a relation between fracturing strain rate and stress rate in the global coordinate
system
"_fij = apiaqj Cpqrs
f a a _
rk sm km (2.78)
20
2.10 Stress-strain relation
The strain rate tensor "_ is assumed to consist of the sum of the strain rates as
described in the previous sections, i.e.
"_ij = "_eij + "_Tij + "_Tij + "_aij + "_cij + "_fij (2.79)
To obtain a relation between the stress rate _ and the total strain rate "_ the
previous expressions will be combined. Substitution of Eqs. (2.15) and (2.78) into
Eq. (2.79) yields
"_ij = Cijkm
e + apiaqj Cpqrs
f a a _ 0
rk sm km + "_ij (2.80)
where
"_0ij = "_Tij + "_Tij + "_aij + "_cij (2.81)
Since the tensor C e is isotropic, its components do not change under rotation of
the coordinate axes, Eq. (2.80) can be written as
"_ij = api aqj Cpqrsark asm _ km + "_0ij (2.82)
where
Cpqrs = Cpqrs
e +C pqrs
f (2.83)
Eq. (2.82) can be used to determine the strain rate when the stress rate is given.
To obtain a nite element formulation, the inverse formulation is of interest.
Therefore, the stiness tensor D is introduced and dened as the inverse of the
compliance tensor C , i.e.
1
D tupq Cpqrs = (tr us + ts ur ) (2.84)
2
Multiplication of Eq. (2.82) by avn awoD vwtu ati auj , use of the relations aik ajk =
ij , aik aim = km , _ km = _ mk , and Eq. (2.84), and rearrangement yields the relation
_ km = Dkmij "_ij 0
km (2.85)
where the material stiness D and the pseudo stress rate _ 0 are given by
Dkmij = ark asm D rspq api aqj (2.86)
0 = Dkmij "_0
_ km (2.87)
ij
21
2.11 Matrix formulation
For computer programming, matrix formulation is in general most convenient. The
proposed constitutive equations are therefore summarized below, using matrix no-
tation.
Since the stress and strain tensors are symmetric, the previous equations can be
expressed by using six components instead of nine. The constitutive equations, in
the present work, are applied to two-dimensional problems, where two of the shear
components are zero, and therefore may be excluded. The tensors _ , "_ and "_0 are
therefore expressed by four-element column matrices, i.e.
T
_ = _ 11 _ 22 _ 12 _ 33 (2.88)
T
"_ = "_11 "_22 2"_12 "_33 (2.89)
T
"_0 = "_011 "_022 2"_012 "_033 (2.90)
It should be noted that the shear strain components have been multiplied by 2.
Using matrix notation, Eq. (2.82) can be written
"_ = G C G T _ + "_0 (2.91)
where 2 3
a11 a11 a21 a21 a11 a21 a31 a31
G = a12 a12 a22 a22 a12 a22 a32 a32
6 7
6
42a11 a12 2a21a22 a11 a22 + a21 a12 2a31 a32
7
5 (2.92)
a13 a13 a23 a23 a13 a23 a33 a33
2 1 3
E+ J1
E 0 E
E + J2 0
1 7
C = 6
6
4
E
0 0 2(1+E ) + 2A12 1 0 5
E 7 (2.93)
E E 0 E + J3
"_0 = "_T + "_T + "_c (2.94)
"_T = T (rT + T (1 rT ))[1 1 0 1]T (2.95)
2 T 3 2
1
ET
T
ET 0 ET 11
3
6 T
"_T = T T 6
6 ET E
1 0 T 7 6
E T 7 22 7 (2.96)
T
2(1+ T ) 76 7
4 0 0 ET 0 5 4 12 5
T
ET
T
ET 0 1
ET
33
"_c =
N
X 1e t
n
n (2.97)
n=1 n
22
2 32
Zt
1
Ecn
cn
Ecn 0 cn
Ecn
d11
dt0
3
n =
t0
6
e n 6
cn
Ecn
1
Ecn 0 cn 7 6 d
Ecn 7
76
22
dt0
7 0
7 dt (2.98)
4
6
0 0 2(1+cn ) 0 d
5 4 dt12
0 5
Ecn
0 cn cn 0 1 d33
Ecn Ecn Ecn dt0
Dening G and D as
G = G 1 (2.99)
D = C 1 (2.100)
and premultiplying Eq. (2.91) by GT DG yields
_ = D"_ _ 0 (2.101)
where the material stiness D and pseudo stress rate _ 0 are given by
D = GT DG (2.102)
_ 0 = D"_0 (2.103)
For axial symmetry, the matrices given by Eqs. (2.102) and (2.103) can be
directly used to establish the element matrices. For plane strain and plane stress,
the matrices have to be reduced to three components, using the conditions "_33 = 0
and _ 33 = 0, respectively, before establishing element matrices.
23
24
Chapter 3
DESCRIPTION OF THEORY
FOR COMPUTATION OF
TEMPERATURE
3.1 Introduction
To facilitate simulation of thermal problems numerical methods are used. The nite
element method has proved to be an ecient tool in thermal analysis and is therefore
applied. In this chapter a nite element formulation of the heat conduction problem
is given. An innite element formulation is also given. This concept is applied to
take into account unbounded regions of e.g. rock.
It should be noted that the index m may normally assume values exceeding 3.
25
Use of the divergence theorem gives the relation
Z Z Z Z
@vm _
qi dV + vm cT dV = vm qs dS + vm QdV (3.3)
@xi
V V S V
26
Figure 3.1: Geometry and degrees of freedom of eight-node isoparametric element
27
The present work deals with plane and axisymmetric problems, using an eight-
node isoparametric element, see Fig. 3.1. The element is obtained by distorting the
rectangular parent element of Fig. 3.2.
The geometry of the element is given by the relations
x1 = Nn xn (3.13)
and
x2 = Nn yn (3.14)
where the interpolation functions, Nn may be found in e.g. Zienkiewicz and Taylor
[37] and are given by
N1 = 41 (1 )(1 )( 1) ; N2 = 21 (1 2) (1 ) ;
N3 = 4 (1 + )(1 ) ( 1) ;
1 N4 = 21 (1 + )(1 2 ); (3.15)
N5 = 14 (1 + )(1 + )( + 1); N6 = 21 (1 2) (1 + );
N7 = 14 (1 )(1 + ) ( + 1); N8 = 21 (1 )(1 2 )
The interpolation functions Nn of Eq. (3.15) are used also for the temperature
in Eq. (3.7).
According to the chain rule, the derivatives of Nn with respect to and can
be expressed in the derivatives with respect to x1 and x2 as
@Nn @Nn
@
@Nn = @x1
J @N n (3.16)
@ @x2
The derivatives of Nn with respect to x1 and x2 can be obtained from Eq. (3.16),
i.e. @Nn @Nn
@x1
@Nn =J 1 @
@Nn (3.19)
@x2 @
where J 1 is the inverse of the Jacobian matrix J of Eq. (3.18). The derivatives of
Nn with respect to and are obtained by dierentiation of the expressions of Eq.
28
(3.15), i.e.
@N1
@ = 14 (1 ) (2 + ); @N = (1 );
@
2
@N3
@ = 14 (1 ) (2 ); @N = 1 (1 2 ) ;
@ 2
4
@N5
@ = 14 (1 + )(2 + ) ; @N = (1 + );
@
6
@N7
@ = 14 (1 + )(2 ); @N = 1 (1 2 );
@ 2
8
(3.20)
@N1
@ = 14 (1 )(2 + ); @N = 1 (1 2 ) ;
@ 2
2
@N3
@ = 14 (1 + )(2 ); @N = (1 + );
@
4
@N5
@ = 14 (1 + )(2 + ); @N = 1 (1 2 ) ;
@ 2
6
@N7
@ = 14 (1 )(2 ); @N = (1 )
@
8
The volume integrals of Eqs. (3.10)-(3.12) are transformed into integrals over
the area of the parent element by the substitution
dV = t det(J )dd (3.21)
for plane problems and
dV = 2x1 det(J )dd (3.22)
for axisymmetric problems, with x1 and x2 as radial and axial coordinate, respec-
tively. The integrals can be evaluated numerically by Gauss quadrature as
Z1 Z1 n X
X n
f (; ) dd = Hi Hj f (i ; j ) (3.23)
1 1 j =1 i=1
where, for n = 2,
H1 = H2 = 1 (3.24)
1 = 1 = p1 (3.25)
3
1
2 = 2 = p (3.26)
3
For other values of n the Hi; i and i values may be found in e.g. Zienkiewicz
and Taylor [37].
The surface integral of Eq. (3.12) is transformed into an integral over the side
of a parent element by the substitution
s
dx 2
2
dy
dS = t + d (3.27)
d d
for plane problems and
s
dx 2
2
dy
dS = 2x1 + d (3.28)
d d
29
for axisymmetric problems.
The integral can be evaluated numerically by Gauss quadrature as
Z1 n
X
f ( )d = Hi f (i ) (3.29)
1 i=1
A prescribed heat
ow into the body studied is represented by the quantity qs
in Eq. (3.12). To consider the in
uence of a prescribed environmental temperature,
the heat
ow qs is assumed to be given by the relation
qs = t (T0 T ) (3.30)
where t is the transfer coecient, T0 is the environmental temperature and T is
the surface temperature.
Substitution of Eq. (3.7) into Eq. (3.30) yields
qs = t T0 t Nn Tn (3.31)
Since the second term of Eq. (3.31) is a function of the nodal temperatures Tn,
the integral of this term should be moved to the left-hand side of Eq. (3.9).
N3 = (1+)(1 1) ; N4 = (1+2(1)(1+ )
) ;
(3.32)
(1
N8 = 2(1 ) )(1+)
30
Figure 3.3: Geometry and degrees of freedom of ve-node singly innite element
Figure 3.4: Geometry and degrees of freedom of three-node doubly innite element
31
Dierentation of the expressions in Eq. (3.20) yields
@N = 2 + ;
1 @N =2 4 ;
@ 12 @ 1
@N =
3
; @N = 1+ ;
4
@ 1 @ 2(1 )
@N = 1+
2(1 ) ;
8
@N = 2+ + ;
1
2
@
@N = 2(1 ) ;
2
2 (3.33)
@ (1 )
2 @ (1 )2
@N = 2 + ;
3
2 @N = 1+ ;
4
@ (1 )
2 @ (1 ) 2
@N = 1
8
@ (1 ) 2
which should be used in Eq. (3.18) instead of the expressions in Eq. (3.20).
The doubly innite element has three nite nodes, as shown in Fig. 3.4.
The mapping functions to be used instead of Eq. (3.15) in Eqs. (3.13) and (3.14)
are, according to Marques and Owen [25], given by
N1 = (1+3()(1 )1) ; N2 = (1 2(1+ )
)(1 ) ;
N8 = (1 2(1+ ) (3.34)
)(1 )
Dierentation of Eqs. (3.34) yields
@N = 2(3+)
(1 ) (1 ) ;
@N =
(1 ) (1 ) ;
1 2 4
@ 2 @ 2
@N = 2(1+) ;
8
@N =
1 2(3+) ;
@ (1 ) (1 )
2
@N = 2(1+ ) ;
2
(3.35)
@ (1 )(1 )
2 @ (1 )(1 )
2
@N =8 4
@ (1 )(1 )
2
which should be used in Eq. (3.18) instead of the expressions in Eq. (3.20).
The nodes at innity are assumed to have a constant temperature. Therefore,
the terms Kmn Tn in Eq. (3.9) which correspond to nodes at innity are known and
can be moved to the right-hand side. Since the temperature at the innite nodes is
constant, the time derivatives are zero. This means that a set of equations, expressed
in the nite nodes only, is obtained.
32
T7 T6 T5 T7 T6 T5
T8 T4 T8 T4
1
2 (T2 + T8)
T1
T1
T2 T3 T2 T3
Figure 3.5: Geometry and degrees of freedom for element for pipe modelling
parameter is in the present work chosen as = 1. With this choice Eq. (3.31) can
be written as
K mn Tn (t) = Pm (3.37)
where 1
K mn = Kmn + Cmn (3.38)
t
1
Pm = Pm + Cmn Tn (t t) (3.39)
t
This procedure is known as a backward dierence method and is unconditionally
stable.
33
In addition to pipe diameter, the transfer coecient of the pipe wall and the
temperature inside the pipe are used in the computation. For determination of the
pipe temperature it is considered that a pipe arrangement may consist of several
pipe parts connected to each other. Each pipe part is assumed to have the length
lp perpendicular to the plane of computation. To each arrangement the input tem-
perature Tin and either output temperature Tout or pipe
ow qp has to be specied.
An approximative value of the heat Qp to be removed from the concrete by the pipe
arrangement can be obtained from the results of the previous time step as
np
X
Qp = qip lp (3.40)
ip =1
34
Chapter 4
DESCRIPTION OF THEORY
FOR COMPUTATION OF
DISPLACEMENTS AND
STRESSES
4.1 Introduction
In the previous chapter the theory for computation of temperatures was described.
The temperature distributions obtained are used as input for the computation of
displacements and stresses. The theory for this computation will be described in
this chapter. Just as for computation of temperatures, the nite element method is
used for computation of displacements and stresses. This chapter presents a nite
element formulation of the displacement and stress problem.
35
The index m may normally exceed 3. Using the divergence theorem yields
Z Z Z
@vim
dV
@xj ji
= vim ti dS + vim fi dV (4.3)
V S V
where t is the traction vector dened by
ti = ji nj (4.4)
where n is the outward unit vector normal to the surface S of the body studied.
Since vim is assumed to be constant with respect to time, Eq. (4.3) yields
Z Z Z
@vim
_ ji dV = vim t_i dS + vim f_i dV (4.5)
@xj
V S V
In Chapter 2 the relation between stress rate _ and strain rate "_ was expressed
as (cf. Eq. 2.85)
_ ij = Djipq "_pq _ ji0 (4.6)
where D is the material stiness and _ 0 is the rate of pseudo stress. The strain rate
"_ is assumed to be given by the kinematic relation
"_pq =
1 @up @uq
2 @xq + @xp (4.7)
where u is the displacement vector.
Eqs. (4.5), (4.6) and (4.7) are combined to form the relation
Z Z Z Z
@vim
Djipq
1 @ u_ p @ u_ q
+ dV = vim t_i dS + vim fi dV + _ @vim 0
_ dV (4.8)
V
@xj 2 @xq @xp
S V V
@xj ji
Since the stress rate and strain rate tensors _ and "_ are symmetric, the material
stiness tensor D has the symmetry properties
Dijpq = Djipq (4.9)
Dijpq = Dijqp (4.10)
Due to these symmetry properties, Eq. (4.8) can be written as
Z Z Z Z
@vim @ u_ p _ _ @vim 0
Dijpq dV = vim ti dS + vim fi dV + _ ji dV (4.11)
@xj @xq @xj
V S V V
To obtain nite element equations, the displacements u are expressed as a func-
tion of the nodal displacements un as
up = Npn un (4.12)
36
Figure 4.1: Geometry and degrees of freedom of eight-node isoparametric element
where Npm are interpolation functions. The summation index in Eq. (4.12) may
normally exceed 3. Since the interpolation functions Npn are assumed to be constant
with respect to time, dierentiation of Eq. (4.12) yields
u_ p = Npn u_ n (4.13)
The interpolation functions are chosen as weighting functions, according to the
Galerkin method, i.e.
vim = Nim (4.14)
Substitution of Eqs. (4.13) and (4.14) into Eq. (4.11) yields
Kmn u_ n = P_m + P_m0 (4.15)
where
@ Nim @ N
Z
Kmn = Dijpq pn dV (4.16)
@x j @x q
V
Z Z
P_m = Nim t_i dS + Nim f_i dV (4.17)
S V
@ Nim
Z
P_m0 = _ ij0 dV (4.18)
@xj
V
The present work deals with plane stress, plane strain and axisymmetric problems
and, just as for temperatures described in Chapter 3, an eight-node isoparametric
element is used, see Fig. 4.1.
The element is obtained by distorting the rectangular parent element shown in
Fig. 4.2. The geometry of the element is given by the relations
x1 = Nn xn (4.19)
37
Figure 4.2: Parent element
and
x2 = Nn yn (4.20)
where the interpolation functions Nn are given by Eq. (3.15).
The interpolation functions Npn for displacements, introduced in Eq. (4.12), are
given by the relations
N1 2n 1 = Nn ; N1 2n = 0; (4.21)
N2 2n 1 = 0; N2 2n = Nn
To obtain a formulation suitable for computer programming, the symmetry prop-
erties of the tensors D and _ 0, and the fact that some terms are zero for plane and
axisymmetric problems, can be utilized. Eqs. (4.16) - (4.18) can thus be expressed
as Z
Kmn = BmT DBndV (4.22)
V
Z Z
P_m = NmT tdS + NmT fdV (4.23)
S V
Z
P_ m0 = BmT _ 0 dV (4.24)
V
where D and _ 0 are dened by Eqs. (2.102) and (2.103) and the matrices t_, f_, and
NnT are given by
t_ = t_1 t_2 T
(4.25)
f = [f1 f2 ]T (4.26)
38
Nn = N1 n N2
T
n (4.27)
The matrix Bn is for plane stress and plane strain given by
h iT
Bn = @ N1 n @ N2 n @ N1 n + @@xN n 2
(4.28)
@x1 @x2 @x2 1
The derivatives of the shape functions are given by Eq. (3.19) and the integrals
in Eqs. (4.22) - (4.24) are evaluated by Gauss quadrature according to Eqs. (3.23)
and (3.29).
39
where the total external loads Pm are given by
Z Z
Pm = NmT tdS + NmT fdV (4.36)
S V
The last term in Eq. (4.35) represents the internal forces, corresponding to the
stress obtained , as computed according to the constitutive equations. Dierentia-
tion of Eq. (4.35) with respect to time yields the rate of change of the out-of-balances
forces Pm , i.e. Z
_Pm = P_m BmT _ dV (4.37)
V
Substitution of Eq. (2.101) and the relation
"_ = Bn u_ n (4.38)
into Eq. (4.37) yields the expression
Kmn u_ n = P_m + P_m0 P_m (4.39)
where Kmn and P_m0 are dened by Eqs. (4.22) and (4.23), respectively.
Evaluation of Eq. (4.39) can be performed by an Euler forward expression,
assuming a linear relation from time t to time t + t, i.e.
Kmn (t)un = Pm + tP_m0 (t) tP_m (4.40)
where Kmn (t) is the tangential stiness at time t, un the change of nodal displace-
ments from time t to time t +t, P_m0 (t) the rate of pseudo load at time t, and Pm
incremental load, given by
tZ+t
Pm = P_m dt (4.41)
t
The rate of change of the out-of-balance forces P_m is assumed to be given by
P_m =
1 P (t) (4.42)
t m
to make the out-of-balance forces reduce to zero, during the time increment, and
guide the solution towards the true one. Substitution of Eq. (4.42) into Eq. (4.40)
yeilds
Kmn (t)un = Pm + tP_m0 (t) + Pm (t) (4.43)
This procedure is described by Stricklin et al. [33], and is called a self-correcting
form, since the solution is guided towards the true one. Making P_m = 0 in Eq.
(4.40) yields a purely incremental stiness method, which yields solutions with a
tendency to drift away from the true solution, because of increasing out-of-balance
forces.
40
Figure 4.3: Illustration of equivalent length in eight-node element with four integra-
tion points.
4.4 Denition of equivalent length
In Eqs. (2.54) and (2.66), when the fracturing strain was dened, an equivalent
length L was introduced. This length is related to the size of the nite element
where the crack develops. In the present work the equivalent length is chosen ac-
cording to the denition proposed by Ottosen and Dahlblom [27]. The equivalent
length is thus dened as the maximum length of the element region of interest, in the
direction normal to the crack plane. For an eight-node element, the region bounded
by the nodal points closest to the integration point in question and the element
mid-point is considered. This denition is illustrated in Fig. 4.3.
41
42
Chapter 5
COMPUTER PROGRAM
43
Figure 5.1: Graphical user interface to HACON
The signicance of most of the data to be input is obvious. The manner of
specifying the geometry may, however, need an explanation. The nite element
mesh is established by a mesh generator which is based on the technique used by
Liu and Chen [7]. According to this technique, the region studied is divided into a
number of superelements, which are rened to nite elements by the mesh generation
algorithm. The superelement mesh is established by dividing the region studied into
four-sided sub-regions. A superelement where all the sides are straight is dened by
the four superelement nodes at the corners.
If any of the sides of a superelement are curved, the superelement is dened by
midside superelement nodes, in addition to the corner superelement nodes. The
superelement node numbers have to be given counter-clockwise, starting with a
corner node. In Fig. 5.2 an example of a subdivision of a region into superelements is
shown. In Tab. 5.1 the corresponding data is given. In addition to the superelement
node numbers, the table also shows the number of element rows and element columns
the superelements are to be subdivided into. The resulting mesh is shown in Fig. 5.3.
It should be noted that the numbers of rows and columns refer to a local coordinate
system in a superelement and that the corner node which is given rst is situated
at the lower-left corner according to this local system. If a region is subdivided
into rectangular superelements, it is therefore recommended that the superelement
44
Figure 5.2: Subdivision of region into superelements.
Super- Superelement nodes Number of Number of
element rows columns
1 1264 3 3
2 2376 3 4
3 4 5 6 9 12 11 10 8 4 3
4 6 7 13 12 4 4
5 16 14 10 11 12 15 18 17 3 4
Table 5.1: Superelement data for region in Fig. 5.2.
node positioned at the lower-left corner is given as the rst of the nodes dening the
superelement. It should also be observed that the subdivision of two neighbouring
superelements has to be the same at the common boundary.
The nite element solutions of the temperature and displacement elds will con-
verge towards the exact solutions, as the element size is decreased. As the time
step is decreased, the numerical solution procedure gives a solution which converges
towards the correct one. The element size and the time step suitable for a compu-
tation depends on the problem studied and the requirements of accuracy. To nd
out what is suitable to use in a specic situation, one may perform computations
with dierent element sizes, or dierent time steps, and study the in
uence on the
result.
When the object studied includes an unbounded region of rock or old concrete
which will be considered in the temperature modelling, superelement nodes are not
positioned at the corners at innity. Instead, superelement nodes are positioned
45
Figure 5.3: Finite element mesh resulting from the geometry in Fig. 5.2 and the
data in Tab. 5.1.
where nite nodes are to be localized. The region will extend to innity beyond
these superelement nodes.
46
presented in diagrams or tables.
The history of maximum principal stress, strength and stress-strength ratio at a
specied point or in a specied region, can be presented in diagrams or tables.
47
48
Appendix A
NOTATIONS
Notations are explained in the text where they rst occur. Most of them are also
given in this appendix.
a1c development parameter for compressive strength
a1E development parameter for elastic modulus
a1t development parameter for tensile strength
a2c development parameter for compressive strength
a2E development parameter for elastic modulus
a2t development parameter for tensile strength
b1c development parameter for compressive strength
b1E development parameter for elastic modulus
b1t development parameter for tensile strength
b2c development parameter for compressive strength
b2E development parameter for elastic modulus
b2t development parameter for tensile strength
C heat capacity matrix, cement content
c specic heat
D material stiness
E elastic modulus
E0 elastic modulus at te = tr
Ecn creep parameter
ET parameter for stress induced thermal strain
f body force
fc0 compressive strength at te = tr
ft tensile strength
ft0 tensile strength at te = tr
GF fracture energy
GF 0 fracture energy at te = tr
Gs slip modulus
Gs0 slip modulus at te = tr
J Jacobian matrix
49
K conductivity, stiness matrix
k thermal conductivity
N interpolation function
n unit vector
P heat
ow matrix, load matrix
P0 pseudo load matrix
Q generated heat
q heat
ow vector
qs heat
ow into body
rT recoverability of thermal strain
T temperature
T nodal temperature matrix
T0 environmental temperature
Tr reference temperature (293 K = 200C)
t time
t traction
t1 material parameter for cement
te equivalent maturity time
tr reference time (672 h = 28 d)
u displacement
u nodal displacement matrix
v weighting function
wC quantity of heat developed
Wc0 quantity of heat developed at complete hydration
w crack width
ws shear displacement in crack
xi coordinate
degree of hydration
1 development parameter for Poisson's ratio
2 development parameter for Poisson's ratio
T coecient of thermal expansion
t transfer coecient
ij Kronecker's delta
" strain
"a autogeneous shrinkage strain
"c creep strain
"e elastic strain
"f fracturing strain
"T thermal strain
"T stress induced thermal strain
E development function for elastic modulus
t development function for tensile strength
parameter used in time stepping procedure,
50
activation energy divided by universal gas constant
0 material parameter for cement
k0 material parameter for cement
k1 material parameter for cement
1 material parameter for cement
Poisson's ratio
1 initial value of Poisson's ratio
2 nal value of Poisson's ratio
cn creep parameter
T parameter for stress induced thermal strain
mass density
stress
_ 0 pseudo stress rate
n creep parameter
0 creep parameter
n creep parameter
t creep parameter
0t creep parameter
51
52
Appendix B
DETERMINATION OF HEAT
DEVELOPMENT
PARAMETERS
53
The data obtained from the experiment is the concrete temperature and the
environmental temperature at dierent times. During the time from ti 1 to ti the
mean value of the concrete temperature, the mean value of the environmental tem-
perature, and the time derivative of the concrete temperature are assumed to be
1
T = (T (ti 1 ) + T (ti )) (B.6)
2
1
Te = (Te (ti 1 ) + Te (ti )) (B.7)
2
T (t ) T (ti 1 )
T_ = i (B.8)
ti ti 1
respectively. A value of the coecient k may be obtained by substituting these
values into Eq. (B.5), i.e.
1
k(ti ) = T_ (B.9)
T Te
When k is to be determined the mean value of the values k(ti ) obtained at
dierent times may be computed, i.e.
k=
1X n
k(t ) (B.10)
i
n i=1
Assuming the value of k has been determined, the mean value of the heat devel-
opment Q (ti) during the time from ti 1 to ti may be obtained from Eqs. (B.1)-(B.3)
with the temperature and its time derivative given by Eqs. (B.6)-(B.8).
According to Eq.(2.9) the total heat developed W is given by
Zt
W = Qd (B.11)
0
Since the mean value of Q during each time interval is obtained as described
above, the integral in Eq. (B.11) may be evaluated as
n
X
W = Q (ti )(ti ti 1 ) (B.12)
i=1
The generated heat per mass unit of cement We is related to W by Eq. (2.7).
Thus,
W
Wc = (B.13)
C
where C is the cement content.
54
B.2 Determination of maturity
The equivalent maturity time is assumed to be given by Eqs. (2.4) and (2.5). Since
the temperature is obtained at specic time intervals the integral of Eq. (2.4) may
be evaluated as n
te = e[ Tr T ] (ti ti 1 )
X 1 1
(B.14)
i=1
with 0
T T
= 0 r a (B.15)
T Ta
where T is the mean value of the concrete temperature during the time interval,
according to Eq. (B.6).
If temperature data from experiments with a certain concrete mix at dierent
casting temperatures are available, the same maturity should be obtained by Eq.
(B.14) for the dierent data sets at each specic value of generated heat Wc. This
means that the same te - Wc -curve should be obtained for the dierent data sets.
The values of 0 and k0 may be optimized to minimize the dierence between the
curves.
using the method of least squares. The values of t1 may be optimized to minimize
the deviation between the experimental data and the theoretical expression given
by Eq. (2.6).
The in
uence of the values of Wc0 and t1 on the deviation between experimental
data and the theoretical expression is not very strong. The heat development of one
55
type of cement may therefore with good agreement be described by several dierent
combinations of Wc0, 1 , t1 and 1 .
56
Bibliography
[1] Anderson, C.A.: Numerical creep analysis of structures, Chapter 8 of Creep and
Shrinkage in Concrete Structures, edited by Z.P. Bazant and F.H. Wittman,
John Wiley & Sons, 1982.
[2] Auperin, M., de Larrard, F., Richard, P. and Acker, P.: (Shrinkage and creep
of high-strength concrete - In
uence of age at loading), Annales de l'Institut
Technique du Batiment et des Travaux Publics,No. 474, pp. 50-76, 1989.
[3] Bazant, Z.P.: Mathematical models for creep and shrinkage of concrete, Chap-
ter 7 of Creep and Shrinkage in Concrete Structures, edited by Z.P. Bazant
and F.H. Wittmann, John Wiley & Sons, 1982.
[4] Betonghandbok - Material (Swedish Handbook for Concrete Construction -
Material), Svensk Byggtjanst, Stockholm, 1994.
[5] Byfors, J.: Plain Concrete at Early Ages, Swedish Cement and Concrete Re-
search Institute, FO 3:80, Stockholm 1980.
[6] Carslaw, H.S. and Jaeger, J.C.: Conduction of Heat in Solids, Second edition,
Oxford University Press, Oxford 1959.
[7] Dahlblom, O.: Constitutive modelling and nite element analysis of concrete
structures with regard to environmental in
uence, Report TVSM-1004, Lund
Institute of Technology, Division of Structural Mechanics, Lund 1987.
[8] Dahlblom, O.: HACON-T - A program for simulation of temperature in harden-
ing concrete, R, D & D-report, Serial number U 1990/31, Vattenfall Utveckling
AB, 1990.
[9] Dahlblom, O.: HACON-S - A program for simulation of stress in hardening
concrete, Report from Vattenfall Hydro Power Genertion, Serial number H
1992/3, Vattenfall Utveckling AB, 1990.
[10] Dahlblom, O.: HACON-H - A program for determination of heat development
parameters for hardening concrete, Report from Vattenfall Hydro Power Gen-
ertion, Serial number H 1992/4, Vattenfall Utveckling AB, 1990.
57
[11] Dahlblom, O. and Ottosen, N.S.: Smeared crack analysis using generalized
ctitious crack model, Journal of Engineering Mechanics, Vol. 116, No. 1, pp.
55-76, 1990.
[12] Damjanic, F. and Owen D.R.J.: Mapped Innite Elements in Transient Ther-
mal Analysis, Computers and Structures, Vol. 19, No. 4. pp. 673-687, 1984.
[13] Emborg, M.: Thermal stresses in concrete structures at early ages, Report
1989:73D, Lulea University of Technology, Division of Structural Engineering,
Lulea 1989.
[14] Freiesleben Hansen, P. and Pedersen, E.J.: Maturity Computer for Controlled
Curing and Hardening of Concrete (in Danish), Journal of the Nordic Concrete
Federation, No. 1, pp. 21-25, 1977.
[15] Hansen, T.C. and Eriksson, L.: Temperature change eect on behaviour of
cement paste, mortar and concrete under load, Journal of ACI 63, pp. 489-504,
1966.
[16] Hillerborg, A.: A model for fracture analysis, Report TVBM-3005, Lund Insti-
tute of Technology, Division of Building Materials, Lund 1978.
[17] Hillerborg, A., Modeer, M. and Petersson, P-E.: Analysis of crack formation
and crack growth in concrete by means of fracture mechanics and nite ele-
ments, Cement and Concrete Research, Vol. 6, pp. 773-782, 1976.
[18] Illston, J.M. and Sanders, P.D.: Characteristics and prediction of creep of satu-
rated mortar under variable temperature, Magazine of Concrete Research, Vol.
26, No. 88, pp. 169-179, 1974.
[19] Illston, J.M. and Sanders, P.D.: The eect of temperature change upon the
creep of mortar under torsional loading, Magazine of Concrete Research, Vol.
25, No. 84, pp. 136-144, 1973.
[20] Jonasson, J-E.: Slipform Construction - Calculations for Assessing Protection
Against Early Freezing, Swedish Cement and Concrete Research Institute, FO
4:84, Stockholm 1985.
[21] Khoury, G.A., Grainger, B.N. and Sullivan, P.J.E.: Transient thermal strain
of concrete: literature review, conditions within specimen and behaviour of
individual constituents, Magazine of Concrete Research, Vol. 37, No. 132, pp.
131-144, 1985.
[22] Liu, Y. and Chen, K.: A two-dimensional mesh generator for variable order
triangular and rectangular elements, Computers and Structures, Vol. 29, No.
6, pp. 1033-1053, 1988.
58
[23] Lofquist, B.: Temperatureekter i hardnande betong (Temperature eects in
hardening concrete), Kungl. Vattenfallsstyrelsen, Teknisk meddelande, Serie B,
Nr. 22, Stockholm 1946.
[24] Malvern, L.E.: Introduction to the Mechanics of a Continous Medium, Prentice-
Hall, Englewood Clis, New Jersey 1969.
[25] Marques, J.M.M.C. and Owen, D.R.J.: Innite Elements in Quasi-Static Ma-
terially Nonlinear Problems, Computers and Structures, Vol. 18, No. 4, pp.
739-751, 1984.
[26] Modeer, M.: A fracture mechanics approach to failure analysis of concrete
materials, University of Lund, Division of Building Materials, Lund 1979.
[27] Ottosen, N.S. and Dahlblom, O.: Smeared crack analysis using a nonlinear
fracture model for concrete, Numerical Methods for Non-Linear Problems, Vol.
3, Edited by C. Taylor, D.R.J. Owen, E. Hinton and F.B. Damjanic, Pineridge
Press, pp. 363-376, 1986.
[28] Parrott, L.J.: A study of transitional thermal creep in hardened cement paste,
Magazine of Concrete Research, Vol. 31, No. 107, pp. 99-103, 1979.
[29] Parrott, L.: Eect of loading at early age upon creep and relaxation of concrete.
RILEM Committee 42-CEA, International report UK5, 1978.
[30] Paulay, T. and Loeber, P.S.: Shear transfer by aggregate interlock, Shear in
reinforced concrete, Volume 1, Special Publications SP-42, pp. 1-15, American
Concrete Institute, Detroit, Michigan 1974.
[31] Petersson, P-E.: Crack growth and development of fracture zones in plain con-
crete and similar materials, Report TVBM-1006, Lund Institute of Technology,
Division of Building Materials, Lund 1981.
[32] Smeplass, S.: Kalor { program documentation (in Norwegian). Sintef-report
STF65 A88031, Trondheim 1988.
[33] Stricklin, J.A., Haisler, W.E and von Riesemann, W.A.: Evaluation of solu-
tion procedures for material and/or geometrically nonlinear structural analysis.
AIAA Journal, Vol. 11, pp. 292-299, 1973.
[34] Thelandersson, S.: Modeling of combined thermal and mechanical action in
concrete, Journal of Engineering Mechanics, vol. 1123, No. 6, 1987.
[35] Taljsten, B.: Temperature development and maturity growth for Swedish Port-
land Cement Type II (in Swedish), Report 1987:035 E, Lulea University of
Technology, Division of Structural Engineering, Lulea 1987.
59
[36] Zienkiewicz, O.C., Emson, C. and Bettess, P.: A Novel Boundary Innite El-
ement, International Journal for Numerical Methods in Engineering, Vol. 19,
pp. 393-404, 1983.
[37] Zienkiewicz, O.C. and Taylor, R.L.: The Finite Element Method, Fourth Edi-
tion, Volume 1, McGraw-Hill, London 1989.
[38] Zienkiewicz, O.C. and Taylor, R.L.: The Finite Element Method, Fourth Edi-
tion, Volume 2, McGraw-Hill, London 1991.
60