Fibre Beam Element
Fibre Beam Element
Fibre Beam Element
IN
CASTEM 2000
J. Guedes, P. Pegon and A.V. Pinto
Applied Mechanics Unit, Safety Technology Institute,
Joint Research Centre, European Commission,
I-21020 Ispra (VA) Italy
The work presented in this report has been realized in the framework of the collaboration contract
num. 5090-92-11 TS ISP F between CEA/DMT/SMTS and JRC/IST/AMU for the development of
the computer code CASTEM 2000.
1 of 55
2 of 55
CONTENTS
1 INTRODUCTION 5
2 THE BASIC ASSUMPTIONS OF TIMOSHENKO BEAM THEORY AND
INTRODUCTION OF THE FIBRE MODELLING 7
2.1 Assumed displacement field 7
2.2 Strain-displacement relationships 8
2.3 Virtual work expression 8
2.4 The elastic case 9
2.5 Generalized fibre modelling 10
2.6 Usual fibre modelling 10
2.7 Cross-sectional warping 11
3 SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED
CONCRETE FIBRE MODELLING 13
3.1 Concrete Constitutive Laws 13
3.1.1 Monotonic loading 13
3.1.2 Cyclic loading 16
3.2 Steel constitutive law 18
3.2.1 Monotonic loading 18
3.2.2 Cyclic loading 19
4 FINITE ELEMENT IDEALIZATION FOR TIMOSHENKO BEAMS 23
4.1 Displacement and strain representation 23
4.2 Stiffness matrix evaluation 24
4.3 Element stress resultant and nodal resultant 24
4.4 The mass matrix 25
5 THE TIMO ELEMENT WITH FIBRES IN CASTEM 2000 27
5.1 Overview 27
5.2 Short remarks about CASTEM 2000 28
5.3 Global and local referentials for structural beam elements 29
5.4 The two-level approach 29
5.4.1 Description of the section. 29
5.4.2 Description of the beam 30
6 IMPLEMENTATION NOTES 33
6.1 Definition of the elements TRIS and QUAS 33
6.2 Models 33
6.3 Identification subroutines 34
6.4 Complex operators working for the beam (POUTRE) formulation 35
6.4.2 Element elastic stress calculation (SIGM operators) 35
6.4.3 Stiffness calculation (RIGI operator) 36
6.4.4 Mass matrix calculation (MASS operator) 36
6.4.5 Initialisation (ZERO operator) 37
6.4.6 Non linear flow evaluation (ECOU operator) 37
6.4.7 Tangent Stiffness (KTAN operator) 38
A FiBRE/TIMOSHENKO BEAM ELEMENT IN CASTEM 2000
3 of 55
4 of 55
INTRODUCTION
1. INTRODUCTION
Seismic analysis of Reinforced Concrete (R/C) structures requires realistic models able to represent
non-linear stiffness and strength characteristics of members under cyclic loading.
Many analytical models have been proposed in the past for the different structural elements. They
can be classified into two main groups: a) detailed models based on the mechanics of solids and the
consequent description of the local material behaviour (microscopic approach) and b) models based
on a simplified idealization of the overall behaviour (macroscopic approach).
Examples of the second group are the heuristic Takeda-type models able to represent uniaxial bending behaviour of R/C beams. In the first group are included the classical finite element models (continuum element support) and the semi-finite element models (discretization at section level) - fibre
type models (beam element support).
While the finite element models are a powerful tool for simulating the non-linear behaviour of complex structural parts (e.g. joints) their application to the common structural assemblages are impractical due to the large computational and memory requirements. On the other hand, the fibre type
modelling (see Fig. 1), taking advantage of the simplified kinematic hypotheses of the Euler-Bernoulli or Timoshenko theories offers a reliable and practical solution for the non-linear analysis of
composite structural elements such as building columns and walls and bridge piers. Moreover, with
this intermediate modelling approach it is in general possible to overcome the numerical difficulties
associated with the non-linear behaviour of concrete.
This report presents a fibre Timoshenko beam element implemented in the general purpose computer
code CASTEM 2000 [1]. This code allows the user to manipulate meaningful data-structures
(objects) by mean of a high level language (GIBIANE). The consequence of this unusual approach
Reinforced concrete
column
Beam assumption
and fibre modelling
Finite element
discretization
5 of 55
INTRODUCTION
(at least in the development and use of large structural analysis codes) for the implementation of
models has been investigated in [2].
The Timoshenko beam theory, assuming that plane sections remain plane after deformation but not
necessarily normal to the beam axial axis, has been considered in order to allow the development of
the so called generalized fibre model in the sense that very complex interactions between normal
force, bending moment and shear can be taken into account. A suitable constitutive model is still to
be derived, however the general formulation and the corresponding numerical data structure have
been created. In the mean time, uniaxial material constitutive relationships (usual fibre modelling)
for concrete and steel have been implemented.
Concrete behaviour is represented by a parabolic curve up to the peak stress point followed by a
straight line in the softening zone. Confinement is taken into account by the modification of the
plane concrete curve and including an additional plateau zone. Cyclic behaviour accounts for stiffness degradation and crack closing phenomena. Tensile resistance has also been considered.
Steel behaviour is represented by a modified Menegoto-Pinto model with a three-stage monotonic
curve (linear, plateau and hardening). Bauschinger effects are taken into account and buckling effects
can be simulated.
It is to be noticed that the present report intends to be a general document able to both describe the
work developed (finite element, constitutive models, etc.) and to assist in the implementation of new
models (implementation manual) in CASTEM 2000. In order to facilitate the access to the document, modelling and implementation aspects have been clearly separated as follows:
section 2, 3 and 4 are dedicated to the Timoshenko beam element and to the fibre modelling
including the material constitutive relations for reinforced concrete elements,
sections 5 and 6 deal with the specific aspects of the implementation in CASTEM 2000. A short
description of the code features is also included,
section 7 presents some numerical applications and Annex 1 one of the GIBIANE input,
finally the main conclusions are drawn and needs for further improvement are identified.
6 of 55
THE BASIC ASSUMPTIONS OF TIMOSHENKO BEAM THEORY AND INTRODUCTION OF THE FIBRE MODELLING
(1)
u z ( x, y, z ) = U z ( x ) + y x ( x )
Note that the planar normal rotations z and y are equal (except for the sign) to the slope of the
Ox axis corrected by the rotations y and z which are due to the transverse shear deformation:
z =
dU y
y
dx
y = z
dU z
dx
(2)
(3)
y
z
x
7 of 55
THE BASIC ASSUMPTIONS OF TIMOSHENKO BEAM THEORY AND INTRODUCTION OF THE FIBRE MODELLING
d y
dz
u x
dU x
=
+z
y
dx
dx
x
dx
d x
d x
u x u y
dU y
+
=
z z
= y z
dx
dx
y x
dx
(4)
d x
d x
dU z
u x uz
+
=
+ y + y
= z + y
dx
dx
dx
z x
Evaluating the appropriate partial derivatives, it is found that the three other components y , z and
yz of the strain tensor are apparently null. In fact, no provision was made at the displacement field
level (1) to account for these lateral deformations since the lateral stresses y , z
and yz are
assumed to be null; the values of these deformations can be eliminated from the expression of the
constitutive law.
( x x + xy xy + xz xz )dS dx ( Uy qy + Uz qz )dx
0S
= 0
(5)
where the first term takes the following form using (2) and (4):
L
d x Ux + z d x y y d x z x + y xy + z xz + d x x ( zxy + yxz ) dS dx
d
0S
x dS
(6.1)
shear forces:
Fy =
xy dS
S
Fz =
xz dS
(6.2)
twisting moment:
Mx =
(6.3)
8 of 55
THE BASIC ASSUMPTIONS OF TIMOSHENKO BEAM THEORY AND INTRODUCTION OF THE FIBRE MODELLING
bending moments:
My =
zx dS
M z = y x dS
(6.4)
Fx d x Ux + Fy y + Fz z + Mx d x x + My d x y + Mz d x z dx
d
(7)
( U y q y + U z q z )dx = 0
0
(8)
E
= -------------------------------------( 1 + ) ( 1 2 )
x
1
1 y
1
z
(9.1)
and:
xy
xz
yz
xy
1 0 0
E
= -------------------- 0 1 0 xz
2(1 + )
0 0 1
yz
(9.2)
From the condition y = z = yz = 0 , the following expression are obtained for y , z and yz :
y = z = x
yz = 0
xy = G xy
xz = G xz
(10)
9 of 55
THE BASIC ASSUMPTIONS OF TIMOSHENKO BEAM THEORY AND INTRODUCTION OF THE FIBRE MODELLING
the generalized stress vector F = ( F x, F y, F z, M x, M y, M z ) , namely the normal force (6.1), the
shear forces (6.2), the twisting moment (6.3) and the bending moments (6.4) take the expression:
Es dS Es zdS
S
Fx
Es zdS Es z
My =
Mz
Es y
Cz
dS
(11)
G s zdS
S
Fy
Ex
dS Es yzdS C y
E s ydS E s yzdS
Gs dS
Mx
Fz
E s ydS
Gs dS
Gs ydS
Ey
Ez
Cx
where no assumption about the position of the intersection R of the reference axis of the beam and
the section (see Fig. 2) nor on the orientation of the local axis Ry and Rz have been made.
(12)
In this report we considered only usual fibre modelling. Examples of such modelling for a reinforced
concrete beam are presented in the next section.
10 of 55
THE BASIC ASSUMPTIONS OF TIMOSHENKO BEAM THEORY AND INTRODUCTION OF THE FIBRE MODELLING
s xy dS
y
s xz dS
z
Fz =
Mx =
y
y
s z xy
(13)
z
s y xz )dS
s Gs dS
y
Mx
Fy
Fz
s G s zdS
z
s G s dS
y
s G s zdS
z
s G s ydS
z
s G s ydS
z 2
Gs ( s y
y 2
s z )dS
Ey
Ez
(14)
Cx
11 of 55
THE BASIC ASSUMPTIONS OF TIMOSHENKO BEAM THEORY AND INTRODUCTION OF THE FIBRE MODELLING
12 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
(15)
-------- = 1.0 + Z ( co )
co
(16)
13 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
co
co
Unconfined
(Z)
co
Confined
(Z*)
co*
the concrete. This phenomenon increases the peak value of the compression strength modifying its
position during the deformation process, as illustrated in Fig. 3. Notice that all the strain and stress
values presented in this report follow the rule: negative sign for compression and positive sign for
tension.
According to the Eurocode specifications [6] and Tassios [7], this effect can be taken into account in
the model through a confinement parameter depending on the section characteristics:
co = co
2
co = co
(17)
0.85
Z = ------------------------------------------------------------------ ( 0.1 + 0.0035 + co )
(18)
in which
= min ( 1 + 2.5 , 1.125 + 1.25 )
where represents the mechanical volumetric ratio of the stirrups equal to:
A s y ( l s )
= -----------------------------------------b c h c co
with A s and l denoting the cross-section of the leg and the total length of the stirrups, s the distance between stirrups along the member axis, b c and h c the dimensions of the confined concrete
core measured from the centre-line of the stirrups and y the yield strength of the stirrups. The
coefficient expresses the effect of both the longitudinal bars and the density of the stirrups on the
degree of confinement of the concrete core:
8
s
s
= 1 ------ 1 -------- 1 --------
3n
2b c
2h c
being n the number of longitudinal bars on the perimeter of the cross-section, engaged in the angle
of a stirrup.
With this parameter it is possible to represent, in a proper way, both the confined and unconfined
concrete, allowing a good description of phenomena like the spalling of concrete.
14 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
Tension law
co*
tm
pt
co*
Compression law
In order to improve the model for confined concrete, a third branch is considered after the softening
compression branch and before failure: a zero slope straight line defining a compression plateau (see
Fig. 4). This additional condition accounts for the residual strength of the concrete core for important
axial post-peak deformations.
The strength for this new branch is a parameter of the model, pt . Park and Priestley[8] suggest
pt = 0.2 co .
3.1.1.2 Tensile stresses
Tensile experimental tests on concrete specimens are not easy to perform and indirect methods may
be considered in measuring the strength of the material, usually less than 20% of the compressive
strength.
For maximum tensile stress, the concrete open cracks and the strength reduces suddenly to almost
zero strength. Because of this non ductile behaviour together with low-strength values, the tensile
strength of the concrete is usually ignored. However, in order to obtain a more realistic response for
structures under earthquake type loading, the tensile behaviour shall be taken into account. The
model herein adopted assumes a tensile stress-strain curve idealized by a straight line, with a slope
equal to the initial compression Young modulus, up to the tensile strength.
When R/C structural elements are considered, a new phenomenon is present due to the bond between
the longitudinal bars and the concrete: the tension-stiffening effect, i.e., the ability of intact concrete between adjacent cracks to carry on tensile stresses, giving to the concrete a fictitious supplementary post-cracking strength (smeared cracking approach).
Following the above considerations, a bilinear stress-strain curve has been adopted. From zero to
maximum tensile strength, point ( t, t) , the model presents a linear elastic behaviour with a slope
equal to the initial compression Young modulus, ( E o = 2.0 ( co co ) ) . The second branch, allowing a softening behaviour after cracking, i.e., after the maximum tensile strength is reached, presents
a straight line from point ( t, t) to a zero stress point at tm , as indicated in Fig. 4.
The constitutive equations are then given by:
= Eo
(19)
15 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
(20)
r = tm t
(21)
Numerical tests performed on R/C elements show that during the initial loading steps tensile stresses
may assume an important role on the structures global behaviour and that a non-smooth response
curve is obtained if sudden cracking is assumed, ( tm = t ) . The second branch of the curve, representing a post-cracking softening behaviour, not only accounts for the tension-stiffening effect but
also contributes to give better convergence characteristics to the solution algorithm.
3.1.2 Cyclic loading
Under dynamic earthquake type loading conditions, the concrete is submitted to loading and unloading processes. Thus, the non-linear model presented in Section 3.1.1, applicable only to monotonic
stress-strain type loading, must be extended to cyclic loading.
Experimental results show that the envelope for the concrete stress-strain curve under cyclic loading
may be considered unique and identical to the monotonic one. Unloading from the envelope reflects
the stiffness degradation due to monotonic damage and the reloading curve from zero stress tends to
the envelope evidencing some strength degradation. The intersection point between the two curves
defines the so-called common point. Stresses above it produce additional plastic strains, while cycles
below the line defined by the locus of the common points result in stable loops.
A small hysteresis effect is present during the cyclic process, as well as a decrease in strength, i.e.,
for the maximum strain ever reached during the loading history the reloading curve assumes a stress
lower than the one given by the previous unloading branch. Along with this phenomenon, the concrete presents, for both the unloading and the reloading branches, decreasing stiffness values as the
maximum strain increases on the stress-strain curve (see Fig. 5).
Under tensile cyclic loading, the envelope for the concrete may also be considered unique and identical to the monotonic curve. Before reaching the tensile strength the concrete presents an almost linear elastic behaviour. After that point, experimental results on R/C elements show that the unloading
and reloading branches follow a straight line and present increasing stiffness degradation for increasing maximum strains on the envelope.
Common points
Envelope
16 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
No experimental tests have yet been performed subjecting the concrete simultaneously to compressive and tensile stresses under cyclic loading. This lack of information prevents a good knowledge of
the crack closing phenomenon, i.e., the way the concrete re-starts to work in compression when, after
crack opening, there is a reversedl strain loading. Nevertheless, some considerations on this subject
will be presented later (see Section 7.3.1).
The cyclic behaviour law adopted in the present model is based on experimental results. The compression monotonic curve will be the envelope of the concrete behaviour law under compression
cyclic loading.
In order to simplify the model, unloading from the envelope follows a law similar to the one proposed by Mercer [9], i.e., a straight line with a slope depending on the maximum strain ever reached
during the loading history, max :
2
( max co )
-2
E d = E o 1.0 ---------------------------------------------------------------------------
1.0 + ( max co ) + ( max co )
(22)
The decrease of the slope with the increase of the maximum strain (equation (22)) tries to represent,
in a simple way, the degradation of the material stiffness. The reloading compression curve, also a
straight line (see Fig. 6), goes from the zero stress point at the plastic strain, pl , to the last loading
point reached on the envelope. No strength degradation is considered.
Concerning tensile stresses, the model considers an envelope going from zero stress point at co (see
Fig. 7) to the tensile strength point on the monotonic curve, ( t, t) . Loading in tension follows a
straight line so that the maximum stress is reached for ( tr1 = pl + t ) . The second branch, a
straight line defining the post-cracking behaviour, presents a slope such that the intersection with the
zero stress line is achieved at ( tr2 = pl + tm ), as shown in Fig. 7.
Strain tm is also a parameter of the model. Barzegar-Jamshidi and Schnobrich [10] suggest that
when the reinforcing steel intersects the cracks at right angles, an appropriate value for tm is the
yield strain of the longitudinal reinforcement steel bars.
Any unloading from the tensile curve follows a slope equal to the secant modulus of the concrete at
the reversal stress point.
With this model, strength and stiffness degradation of the concrete under tensile stresses is implicitly
taken into account. Unloading after tr2 follows a zero stress path. When the maximum tensile
3
7
9
= f()
Figure 6- Numerical model for the concrete behaviour under cyclic loading.
17 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
Envelope
= f()
2
pl
3
4
co
tr1
1
7
tm
tr2
t
tm
strength is reached, no more tensile stresses can be supported by the concrete in subsequent unloadings from compression.
If no tensile stresses are considered in the model, or its maximum value has already been reached in
previous cycles, the plastic strain is determined by the zero stress point from the compression
unloading curve. Otherwise, if during the unloading no cracks were present in the concrete and at the
end the strain goes beyond tr2 , the plastic strain defined by the compression unloading curve is corrected by t , ( pl = pl + t ). Fig. 7 illustrates a possible strain history starting from an unloading
compression curve.
18 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
In terms of modelling, the monotonic tensile stress-strain curve proposed for the reinforcing steel
bars is a five parameters law, including the following three zones: the linear elastic, defined by the
Young modulus, E , and the yield strain, sy ,
= E
(23)
the yielding plateau from the yield strain, sy , up to the hardening strain, sh , and defined by the
yield stress, sy ,
= sy
(24)
and finally the hardening curve up to the ultimate stress-strain point, ( su, su) , represented by a
fourth degree polynomial:
su 4
= su ( su sy ) --------------------
su sh
(25)
The same constitutive law is adopted for compression stresses if no buckling effect is considered for
the reinforcing bars. Otherwise, some modifications must be introduced after the linear elastic zone.
This issue will be referred to later on when presentating the reinforcing steel cyclic behaviour.
3.2.2 Cyclic loading
3.2.2.1 basic model
According to experimental results, steel bars under cyclic loading present an unloading straight line
pattern with a slope equal to the initial Young modulus. If after loading beyond the yield strain
reversed stresses are applied to the bars, the stress-strain response curve becomes nonlinear at a
stress much lower than the initial yield stress. This phenomenon, the so-called Bauschinger effect,
depends on the strain history.
Following the experimental results, the monotonic curve together with the explicit formulation proposed by Giuffr and Pinto, later on implemented by Menegotto [12], were chosen to represent the
behaviour of the reinforcing steel bars under cyclic loading.
If small unloading cycles are imposed to the bars, no hysteretic effect is evidenced and a straight
line, with a slope equal to the elastic Young modulus, can be adopted for both unloading and reloading curves. Under virgin strains, reloading follows the monotonic curve.
Otherwise, if after unloading from a post-yielding position the strain satisfies the relation:
sy
max < --------3.0
the Giuffr curve is activated and an important hysteretic effect may occur. Any subsequent loading
or unloading curve follows this formulation.
The selection of this model to represent the steel cyclic behaviour was based on its simplicity,
numerical efficiency and good agreement with experimental results from cyclic tests on reinforcing
steel bars. The simplicity is well expressed by the fact that a single equation is enough to represent
both loading and unloading stress-strain curves:
(1 b)
= b + ----------------------------------
1
R
R
( 1 + ( ) )
(26)
19 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
where
s r
= -------------o r
s r
= ----------------o r
a1
R = R o -------------a2 +
As illustrated in Fig. 8, this equation defines a family of transition curves between two asymptotes
with slopes E o and E h and having ( o, o) as a common point. The pair of values ( r, r) are the
coordinates of the last reversal loading point. On the other hand, the b factor represents the slope
hardening ratio, i.e., the ratio between the hardening slope, E h , and the initial slope, E o , and R is a
parameter defining the shape of the transition branch of the curve. This parameter allows a good representation of the Bauschinger effect and its value depends on , i.e., the difference between the
maximum value of the strain ever reached on the loading direction and o , normalized by ( o r ) .
Numerical tests performed on several structures showed that when the tangent stiffness is considered
in the resolution algorithm the convergence process can be very sensitive to the adopted slope hardening ratio. A small b factor, together with a softening behaviour in some of the concrete fibres, may
introduce singularities in the global stiffness matrix.
Nevertheless, the hardening slope can be estimate by:
su sy
E h = -------------------- su sy
(27)
Parameters a 1 and a 2 , together with R o , should be obtained from experimental results. For all the
numerical tests performed, the values suggested by Menegotto [12] were adopted:
R o = 20.0
a 1 = 18.5
a 2 = 0.15
To start the cyclic loading model, the stress and strain values for both reversal points are initialized
by the yield parameters.
2x(o-r)2
(o,o)2
(r,r)1
(o,o)
Eh
= f()
R(2)
Eo
Eo
R(1)
(r,r)2
Eh
(o,o)1
(r,r)
1x(o-r)1
20 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
b c = a ( 5.0 L D )e
Eo
b------------------ sy
(28)
sy
in which = 4.0 ----------- and is the major difference between the reversal strain values found in
LD
both directions during the loading history.
On the other hand, during reloading after reversal from compression, a reduced Young modulus is
adopted for the first asymptote:
E r = E ( a 5 + ( 1.0 a 5 )e
( a 6 )
2
(29)
in which a 5 = 1.0 + ( 5.0 L D ) 7.5 and = ( o r ) , and consequently a new b factor for
tensile stresses must be adopted,
b t = ( bE ) E r
(30)
If 5 is the strain at which the compression curve diverges more than 5% from the tensile one
towards lower values, an additional stress given by
b bc
s = s E o -----------------1.0 b c
(31)
= f()
E
s
bcxE
bxE
Er
21 of 55
SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING
where
11.0 L D
s = o 5 = ----------------------------------------c(L D)
10 ( e
1.0 )
(32)
must be added to the yield value so that the new hardening asymptote could be correctly positioned. This shift allows the response curve to follow first the unbuckled curve and then the buckling
compression branch (see Fig. 9).
Expressions (27) to (31) are valid for ( 5.0 < L D 11.0 ) . Parameters a , c and a 6 are experimental
values to be defined for each steel grade. In all the numerical tests presented in this work the values
suggested by Monti and Nuti [11] were adopted:
a = 0.006
c = 0.500
a 6 = 620.0
22 of 55
U ( x ) = N ( x )U + N ( x )U
1
(33)
2
where U and U are the nodal displacements and N and N are the conventional linear shape
functions:
x2 x
1
N ( x ) = ------------l
x x1
2
N ( x ) = ------------l
(34)
( x ) = N ( x ) + N ( x )
1
(35)
Denoting by the set of nodal variables ( U , , U , ) and using (2), (8), (33), (34) and (35)
the expressions for the generalized strains in the element can be readily evaluated:
Ex =
Ey =
dUx
1 1 1 2
s
= --- U x + --- Ux = B x
l
l
dx
dU y
1 1 1 1 1 2 1 2
s
z = --- U y --- z + --- U y --- z = B y
l
2
l
2
dx
(36)
dU z
1 1 1 1 1 2 1 2
s
Ez =
+ y = --- U z + --- + --- U z + --- y = B z
2
l
2 y l
dx
C =
d
1 1 1 2
b
= --- + --- = B
l
l
dx
= x, y, z
In fact, the expressions of the generalized transverse shear strains, say E y and E z have been modified (elimination of the linear terms according to Donea and Lamain[13]), in order to a priori remove
the problem of shear locking which occurs for equal approximation elements. A discussion of this
problem has been given in a previous report dealing with the usual linear Timoshenko element [14].
The equations (36) take the following symbolic form
E = B
(37)
23 of 55
E
K
Edx
(
U
q
+
U
q
)dx
=
0
y y
z z
(38)
Es dS
K 15 =
K 22 =
Es zdS
K 16 = E s ydS
s Gs dS
K 24 = s G s zdS
K 33 =
s Gs dS
z
K 34 =
K 42 =
y
s G s zdS
s Gs ydS
z
K 43 =
z
s G s ydS
K 51 =
Es zdS
(39)
K 44 =
z 2
ES ( s y
K 55 =
Es z
dS
K 56 = E s yzdS
K61 = E s ydS
y 2
s z )dS
K 65 = E s yzdS
K 66 =
Es y
dS
(40)
The contribution K of the element to the global stiffness K takes as usual the following form:
x2
e
K =
[B]
K B dx
(41)
x1
2 2
(42)
which, unlike the Euler-Bernouilli cubic Hermitian element, only has lateral nodal point forces and
no nodal moments.
24 of 55
transverse shear forces F y and F z are also constant in the element according to the expression (36)
of the transverse generalized strains E y and E z .
Knowing this generalized stress state ( F x, F y, F z, M x, M y, M z ) it is possible to determine the contri1
1
2
2
bution of the element to the nodal resultant ( F , M , F , M ) of the element, using again the virtual
work principle (7), the discretization (36) of the generalized strain and a one-point Gauss-Legendre
rule:
1
F = F
= x, y, z
F = F
Mx = Mx
Mx = Mx
l
1
M y = --- F z M y
2
l
1
M z = --- F y M z
2
(43)
l
2
M y = --- F z + M y
2
l
2
M z = --- F y + M z
2
d ux
d uy
d u z
W I = u x 2 + u y 2 + u z 2 dSdx
dt
dt
dt
S
(44)
using the expression (1) of the assumed displacement field and introducing the generalized displacement vector U defined by (3), W I in (44) may be rewritten under the form:
L
WI =
dU
dt
(45)
dx
s dS
M15 =
M 22 =
s zdS
S
s s dS
M 24 = s s zdS
M 33 =
s s dS
z
M 34 =
M 42 =
y
s s zdS
M 51 =
M 43 =
z
s s ydS
M 55 =
M 61 = s ydS
s s ydS
z
s zdS
M 16 = s ydS
M 44 =
z 2
S ( s y
s z
dS
M 56 = s yzdS
M 65 = s yzdS
(46)
y 2
s z )dS
M 66 =
s y
dS
25 of 55
Introducing in (45) the discretization (33)-(34)-(35) leads to the addition of the linear term
2
2
M ( d dt ) to the first member of the discrete equilibrium equation (40). Performing exact integration for all terms:
x2
x2
x2
l
( N N )dx = ( N N )dx = --3-
(N
x1
x1
x1
2
l
N )dx = --6
26 of 55
( U, ) E
F ( F, M )
BEAM level
SECTION level
27 of 55
Evaluation of the generalized strain E at the integration point of each beam element, from the
1
1
2
2
nodal generalized displacement ( U , , U , ) .
Use of the beam model in order to evaluate the strain tensor , and in particular its normal component x at the level of each fibre, located at the Gauss integration points of the elements describing the section.
Use of the constitutive relationship in order to evaluate the stress tensor at the level of each
fibre, and in particular its normal component x .
Integration over the section of the relevant stress components in order to compute the generalized
stress F for the section (the use of Gauss numerical integration rules allows to obtain precise evaluations for these quantities).
1
28 of 55
tensor. For the same geometrical sub-zone, the situation is quite different when changing the formulation of the element. Consider for instance the finite element COQ8, associated with the same 8node geometric element but corresponding now to a 3-D structural thick shell formulation: not only
the number (and the name) of the active components are different (5 components for the stress tensor
in the local axes), but the integration points are arranged in such a way that an integration over the
thickness of the shell is also allowed.
To create the field-by-elements and to further check the internal data consistency, CASTEM 2000
uses internal identification routines, one for each specific field-by-element type: these routines
return, for a given model and a given finite-element formulation, the number and the name of the
components as well as the type of the value.
CASTEM 2000 was developed using the ESOPE language. ESOPE [15] is in fact an extension of
FORTRAN 77, allowing the manipulation (creation, re-sizing, activation, inactivation, destruction)
of data-structures.
29 of 55
mesh is to be generated. Note also that the degree of eccentricity of a given section is simply controlled by the position of the section with respect to the origin of the XOY reference.
The fibres state (e.g. the local stress ) is described by standard objects of the type field-by-element associated with the object of the type model.
A new formulation and accordingly new finite elements are required for at least two reasons:
The number of active components of the stress tensor ( x , xy and xz ) on the section reduces to
3.
the fibres being located at the Gauss integration points of the elements, the useful integrations
over the section, say for deriving the generalized stress state or the section stiffness matrix, are
performed numerically looping over the elements of each sub-zone of the section model and using
the standard Gauss integration rules. Due to the nature of the expression to be integrated (see e.g.
(39)), the full integration (that is, an integration which gives exact results for an elastic case),
should be performed with a number of integration points which may differ from the number used
for instance in the continuum element associated with the same geometric element. This is in particular the case for the TRIS element which requires more than the one integration point in consistent use for the TRI3 (continuum) element.
Only two new elements are useful for the following two reasons: first it is possible to mesh any 2-D
section by mean of triangular or quadrilateral elements, and second neither nodal quantities nor
shape functions being really used, unless for computing the coordinates of the integration points,
there is no need to introduce any element of order higher than one (linear element).
The effective specification of the modelling of the section is completed after the specification of the
material and geometrical characteristics (object of the type field-by-element of the specific type
characteristics). The material characteristics are the usual elastic parameters ( E , and ) and
eventually the non-linear parameters of the constitutive laws of the fibres. The geometrical characteristics are the two cross-sectional warping parameters of the fibre.
In order to check the physical consistency of a given section, an operator called MOCU has been
developed. This operator uses as input the model and the characteristics of the section which defines
a beam, and three objects of the type list-of-real-numbers of same length associated with two
sequences of curvatures C y and Cz and one of normal force F x . These sequences define a bi-axial
loading path for the considered beam. The main output of MOCU are the related two sequences of
bending moments M y and, M z , and the sequence of normal deformation E x .
5.4.2 Description of the beam
A frame of beams with complex sections is described by an object of the type model, associated to
a mesh of 3-D linear structural beam element TIMO[14]. This element, associated with the beam
(POUTRE) formulation, involves one integration point corresponding to the mid-section of the element.
The field-by-elements relevant for the modelling, namely the material or geometric characteristics,
strain, stress and internal variables field are, as usual, a collection of arrays associated to each component of the each field. These arrays contain information relative to the section of integration.
Instead of being simple objects belonging to the real-number, vector, curve or list-of-real-numbers classes, the content of the array may be the complex objects of the type model or field-byelement which effectively describe the section of integration.
The resort to these more complex informations are useful for the TIMO element with fibres:
the material characteristics of each element of a sub-zone are the model (component MODS)
and the field-by-element of characteristics (component MATS) of the section of integration,
30 of 55
the internal variables of each TIMO element are the objects of the type field-by-element and of
the specific type stress (component CONS) and internal variables (component VAIS) on
the section.
These informations are extensively used when a sub-zone of the beam model is non-linear. In order
to limit the computational effort in the case of elastic sections, an optional component have been
added to the field-by-element of material characteristic: the component name is MAHO and the
associated value is a list-of-real-number containing the section constitutive matrix K (39) (Hook
matrix). This matrix can be computed using the HOOK operator on an equivalent beam element
and by extracting from the resulting field-by-element of the specific type Hook-matrix the value
of the component MAHO at the integration point of this element.
It would have been possible to address the problem of the elastic section in a maybe more general
way by introducing a new beam element and consequently a new formulation. This solution has not
been retained in order to limit the duplication of coding effort between this new element (and its
associated formulation) and the already existing TIMO element with beam (POUTRE) formulation.
31 of 55
32 of 55
IMPLEMENTATION NOTES
6. IMPLEMENTATION NOTES
6.1 Definition of the elements TRIS and QUAS
The name of the elements TRIS and QUAS with the name of the related formulation SECTION
are defined in BDATA. TRIS and QUAS are declared in the NOMTP array, whereas SECTION is
declared on the array NOMFR.
The new relations element/formulation and element/geometry are introduced in the two subroutines
NUMMFR and NUMGEO. The value of the index MELE associated with the element TRIS and
QUAS are respectively 104 and 105 and the value of the index MFR associated with the SECTION formulation is 47.
The other informations relative to the elements (e.g. number of integration points, associated geometric element) are introduced, as usual, in the subroutine ELQUOI. The array INFELE which is
constructed by this subroutine is consulted before every computation at element level.
Recall that for the TIMO element: MELE=84 and MFR=7.
6.2 Models
For the SECTION formulation, the models are assumed to be at least elastic isotropic (ELASTIQUE ISOTROPE). The non linear models for the section are introduced by the word PLASTIQUE followed by the name of the non-linear fibre model. The two models already presented (see
Section 3) are available under the name ACIER_UNI for the steel and BETON_UNI for the concrete.
For the beam (POUTRE) formulation, the models are assumed to be at least elastic section
(ELASTIQUE SECTION). Only one non linear model is obviously available: its name is PLASTIQUE SECTION.
At the implementation level, the highest entry point for the modelling, directly associated with the
operator MODE, is the subroutine MODELI. This routine would be informed of all the default in
use. MODELI calls the subroutine MODEL2 which is responsible for the mechanical (MECANIQUE) modelling in use here. The names of the new elements that can be associated with such
modelling should be added here. This routine checks for the existence of elastic modelling using the
subroutine MODELA and then for a non-linear modelling using the subroutine MODNLI, which recognizes the PLASTIQUE behaviour, and eventually branches to the routine responsible for such a
behaviour, MODPLA, which analyzes the possible names of the PLASTIQUE models.
At the end of MODEL2 the name of the elements TRIS and QUAS has been added in the array
LESTEF,
In MODELA a new slot has been added for holding the elastic SECTION behaviour (relative
index IELAS=7).
In MODPLA 3 new slots have been added for the 3 new plastic constitutive behaviours:
BETON_UNI (relative index INPLAS=20), ACIER_UNI (INPLAS=21) and SECTION
(INPLAS=22).
According to these modifications/additions, the manual of the operator MODE has been also
changed.
The routine NOMATE associates to each sub-zone of a given model a character string CMATE, characteristics of the elastic part (say isotropic, anisotropic or... SECTION), and an absolute index
INATU for the non linear material. This routine has been expanded in order to account for the previous modelling: at the SECTION level CMATE=ISOTROPE (isotropic) and INATU=0 (elastic
33 of 55
IMPLEMENTATION NOTES
Geometric characteristics: the names of the two cross-sectional warpings s (ALPY) and s
(ALPZ) are introduced in IDCARB.
Stress: the names of the active stress components x (SMXX), xy (SMXY) and xz
(SMXZ) are introduced in IDCONT.
Internal variables: the names of the internal variables of the two fibre models are introduced in
IDVARI: 6 internal variables for BETON_UNI and 12 for ACIER_UNI. In both models two
internal variables are TANG (current uniaxial stiffness) and EPSO (current uniaxial strain).
The values of all these fields belong to the real-numbers class.
For the beam (POUTRE) formulation some new features are introduced for the following fieldby-elements:
Material characteristics: The material characteristics are defined in the subroutine IDMATR.
They are the section model (MODS, type model), the section characteristics (MATF, type
field-by-element) and, as an optional component, the section constitutive matrix K (39)
(MAHO, type list-of-real-number)
Internal variables: the internal variables for a non linear section (IDVARI) are the stress state of
the section (CONS, type field-by-element) and the associated internal variables (VAIS, type
field-by-element).
34 of 55
IMPLEMENTATION NOTES
The manual of the operator MATE has been changed according to the new possible characteristics.
35 of 55
IMPLEMENTATION NOTES
Use the SIGM operator with the field-by-element of geometrical and material characteristic and
the field-by-node of (generalized) displacement.
Note that the last syntax is the most widely used for instance when performing the elastic predictor
step in the non-linear model computations. However, even in this case, the Hook matrix (which is
not known as input) is internally computed, which means, for the TIMO element with fibres and as
shown in the previous subsection, the need for a work at the SECTION level. It was then important
to design a possibility to escape this call and this was the reason why the MAHO component was
optionally allowed in the field-by-element of material characteristics.
At the beam level the driver routine HOOK2P has been modified for the TIMO element with fibres in
order to account for complex material characteristics. The driver calls HOOK2D. When the first
syntax of SIGM is used or when the value of the component MAHO is available, the Hook matrix
is directly loaded. In the other case, a work at the section level is performed (call to FRIGIE and
FRIGI2) and the matrix is loaded using DOHTIF.
Since in contrast with the usual (isotropic) TIMO element, the TIMO element with fibres leads
generally to an Hook matrix which is not diagonal, two routines TIFSTR (driving the transformation
global/local for the displacement) and TIFOLO (computing the local element stress resultant) have
been added. These routines are called in SIGMA3 when CMATE=SECTION.
Recall that the components of the field-by-element of generalized stress have the following name:
EFFX= F x , EFFY= Fy , EFFZ= F z , MOMX= M x , MOMY= M y and MOMZ= M z .
6.4.3 Stiffness calculation (RIGI operator)
The operator RIGI calculates the elastic contribution of all the elements of a mesh to the global
stiffness matrix.
As for the SIGM operator, RIGI works with two syntaxes: the first one using a field-by-element
of Hook matrix and the second one a field-by-element of material characteristics. Although this is
obviously less strategic than for the SIGM operator (usually the elastic stiffness is computed only
one time when a model is set), it is possible to prevent RIGI from activating calls at the SECTION level in the case of the second syntax when the component MAHO is available.
As usual the driver routine (RIGI1) is modified in order to account for complex material characteristics. RIGI1 calls RIGI4 which either directly loads the Hook matrix or computes it (FRIGIE,
FRIGI2 and DOHTIF). The element stiffness matrix is then computed by mean of two new routines: TIFILO (computing the local element stiffness matrix) and TIFRIG (driving the transformation local axis/global axis).
Note that the product of the stiffness matrix computed with RIGI and any global displacement field
should be equal to the global resultant, computed by BSIG, of the element elastic stress, computed
before using SIGM and associated with the same displacement field.
6.4.4 Mass matrix calculation (MASS operator)
The operator MASS computes the contribution of each element to the global mass matrix. This
operator is generally used at most one time after having setting up a model.
For the TIMO element with fibres, the driver routine MASS1 is modified (complex material characteristics). The routine MASS3 then activates the SECTION level (call to FRIGIE and FRIGI2)
in order to determined the section mass matrix which is loaded by DOHTIF. Two routines TIFALO
(computing the local element mass matrix) and TIFMAS (driving the transformation local axis/global axis) have been added.
36 of 55
IMPLEMENTATION NOTES
new
old
+ (
old
old
, )
new
(47)
old
new
old
+ (
old
old
, )
(48)
The SECTION level is then implemented using the subroutines FCOUL1, FCOUL2 and specialized routines for each kind of fibre. In fact FCOUL1 and FCOUL2 work for the section as
ECOUL1 and ECOUL2 work in general. FCOUL1 loops on the various sub-zones of the section
model and prepares the constant data for the zone, FCOUL2 loops on the elements and on the integration point(s) of each zone. The main difference between FCOUL2 and ECOUL2 comes from
the introduction of the coordinates of the integration points in order to retrieve, using (4) and the
generalized strain increment, the strain increment at the current point. In the same way, after having computed the new local section stress state using a specific constitutive model, the numerical
contribution to the generalized stress integrals (6.1)-(6.4) is derived using the standard Gauss formula. A computed GOTO in FCOUL2 allows to branch to each specific model depending on the
relative section model index INFIBR (ranging from 0). The transformation between the absolute
model index (set by NOMATE) and the relative section model is performed by the new subroutine
37 of 55
IMPLEMENTATION NOTES
TEMANF. The relative numbering of the non linear section models should be compliant with the
one introduced for the model name definition (MODPLA), the material characteristics (IDPLAS)
and the internal variables (IDVARI). As in ECOUL2, the slot 0 is reserved to the elastic model.
Three specific section models are available. The first one is purely elastic and the corresponding
subroutine is FIBELA. The second model concerns the concrete behaviour and the corresponding
subroutine is FIBETO which calls in turn NEWBET. NEWBET is a pure FORTRAN subroutine
which implements the 1-D concrete fibre model described in Section 3.1. The third model deals
with the steel behaviour. The associated subroutine is FIBSTE which calls in turn STEEL1.
STEEL1 is a pure FORTRAN subroutine which implements the 1-D steel fibre model described in
Section 3.2. FIBETO and FIBSTE are not only simple interface routines with the CASTEM 2000
data-structures and the dedicated routines which perform the modelling, they also complete elastically the behaviour with respect to the shear components xy and xz of the stress (see the discussion about generalized and simple fibre modelling in Section 2.5 and Section 2.6)
6.4.7 Tangent Stiffness (KTAN operator)
Due to the complex nature of the global section behaviour, the convergence of the algorithms introduced for performing the step-by-step non linear analysis requires a frequent evaluation of the tangent stiffness matrix. This evaluation is to be performed by the KTAN operator.
Due to the particular nature of the fibre models herein introduced (usual fibre modelling of chapter
2.6), the tangent constitutive matrix in the current configuration depends only on the tangent value
tang
tang
E
of the constitutive curve of all the fibres in this configuration. Substituting E by E
in (39),
tang
it is possible to derive the tangent section matrix K
which is then used for deriving the global tangent stiffness.
As in the implementation of the ECOU operator, the two levels of modelling are treated in a separate way.
The beam level requires the usual modifications of the subroutine KTANGA: convenient treatment of the special kind of component of the field-by-element of the specific type characteristics and internal variables for the TIMO element with fibres. In the same routine are
performed the nested loops on the elements and Gauss points of each TIMO with fibres subzone. The work inside the loops is arranged as follow: First is computed (FRIGTA) the set of the
12 section integrals which allows to derive (DOHTIF) the tangent section constitutive matrix
tang
K
(tangent Hook matrix). The computation then proceed as for the evaluation of the elastic
stiffness matrix (use of TIFRIG and TIFILO as in Section 6.4.3).
The SECTION level is implemented using the subroutine FRIGTA and FRIGT2. FRIGTA performs the loop on the various sub-zones of the section models and prepares the constant data for
the zone. FRIGT2 loops on the elements and on the integration points of each zone. Depending
tang
on the model, the value of the internal variable containing E
is retrieved and the numerical
contributions to the section integrals (39) are computed.
38 of 55
IMPLEMENTATION NOTES
also be considered. At the implementation level, the subroutine MOCUR first reads the input, computes the elastic constitutive matrix (FRIGIE and FRIGI2) and then equilibrates the section for
each point of the loading path using a simple secant method and the evaluation of the generalized
stress on the section (FCOUL1 and FCOUL2). A manual of the operator MOCU is also available.
39 of 55
IMPLEMENTATION NOTES
40 of 55
NUMERICAL APPLICATIONS
7. NUMERICAL APPLICATIONS
7.1 Introduction
In order to validate the fibre model, as well as the presented fibre constitutive laws, and to illustrate
its potential, a set of monotonic and cyclic numerical tests on R/C structural elements has been performed. According to the kind of loading, two main groups are considered: monotonic and cyclic.
Two tests are included in the first group:
The first one simulates the response of a rectangular section for seven different reinforcement
ratios. The concrete and steel stress-strain behaviour laws are defined according to the theoretical
curves assumed by Park and Paulay [16]. It must be noticed that the aim of this test is not to validate the constitutive laws but to check the fibre algorithm procedures in CASTEM 2000.
The second test, also performed on a rectangular R/C structural element, intends to show the
potential of the fibre model in the sense that, if a proper constitutive law is adopted, phenomena
like cracking can be well represented. The importance of the concrete tensile strength in the global response of R/C elements under flexural moments is also analysed.
In the second group, where almost all the cases could be included, two numerical tests are presented:
The first one, performed on a R/C cantilever T beam with B500 steel, tries to simulate the LNECs
experiment on specimen S1V3 as described by Pipa and Carvalho [17]. From the analysis of the
results, some considerations concerning the crack closing phenomenon are presented.
Classical global behaviour models usually fail when biaxial type loading is considered. In order to
validate the fibre model under this type of loading, the numerical simulation of one of the tests
referred by Bousias [18] is presented as second test. Some improvements are then suggested in the
Giuffr-Pinto model so that strength degradation at high strain level could be properly taken into
account.
41 of 55
NUMERICAL APPLICATIONS
0.0375
0.0375
0.0375
0.0250
0.0250
0.0125
0.0125
0.0250
0.0125
0.0125
0.0125
Moment/(b*d2) (MN/m2)
11.0
As
d/10
1
2
8.8
3
6.6
4
5
As
d/10
4.4
b
6
7
2.2
As
= -----bd
Curvature (1/m)
0.0
.00
.56
1.12
1.68
2.24
As
= ------bd
2.80
X1.E-2
42 of 55
NUMERICAL APPLICATIONS
Moment (KNm)
20.0
4 -> r = 22.7
1 -> r = 1.0
2 -> r = 1.5
16.0
3 -> r = 2.0
3
12.0
316mm
2
.031m
1
8.0
.25m
.019m
4.0
.25m
Curvature (1/m)
.0
.00
1.00
2.00
3.00
4.00
5.00
6.00
X1.E-3
Figure 12- Influence of the tensile parameters in the section global response.
Thus, when a section is subjected to an increasing curvature loading path, the moment-curvature
response before cracking is represented by an almost straight line. During this first step the R/C
behaves as linear elastic and the three results are identical. When the first set of the fibres Gauss
points reaches the tensile strength, cracks start opening and the curve inflects tending to the response
curve of the zero tensile strength case.
These results illustrate the importance of adopting a tensile softening behaviour after cracking, since
a non ductile behaviour law, ( tm = t ), gives the discontinuous response represented by curve 1.
This curve reflects the fact that when cracks open the tensile strength goes immediately to zero,
changing instantaneously the internal equilibrium of the section. Part of the tensile force is then
transferred to other concrete and steel fibres and another part is compensated by a movement of the
neutral axis, allowing a decrease of the global compression forces.
This problem, affecting the convergence of the nonlinear algorithm, can, of course, be minimized by
considering a greater number of fibre elements, but this will also increase the computation costs.
Thus, post-cracking ductile behaviour has been adopted in the model, allowing not only a more
smooth response, as presented by curves 2, 3 and even 4 in Fig. 12, but also the possibility to account
for phenomena such as tension-stiffening in the final response.
It should be noticed that Barzegar [10] suggestion produces a ratio ( r = 22.7 ) which seems to be too
high.
43 of 55
NUMERICAL APPLICATIONS
412mm
46mm
0.10m
312mm
6mm // 7mm
0.20m
0.20m
0.40m
The test performed on the specimen S1V3 is simulated by the fibre algorithm. All the concrete was
considered confined with a compressive strength ( co = 35.0 MPa ), at ( co = 0.2 %), and a slope
for the descending part of the envelope ( Z = 100.0 ) . A confinement factor ( = 1.33 ) and a Poisson ratio ( = 0.25 ) were also adopted. No compression plateau nor tensile strength was considered
in the concrete model.
For the steel we adopted the following characteristics: yield stress ( sy = 550.0MPa ), Young modulus ( E = 203.0GPa ), Poisson ratio ( = 0.25 ) , hardening and ultimate strain equal to 0.22 % and
9.0 %, respectively, ultimate stress ( su = 610.0MPa ) and hardening ratio ( b = 0.03 ).
Fig. 13 shows the cross section of the cantilever T beam and Fig. 14 the discretization adopted for
the numerical application. The structure was divided into 7 elements with increasing length from the
base to the free top, as follows: 0.15, 0.15, 0.15, 0.15, 0.30, 0.30 and 0.30m . The loading history, a
horizontal displacement, was imposed at the free top of the beam. No axial compression forces were
considered.
Force Y (KN)
60.0
40.0
Fibre
Experimental
20.0
.0
-20.0
-40.0
-60.0
Disp. UY (m)
-80.0
-.15
-.12
-.09
-.06
-.03
.00
.03
.06
.09
.12
.15
44 of 55
NUMERICAL APPLICATIONS
31.0mm
0.25m
816mm
Y
19.0mm
15.0mm
8mm // 70mm
0.25m
The response curve from the fibre model, illustrated in Fig. 14 by the free top force-displacement
diagram, presents a rather good agreement with the experimental results. The discontinuities
observed in the dashed line curve result from the buckling of some of the bottom steel bars caused by
the opening of the stirrups.
The major differences between both diagrams appear during the crack closing. The pinching effect is
in the case of the fibre model much more pronounced. Although some modifications could be performed in the future in order to improve the behaviour of the concrete under cyclic loading, the
present model seems to fail on representing in a realistic way the crack closing phenomenon.
7.3.2 Biaxial Loading on a R/C Square Section
The structure tested by Bousias [18], a R/C column element with 1.49m length, framing into an
enlarged end block of reinforced concrete, intended to simulate part of a column between the foundation and the point of inflection. From the twelve specimens tested, the specimen S9 was chosen. The
cross section, a square doubly symmetrical one, is presented in Fig. 15.
For the loading history, the axial force and the horizontal displacements imposed at the free top of
the column during the experimental test were adopted. These values were read directly from a data
base and are presented in Fig. 16. The displacement path is what is usually called a shrinking path,
since, after an important displacement imposed to the structure, there is a progressive movement, in
both directions, towards the initial zero displacement position. The loading follows a four squares
Disp. UZ (mm)
Force X (KN)
110.
140.
2
4
90.
70.
1
.0
50.
7
30.
5
Disp. UX (mm)
Disp. UY (mm)
10.
-140.
-2.0
.0
2.0
4.0
-140.
.0
140.
45 of 55
NUMERICAL APPLICATIONS
displacement path centered at the origin with half-side lengths equal to 100, 80, 60 and 40mm in
each direction.
The characteristics adopted for the materials are the mean values found from the laboratory tests performed on concrete and steel specimens.
Around all the section a concrete cover of ( c = 19mm ) thickness with a compressive strength
( co = 30.0MPa ) , at ( co = 0.2 %), and a slope for the descending part of the envelope
( Z = 100.0 ), was considered. For the concrete core a confinement factor ( = 1.33 ) was introduced. For the tensile parameters, a strength ( t = 3.0MPa ) and a ratio ( r = 3.0 ) were used. The
model was completed assuming a Poisson ratio ( = 0.25 ) and no compression plateau.
For the steel we adopted the following average characteristics: yield stress ( sy = 440.0MPa ) ,
Young modulus ( E = 203.0GPa ), Poisson ratio ( = 0.25 ), hardening and ultimate strain equal to
0.8 % and 13.0 %, respectively, and ultimate stress ( su = 760.0MPa ) . For the hardening ratio
( b = 0.13 ) was adopted.
In terms of finite element discretization, the structure was divided into 7 elements with increasing
length from the base to the free top, as follows: 0.05, 0.10, 0.15, 0.20, 0.25, 0.35 and 0.39m. The discretization adopted for the section is represented in Fig. 15.
The response curves from the fibre model, showed in Fig. 17 and Fig. 18, present a good agreement
with the experimental results. The reduction of strength along the cycles, not evidenced in the fibre
results, is due to the degradation of the concrete and steel properties. After the first square path
almost all the concrete is destroyed and the response curve is mainly controlled by the steel.
The high strength degradation verified in the section, which results from the high strain levels
imposed to the column, is not well represented by the fibre model; Since no strength degradation is
considered in the steel model, the fibre response curves from the second and third square loading
cycles tend to the same asymptote. In what concerns the present test , this characteristic of the Giuffr-Pinto model is responsible for the most important differences found between the experimental
and the numerical results.
Force Y (KN)
60.0
Fibre
Experimental
40.0
20.0
.0
-20.0
-40.0
Disp. UY (m)
-60.0
-.10
-.08
-.06
-.04
-.02
.0
.02
.04
.06
.08
.10
46 of 55
NUMERICAL APPLICATIONS
Force X (KN)
80.0
Fibre
Experimental
60.0
40.0
20.0
.0
-20.0
-40.0
Disp. UX (m)
-60.0
-.10
-.08
-.06
-.04
-.02
.0
.02
.04
.06
.08
.10
47 of 55
NUMERICAL APPLICATIONS
48 of 55
ACKNOWLEDGEMENTS
The authors would like to thank Alain Millard and Laurent Le Ber for their contribution (and hopefully a touch of rational elegance) to the last step of this work: the effective implementation of the
fibre approach in the official version of CASTEM 2000.
49 of 55
REFERENCES
REFERENCES
[1] CASTEM 2000, Guide dutilisation, CEA, France, 1990.
[2] Pegon, P., Model Implementation in CASTEM 2000: Some General Considerations and a Simple Practical Realization, JRC-Technical Note No. I.93.06, 1993, JRC-Special Publication No
I.94.05, 1994 and CEA-Report DMT/94-195, 1994.
[3] Oller, S., Barbat, A.H., Oate, E., Hanganu, A., A Damage model for the seismic analysis of
building structures, Proc. 10th World Conf. on Earthquake Engng., Ed. A.A. Balkema, Rotterdam, Vol. 5, pp 2593-2598, 1992.
[4] Combescure, D., Pegon, P., A Fiber model accounting for transverse shear in CASTEM 2000,
JRC-Special Publication, to appear in 1994.
[5] Hognestad, E., A Study of Combined Bending and Axial Load in Reinforced Concrete, Bulletin Series 339, Univ. of Illinois Exp. Sta., November 1951.
[6] Eurocode No 8 - Structures in Seismic Regions - Design, Part 1, General and Building, May
1988.
[7] Tassios, T. and Lefas, N., Ductility of confined reinforced concrete elements, Report NTUA (in
Greek), Athens, Greece, 1983.
[8] Park, R., Priestley, M., Ductility of Square-Confined Concrete Columns, Journal of the Structural Division, ASCE, Vol. 108, No ST4, pp. 929-950, April 1982.
[9] Mercer, C. and Martin, J., A Beam Element for Cyclically Loaded Reinforced Concrete Structures, Thecnical Report No. 98, FRD/UCT Centre for Research in Computational and Applied
Mechanics, University of Cape Town, November 1987.
[10] Barzegar-Jamshidi, F. and Schnobrich, W., Nonlinear Finite Element Analysis of Reinforced
Concrete Under Short Term Monotonic Loading, Civil Engrg. Studies, SRS 530, Univ. of Illinois at Urbana-Champaign, Urbana, Illinois, November 1986.
[11] Monti, G. and Nuti, C., Nonlinear Cyclic Behaviour of Reinforcing Bars Including Buckling,
Journal of Structural Engineering, Vol. 118, No 12, December 1992.
[12] Menegotto, M. and Pinto, P., Method of Analysis for Cyclically Loaded Reinforced Concrete
Plane Frames Including Changes in Geometry and Nonelastic Behaviour of Elements under
Combined Normal Force and Bending, IABSE Symposium on Resistance and Ultimate
Deformability of Structures Acted On by Well-Defined Repeated Loads, Final Report, Lisbon,
1973.
0
[13] Donea, J. and Lamain, L.G., A modified representation of transverse shear in C quadrilateral
plate elements, Int. J. Num. Meth. Engng., Vol. 63, pp 183-207, 1987.
[14] Pegon, P., A Timoshenko simple beam element in CASTEM 2000, JRC-Technical Note No.
I.93.05, 1993, JRC-Special Publication No I.94.04, 1994 and CEA-Report DMT/94-200, 1994.
[15] ESOPE/GEMAT, MANUEL dutilisation, CEA, France, 1989.
[16] Park, R. and Paulay, T., Reinforced Concrete Structures, John Wiley & Sons, 1975.
[17] Pipa, M., Carvalho, E., Cyclic Tests on R/C Cantilever T Beams with B500 Steel - Cooperative
Research Program on Seismic Response of R/C Structures, LNEC, Lisboa, June 1993.
[18] Bousias, S., Load Path Effects in Reinforced Concrete Column - Biaxial Bending with Axial
Force, Technical Note No I.92.63, Commission of the European Communities, Joint Research
Centre, Ispra Site, Italy, June 1992.
50 of 55
51 of 55
pc12 = 20.0d-2
0.0d-2;
carc5 = pc10 d 10 pc11 tran 1 pc12;
*
mesc = coul carc5 bleu;
*
elim (mesf et mesu et mesc) 0.001;
*-----------------------------------------*
TOTAL MESH (SECTION)
*-----------------------------------------*
titre LNEC TESTS - "T" SECTION;
*
trac (mesf et mesu et mesc);
*
*-----------------------------------------*
DEFINITION OF THE TESTED BEAM (BEAM)
*-----------------------------------------opti dime 3;
*
pp0 =
0.0d-2 0. 0.;
pp1 = 15.0d-2 0. 0.;
pp2 = 40.0d-2 0. 0.;
pp3 = 75.0d-2 0. 0.;
pp4 = 110.0d-2 0. 0.;
pp5 = 150.0d-2 0. 0.;
*
llb = pp0 d 1 pp1 d 1 pp2 d 1 pp3 d 1 pp4;
llp = pp4 d 1 pp5;
*-----------------------------------------*
DEFINITION OF THE AUXILIARY BEAM (FOR THE MAHO COMPONENT)
*-----------------------------------------pa0= 0 0 0; pa1=1 0 0; lla=pa0 d 1 pa1;
*-----------------------------------------*
CARACTERIZATION OF THE STEEL AND CONCRETE MODELS (SECTION LEVEL)
*-----------------------------------------*-----------------------------------------*
Quadrangular Elements QUAS
*-----------------------------------------modf=mode mesf mecanique elastique plastique acier_uni quas;
modu=mode mesu mecanique elastique plastique beton_uni quas;
modc=mode mesc mecanique elastique plastique beton_uni quas;
*-----------------------------------------*
Triangular Elements TRIS (commented out)
*-----------------------------------------* modf=mode mesf mecanique elastique plastique acier_uni tris;
* modu=mode mesu mecanique elastique plastique beton_uni tris;
* modc=mode mesc mecanique elastique plastique beton_uni tris;
*-----------------------------------------*
Steel
*-----------------------------------------matf=mate modf YOUN 2.06e5 NU
0.30
STSY 548.0 EPSU .098
STSU 625.0 EPSH 0.023
FALD 4.900 A6FA 620.0 CFAC 0.5 AFAC 0.008
ROFA 20.0
BFAC 0.0116 A1FA 18.5 A2FA 0.15;
carf=carb modf ALPY 1.0
ALPZ 1.0;
*-----------------------------------------*
Unconfined concrete
*-----------------------------------------matu=mate modu YOUN 0.35e5 NU
.25
STFC 35.0
EZER .002
STFT 0.0
ALF1 .22429 OME1 .12953
52 of 55
ZETA 100.0
ST85 .00
TRAF 2.0;
caru=carb modu ALPY .66
ALPZ .66;
*-----------------------------------------*
Confined concrete
*-----------------------------------------*
Initial concrete Young modulus = 2 * STIFC / ( BETA * EZERO )
*-----------------------------------------matc=mate modc YOUN 0.33e5 NU
.25
STFC 35.0
EZER .002
STFT 0.0
ALF1 .22429 OME1 .12953
ZETA 0.0
ST85 10.0
TRAF 2.0;
carc=carb modc ALPY .66
ALPZ .66;
*-----------------------------------------* Complete model and material at the SECTION level
*-----------------------------------------modq=modf et modu et modc;
macq=matf et matu et matc et carf et caru et carc;
*-----------------------------------------* Determination of the hook matrix at the BEAM level
*-----------------------------------------moa=mode lla mecanique elastique section timo;
maa=mate moa MODS modq MATS macq;
hoa=hook moa maa;
lrig=extr hoa MAHO 1 1 1;
*-----------------------------------------*
CHARACTERIZATION OF THE BEAM
*-----------------------------------------modfi=mode llb mecanique elastique section plastique section timo;
matfi=mate modfi MAHO lrig MODS modq MATS macq;
*
modpo=mode llp mecanique elastique section timo;
modpo=mode llp mecanique elastique section plastique section timo;
matpo=mate modpo MAHO lrig MODS modq MATS macq;
*
modto=modfi et modpo;
matto=matfi et matpo;
*-----------------------------------------*
STIFFNESS MATRIX
*-----------------------------------------rigi1 = rigi modto matto;
*
rigb1 = bloq DEPL pp0;
rigb2 = bloq ROTA pp0;
rigdz = bloq UZ pp5;
*
rigt = rigi1 et rigb1 et rigb2 et rigdz;
*-----------------------------------------*
LOADING (CYCLIC DISPLACEMENT ON THE TOP OF THE CANTILIVER BEAM)
*-----------------------------------------*
Loading program
*-----------------------------------------dhoz = prog -.0002 pas -.002 -10.0e-3
pas -.002 -30.0e-3
pas .002
30.0e-3
pas -.002 -40.0e-3;
npdy = dime dhoz;
xxxx = prog 1. pas 1. npas (npdy-1);
evolz = evol manu xxxx dhoz;
*-----------------------------------------*
Horizontal displacement on the top of the beam
*------------------------------------------
53 of 55
deslo = 1.0;
fdepz = depi rigdz deslo;
*
cargz = char fdepz evolz;
chart = cargz;
*-----------------------------------------*
RESOLUTION THE NON-LINEAR PROBLEM (USE OF NONLIN)
*-----------------------------------------tab1 = table;
tab1.PLASTIQUE = vrai;
tab1.MOVA = RIEN;
tab1.MAXITERATION = 20;
tab1.KTA = vrai;
*
nonlin rigt matto chart xxxx modto tab1;
*-----------------------------------------*
ANALYSIS OF THE RESULTS
*-----------------------------------------*
Curve (shear top force)/(horizontal top displacement)
*-----------------------------------------recont=tab1.RESUCONT;
redepl=tab1.RESUDEPL;
indj=index redepl;
ftop=prog 0;dtop=prog 0;
j=0; repe loopt (dime indj); j=j+1;
timj=indj.j;
ftop=inse ftop (extr (bsig modto recont.timj matto) pp5 FZ) (j+1);
dtop=inse dtop (extr redepl.timj pp5 UZ) (j+1);
fin loopt;
*
tt = table;
tt.1 = MARQ CARR;
tt.2 = MARQ CROI;
fdtop=evol manu DISPL. (m) dtop FORCE (MN) ftop;
*-----------------------------------------*
Reference curve (for checking the results)
*-----------------------------------------dtopref=prog
.00000000000000E+00 -2.00000000000000E-04 -2.16000000000000E-03
-4.12000000000000E-03 -6.08000000000000E-03 -8.04000000000001E-03
-1.00000000000000E-02 -1.20000000000000E-02 -1.40000000000000E-02
-1.60000000000000E-02 -1.80000000000000E-02 -2.00000000000000E-02
-2.20000000000000E-02 -2.40000000000000E-02 -2.60000000000000E-02
-2.80000000000000E-02 -3.00000000000000E-02 -2.80000000000000E-02
-2.60000000000000E-02 -2.40000000000000E-02 -2.20000000000000E-02;
dtopref=dtopref et (prog
-2.00000000000000E-02 -1.80000000000000E-02 -1.60000000000000E-02
-1.40000000000000E-02 -1.20000000000000E-02 -1.00000000000000E-02
-8.00000000000000E-03 -6.00000000000000E-03 -4.00000000000000E-03
-2.00000000000000E-03 6.07153216591883E-18 2.00000000000001E-03
4.00000000000001E-03 6.00000000000001E-03 8.00000000000001E-03
1.00000000000000E-02 1.20000000000000E-02 1.40000000000000E-02
1.60000000000000E-02 1.80000000000000E-02 2.00000000000000E-02);
dtopref=dtopref et (prog
2.20000000000000E-02 2.40000000000000E-02 2.60000000000000E-02
2.80000000000000E-02 3.00000000000000E-02 2.80000000000000E-02
2.60000000000000E-02 2.40000000000000E-02 2.20000000000000E-02
2.00000000000000E-02 1.80000000000000E-02 1.60000000000000E-02
1.40000000000000E-02 1.20000000000000E-02 1.00000000000000E-02
8.00000000000000E-03 6.00000000000000E-03 4.00000000000000E-03
54 of 55
55 of 55