Fibre Beam Element

A FiBRE/TIMOSHENKO 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

Special Publication Nr. I.94.31 July 1994

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

A FiBRE/TIMOSHENKO BEAM ELEMENT IN CASTEM 2000

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

6.5 Operators working for the ‘SECTION’ formulation 38 6.5.1 Biaxial bending of a section (‘MOCU’ operator) 38 6.6 How to implement a new fibre model 39 7 NUMERICAL APPLICATIONS 41 7.1 Introduction 41 7.2 Monotonic tests 41 7.2.1 Doubly reinforced rectangular cross sections 41 7.2.2 Influence of the concrete tensile strength 42 7.3 Cyclic tests 43 7.3.1 Uniaxial Loading on a Cantilever T Beam 43 7.3.2 Biaxial Loading on a R/C Square Section 45 8 SUMMARY AND CONCLUSION 49 ACKNOWLEDGEMENTS 49 REFERENCES 50 ANNEX 1- A SAMPLE GIBIANE INPUT FOR FIBRE APPLICATION 51


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

Figure 1- The scope of the fibre modelling


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

2. THE BASIC ASSUMPTIONS OF TIMOSHENKO BEAM THEORY AND INTRODUCTION OF THE FIBRE MODELLING 2.1 Assumed displacement field Consider the beam of Fig. 2. It is assumed that planes which are normal to the axis Ox before deformation remain plane, but not necessarily normal to the Ox axis after deformation. This implies that the displacement field u = ( u x, u y, u z ) can be expressed in terms of the generalized displacement of any section, that is the displacement U = ( Ux, U y, U z ) of the Ox axis and the rotations Θ = ( Θ x, Θ y, Θ z ) of the plane. U ( P ) = U ( R ) + Θ ∧ RP or, using components: u x ( x, y, z ) = U x ( x ) – yΘ z ( x ) + zΘ y ( x ) u y ( x, y, z ) = U y ( x ) – zΘ x ( x )

(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)

Let U denote the ‘vector’ of generalized displacement U = ( U, Θ )

(3)

y Θz Θx

O

z

Θy

R

x

P

Figure 2- Timoshenko beam


7 of 55


2.2 Strain-displacement relationships Assuming small deflection theory, the axial and shear components of the strain tensor are obtained from (1)-(2) as follows: εx = γxy = γxz =

dΘ y dΘz ∂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.

2.3 Virtual work expression Consider a Timoshenko beam of length L , of reference axis Ox and section S ( x ) . The beam is subjected to a lateral distributed loading q of components qy and q z . If the beam undergoes a set of virtual displacements δu x , δu y and δu z and rotations δΘ x , δΘ y and δΘ z , the virtual work expression is: L

L

∫ ∫ ( δεx σx + δγxy τxy + δγxz τxz )dS dx – ∫ ( δUy qy + δUz qz )dx 0S

= 0

(5)

0

where the first term takes the following form using (2) and (4): L

 d



∫ ∫   d x δUx + z d x δΘ y – y d x δΘz σx + δβy τxy + δβz τxz + d x δΘx ( – zτxy + yτxz ) dS dx d

d

d

0S

Introducing now the following quantities: • normal force: Fx =

∫ σ x dS

(6.1)

S

• shear forces: Fy =

∫ τxy dS S

Fz =

∫ τxz dS

(6.2)

S

• twisting moment: Mx =

∫ ( – zτxy + yτxz )dS

(6.3)

S


8 of 55


• bending moments: My =

∫ zσx dS

M z = – ∫ yσ x dS

S

(6.4)

S

the virtual work takes the expression: L

∫  Fx d x δUx + Fy δβy + Fz δβz + Mx d x δΘx + My d x δΘy + Mz d x δΘz dx d

d

d

d

0

(7)

L

– ∫ ( δU y q y + δU z q z )dx = 0 0

It is thus possible to associate to the components of the ‘generalized’ stress vector F = ( F x, F y, F z, M x, M y, M z ) the components of the ‘generalized’ strain field E defined by dU x dΘ x dΘ y dΘ z E = ( E x, E y, E z, C x, C y, C z ) =  , β y, β z, , , dx dx dx dx 

(8)

2.4 The elastic case For an isotropic elastic material for which E is the Young modulus and ν the Poisson ratio, the constitutive law takes the simple form: σx σy σz

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

The stress-strain relationships (9.1)-(9.2) then simplify under the form: σ x = Eε x

τ xy = Gγ xy

τ xz = Gγ xz

(10)

where G = E ⁄ 2 ( 1 + ν ) is the shear modulus Assume now that these material properties are not homogeneous on the section S and denote by E s and G s the elastic properties of any ‘fibre’ of the section. Using (10) and (2)-(4), the components of


9 of 55


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

S

∫ Es zdS ∫ Es z

My =

S

Mz

2

S

∫ Es y

S

0

Cz dS (11)

– ∫ G s zdS S

0

=

2

S

S

Fy

Ex

dS – ∫ Es yzdS ⋅ C y

– ∫ E s ydS – ∫ E s yzdS

∫ Gs dS Mx

S

S

S

Fz

– ∫ E s ydS

∫ Gs dS

∫ Gs ydS

S

S

– ∫ Gs zdS ∫ G s ydS ∫ G s ( y + z )dS 2

2

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.

2.5 Generalized fibre modelling The elastic constitutive relation (9.1)-(9.2) may be considered as the reversible part of a more complex constitutive relation of the generic form σ = f ( ε, ... ) In particular, any constitutive law developed for performing refined analysis (e.g. plasticity or damage modelling) can be used at any point of the section. Though it may be questionable to use such a kind of 3-D constitutive modelling in a structural beam element (and not in a 3-D continuum element), a generalized fibre modelling is sometimes used when very complex interactions between normal force, bending and transverse shear are expected (see for instance [3]). The Mazars model is currently used in CASTEM 2000 in order to evaluate the potential of this approach [4].

2.6 Usual fibre modelling The usual fibre modelling, which accurately accounts only for interaction between normal force and bending, uses (9.1)-(9.2) as the elastic part of a model in which a complex relationship is introduced only between the axial component σ x of the stress field (the fibre stress) and the axial deformation ε x (the fibre strain) σ x = f ( ε x, ... )

(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


2.7 Cross-sectional warping It should be noted that, in order to account for possible cross-sectional warping, the expressions (6.2)-(6.3) of the shear forces and the twisting moment are changed according to Fy =

∫ αs τxy dS y

∫ αs τxz dS z

Fz =

S

Mx = y

S

∫ (–

y α s zτ xy

(13)

z α s yτ xz )dS

+

S

z

where α s and α s are parameters lower or equal to 1. In elasticity, the expression (11) of the shear forces and the twisting moment is then changed to:

∫ αs Gs dS y

Mx

y

S

Fy Fz

– ∫ α s G s zdS

0

S

0

=

∫

∫

z α s G s dS

S

–∫ S

y α s G s zdS

∫

z α s G s ydS

S


z α s G s ydS

S

∫

z 2 Gs ( αs y

+

y 2 α s z )dS

Ey ⋅ Ez

(14)

Cx

S

11 of 55



12 of 55

SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING

3. SOME TYPICAL CONSTITUTIVE RELATIONS FOR USUAL REINFORCED CONCRETE FIBRE MODELLING 3.1 Concrete Constitutive Laws 3.1.1 Monotonic loading 3.1.1.1 Compression stresses As the concrete is subjected to loading of increasing intensity, it undergoes different phases of damage, from micro-cracking up to ultimate failure. Fig. 3 shows a typical response curve of the concrete under increasing deformation. After an almost linear zone up to about one-half of the compression strength, the concrete presents a degradation of the stiffness due to micro-cracking and consequently nonlinear behaviour is pronounced. The peak stress zone in the stress-strain curve, where transversal strains increase rapidly owing to internal cracking parallel to the loading direction, is relatively sharp for high-strength concrete and has a flat top for low-strength concrete. The strain at the maximum stress is approximately ( ε = 0.2 %). After the peak, the concrete can still carry on compression stresses although it presents high decreasing strength ratios for increasing strain. The evolution of the post-peak curve depends mainly on the degree of confinement (e.g. presence of stirrups on R/C structural elements). This phenomenon is explained by the triaxial behaviour of the concrete: due to the Poisson ratio, during an uniaxial compression loading the transversal strains increase. Although the volume of the specimen starts to decrease, at stresses near the compression strength, the transversal strains become so high that the volume of the specimen will actually increase. If transversal bars exist, the expansion of the concrete core is prevented, subjecting the material to lateral compression forces that improve the ductile performances of the concrete. In this case, the post-peak curve presents not only a lower decreasing strength ratio but also a residual strength for important strains, reflecting the coupling between the confinement effect and the interlocking of the aggregate. Analytical modeling of concrete behaviour can be performed using several approaches ranging from the classical Plasticity Models to the recent Continuum Damage Mechanics Models. However, due to the complex behaviour under cyclic loading, a more flexible approach has been adopted herein. Thus, a two branches law of Hognestad type [5] was considered for the non-linear behaviour of the concrete under monotonic compression loading. The first branch, a parabolic function, defines the ascending part of the curve and goes from zero to the maximum compression stress point, (ε co, σ co). The second branch, a descending straight line, represents the concrete softening behaviour after maximum strength until failure: ε σ ε -------- = -------  2.0 – ------- ε co σ co ε co

(15)

σ -------- = 1.0 + Z ( ε – ε co ) σ co

(16)

for 0 < ε < ε co ,

for ε ≥ ε co till failure. The slope of the second branch, Z , depends on the degree of confinement of the concrete: the stirrups, preventing the normal expansion of the concrete core, determine a triaxial compression state in


13 of 55


σ εco∗ εco

ε

Unconfined (Z)

σco

Confined (Z*)

σco*

Figure 3- Confinement effect in the concrete.

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


σ Tension law

σt

σ εco* ε

εt

εtm

ε

σpt σco*

Compression law

Figure 4- Concrete response curves for monotonic loading.

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


for 0 < ε < ε t , and r – ( ε ⁄ εt ) σ = σ t  ----------------------- r–1

(20)

r = ε tm ⁄ ε t

(21)

for ε t ≤ ε < ε tm , with

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

Figure 5- Concrete response under compression cyclic loading.


16 of 55


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

6

1

4

ε

8

σ = f(ε)

2

5

Figure 6- Numerical model for the concrete behaviour under cyclic loading.


17 of 55


σ Envelope

σ = f(ε) 2

εpl

σt

3 4

6

εco

εtr1

εt

1 7

εtm

5

εtr2 εt

ε

εtm

Figure 7- Numerical model for the concrete under tensile stresses.

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.

3.2 Steel constitutive law 3.2.1 Monotonic loading Typical monotonic stress-strain curves for steel in general come from bars loaded in tension. After a linear elastic zone, the steel exhibits a yielding plateau, i.e., a point beyond which the stress is almost constant. In the plateau zone the material suffers important and irreversible changes in its internal matrix, and the Poisson ratio, generally in the range 0.2 to 0.3 , becomes equal to 0.5 , meaning that during this stage the steel bars do not change their volume. The stress defining the plateau, referred to as the yield stress, may present two different values: the upper and the lower yield stress, the last one being the value usually adopted to represent it. The existence of the upper value is determined not only by the geometric characteristics of the bar but also by the speed of the test. The length of the plateau is mainly determined by the strength and internal characteristics of the steel: steel bars with higher strength and carbon levels, usually present a shorter plateau. After this stage, the steel presents a strain-hardening range up to the maximum strength, followed by a strain-softening curve pointing to failure. Either by using a cold worked or a high-strength and carbon steel, the tensile monotonic curve may not present the yielding plateau, going from the linear elastic zone immediately to the strain-hardening curve. Although the behaviour of steel bars under compression forces can be considered identical to the tensile one, if the tie spacing is not small enough steel reinforcing bars in R/C elements can evidence inelastic buckling. Experimental tests [11] show that when the ratio between the length and the diameter of the bar exceeds 5 , this phenomenon occurs. In this case, the yielding plateau and the strainhardening zones may even disappear from the response, being substituted by a softening behaviour curve immediately after the linear elastic zone.


18 of 55


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


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

Figure 8- Numerical model for the steel under cyclic loading.


20 of 55


3.2.2.2 Inelastic buckling The proposed model also accounts for the inelastic buckling effect. Tests made by Monti and Nuti [11] on several steel bars showed that, for ratios between the length L and the diameter D less or equal to 5 , the compression curve is similar to the tensile one, i.e., no buckling effect is observed. According to the same author, if L ⁄ D ratios are greater than 5 , modifications should be introduced in the model, namely the unloading curve from tensile stresses and the reloading curve after reversal from compression (see Fig. 9). In this case the post-yield compression zone has a softening behaviour and the new b factor is computed by:

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 L⁄D 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

εo

ε5

ε

γs bcxE

bxE Er

σs∗

Figure 9- Inelastic buckling on the steel cyclic behaviour.


21 of 55


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

FINITE ELEMENT IDEALIZATION FOR TIMOSHENKO BEAMS

4. FINITE ELEMENT IDEALIZATION FOR TIMOSHENKO BEAMS 4.1 Displacement and strain representation In the linear beam element, the displacement field U = ( Ux, U y, U z ) is approximated by the relationship: 1

1

2

U ( x ) = N ( x )U + N ( x )U 1

2

2

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)

and l = x 2 – x 1 is the length of the element. Similarly, the rotations Θ = ( Θ x, Θ y, Θ z ) are represented by: 1

1

2

Θ ( x ) = N ( x )Θ + N ( x )Θ 1

2

(35)

2

where Θ and Θ are the nodal rotations. 1

1

2

2

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


4.2 Stiffness matrix evaluation In the elastic case, the introduction in the virtual work principle (7) of the expressions (11) and (14) relating the generalized stresses and the generalized strains (8) leads to:  L    δE K Edx – – ( δU q + δU q )dx = 0 y y z z ∫  ∫   0  0

L

T

(38)

where the section constitutive matrix K is given by: K 11 =

∫ Es dS

K 15 =

S

K 22 =

∫ Es zdS

K 16 = – ∫ E s ydS

S

∫ αs Gs dS

S

K 24 = – ∫ α s G s zdS

y

y

S

K 33 =

S

∫ αs Gs dS z

K 34 =

S

K 42 = – ∫

y αs G s zdS

∫ αs Gs ydS z

S

K 43 =

S

∫

z α s G s ydS

S

K 51 =

∫ Es zdS

(39)

K 44 =

∫

z 2 ES ( αs y

S

K 55 =

S

∫ Es z

2

dS

K 56 = – ∫ E s yzdS

S

K61 = – ∫ E s ydS

+

y 2 αs z )dS

S

K 65 = – ∫ E s yzdS

K 66 =

∫ Es y

2

dS

the other elements being null. Using now a discretization of the domain [ 0, L ] in elements and introducing the linear approximation of the type (33), (34) and (35), (38) is equivalent to the solution of a linear system. KΦ = f

(40)

e

The contribution K of the element to the global stiffness K takes as usual the following form: x2 e

K =

∫ [B]

T

K B dx

(41)

x1

which can be evaluated exactly using a one-point Gauss-Legendre rule. Integrating now the second term in (38), the consistent nodal force vector for the element is derived: qy l qz l qy l qz l e f =  0, -------, -------, 0, 0, 0, 0, -------, -------, 0, 0, 0  2 2  2 2

(42)

which, unlike the Euler-Bernouilli cubic Hermitian element, only has lateral nodal point forces and no nodal moments.

4.3 Element stress resultant and nodal resultant Using the expressions (36) of the generalized strain and, in the elastic case, the relations (11) and (14), it is easy to observe that the axial force F x and the moments M x , M y and M z are constant. The


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

2

F ξ = – Fξ

ξ = x, y, z

Fξ = Fξ

1

2

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

4.4 The mass matrix Since the Timoshenko element is to be used also for performing dynamic analysis, it may be interesting to derive its element mass matrix. The virtual work W I of the inertial forces which have to be added to (5) takes the form: 2

L

2

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)

0

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 =

2

∫ δU

T

M

dU dt

0

2

(45)

dx

where the section mass matrix M is given by: M 11 =

∫ ρs dS

M15 =

S

M 22 =

∫ ρs zdS S

∫ αs ρs dS

S

M 24 = – ∫ α s ρ s zdS

y

y

S

M 33 =

S

∫ αs ρs dS z

M 34 =

S

M 42 = – ∫

y αs ρ s zdS

S

M 51 =

M 43 =

∫

z α s ρ s ydS

M 55 =

S

M 61 = – ∫ ρ s ydS

∫ αs ρs ydS z

S

S

∫ ρs zdS

M 16 = – ∫ ρ s ydS

M 44 =

∫

z 2 ρS ( αs y

+

S

∫ ρs z

2

dS

M 56 = – ∫ ρ s yzdS

S

M 65 = – ∫ ρ s yzdS


(46) y 2 α s z )dS

S

M 66 =

∫ ρs y

2

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

1

1

2

2

x1

1

2 l ⋅ N )dx = --6

e

the element mass matrix M may be readily derived.


26 of 55

THE TIMO ELEMENT WITH FIBRES IN CASTEM 2000

5. THE TIMO ELEMENT WITH FIBRES IN CASTEM 2000 5.1 Overview The implementation of the TIMOshenko element with fibres could be realized following a 2-D multi-layer approach, for which selected integration points are used for performing the integration over the section (or part of the section defined by a material). This approach would generalize the one which has been introduced in CASTEM 2000 for treating thick shells. Though being natural at a first glance, this approach suffers from severe drawbacks: • the difficulty to account for the a priory unknown complexity of the section that have to be described, • the resulting cryptic character of any solution aiming to have some generality, • the difficulty to post-treat the relevant informations at the level of the section. In order to overcome these difficulties, a rather different approach is adopted here. In fact, the sections are described by a model with various sub-zones corresponding to various materials (say steel fibre, concrete fibre,...) and associated with the meshes of these various components. This description allows in turn to define the beam modelling. The description of the sections (shape, material characteristics) is then general, easy to generate and to post-treat. The evaluation of the stress resultant for each beam element proceeds as follow (see illustration in Fig. 10):

( U, Θ ) ⇒ E

F ⇒ ( F, M )

⇒

⇒

‘BEAM’ level

ε⇒σ

‘SECTION’ level

Figure 10- The fibre element


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

1

2

2

• Computation of the stress resultant ( F , M , F , M ) for the beam element.

5.2 Short remarks about CASTEM 2000 In order to understand better the next sections, it may be interesting to present an outline of the modelling concepts of CASTEM 2000. In-depth descriptions may be found in [1] or [2]. CASTEM 2000 uses an internal language (GIBIANE) which allows the user to create, manipulate and destroy objects belonging to a relevant class by means of operators. The implementation of the object classes is realized by means of generic data-structure types, the implementation of the operators by means of subroutines. The geometry of a problem is specified creating objects of the type ‘mesh’. The corresponding datastructure is subdivided in zones corresponding to an identical parent geometric element. A geometric element is different from a finite element since the same geometry may be associated for example withcontinuum or structural elements. It is said that the geometrical element does not have the same formulation. The modelling of a mechanical problem is specified creating objects of the type ‘model’. The corresponding data-structure is subdivided in zones corresponding to identical modelling and identical parent finite elements. They are called the sub-zones of the model. CASTEM 2000 being a displacement based finite-element code, the quantities important for the modelling (say material characteristics, strain, stress, internal variables, etc...) have to be known at the integration points of the mesh. These informations are hold by objects of the type ‘field-by-element’. These objects have specific types (‘characteristics’, ‘strain’, ‘stress’, ‘internal variable’, etc...) depending on their content. Since almost all the operators creating an object of the type ‘field-by-element’ involve a model, a field-by-element is subdivided in sub-zones according to a model. The operations involving a model and the related field-by-elements are generally implemented by means of 3 nested loops: the first loop on the sub-zone of the specified model, the second on the elements of the zone and the third on the integration points of the element. In terms of the data-structure, a field-by-element on each sub-zone is a collection of two-entry arrays containing the values of the field, associated to each components of the field. The two entries are the number of elements in the sub-zone and the number of integration points. The values are either real numbers or more complex typed data-structures. The number and the name of the components depends on the specific type of the field-by-element, on the name of the model and on the formulation of the finite element. For example, consider a field-by-element of the specific type ‘stress’ in a sub-zone of parent finite element ‘QUA8’ (continuum 2-D plane strain element corresponding to a quadratic 8-node quadrilateral geometric element with 3x3 Gauss-integration points). It is a collection of two-entry arrays corresponding to the active components σ xx (component ‘SMXX’), σ yy (component ‘SMYY’), σ zz (component ‘SMZZ’) and σ xy (component ‘SMXY’) of the stress


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.

5.3 Global and local referentials for structural beam elements In order to be able to compute the equilibrium of a frame structure, CASTEM 2000, as many codes do, distinguishes between a global referential ( OXYZ ) in which the equilibrium (or the stiffness or the mass) is assembled and a local referential ( Oxyz ) in which the computations are mainly performed. A local referential is associated to each element in the following way: • O position of the first node of the element; • Ox axis passing through the two nodes of the element and oriented from the first to the second node. • Oy axis orthogonal to Ox , in the plane Oxv where v is a vector supplied by the user (default OY ) and verifying Oy ⋅ v > 0 . • Oz axis defined by Oz = Ox ∧ Oy , for which ( Oxyz ) is a direct orthogonal referential. Then for all computations at element level, CASTEM 2000 computes the local referential ( Oxyz ) , eventually projects global nodal information (say the generalized displacement field ‘UX’, ‘UY’, ‘UZ’, ‘RX’, ‘RY’, ‘RZ’ for U x , U y , U z , Θ x , Θ y , Θ z ) expressed in ( OXYZ ) ) in the local referential, performs element computations in the local referential and eventually projects back the results in the global one before assembling. The same approach is of course used for the ‘TIMO’ element with fibres so that many functions used in the implementation of the usual Timoshenko element (‘TIMO’) or the Hermitian element (‘POUT’) can be reused.

5.4 The two-level approach As outlined in the previous overview section, a two-level approach which uses the existing data structures and operators of CASTEM 2000 in order to specify the behaviour of each individual section and in order to post-treat them has been adopted: the first level is the description of the section, the second is the description of the beam. 5.4.1 Description of the section. Each section is described by a object of the type ‘model’, associated to mesh(es) of simple 2-D linear elements (the 3-node triangle ‘TRIS’ and the 4-node quadrilateral ‘QUAS’) associated with a new formulation called ‘SECTION’. The local referential yOz of the section is assumed to coincide with the global referential XOY of the mesh. This assumption allows to generate the mesh of the various sections using the 2-D options of CASTEM 2000, and to switch to the 3-D level only when the beam


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


case), INATU=39 (‘BETON_UNI’) or INATU=40 (‘ACIER_UNI’); at the ‘beam’ level CMATE=’SECTION’ and INATU=0 (elastic case) or INATU=41 (‘SECTION’). Finally, to perform the computations at the ‘SECTION’ level, it is interesting to define a relative index NUFIBR ordering the fibre material. A new subroutine TEMANF associates to the absolute index INATU the relative index NUFIBR: NUFIBR=1 for ‘BETON_UNI’, NUFIBR=2 for ‘ACIER_UNI’ and NUFIBR=0 in the other cases and in particular the elastic one.

6.3 Identification subroutines As usual, when a new formulation is introduced, the number, the name and the type of values of the components of the field by elements that can be associated with this formulation should be declared in the related identification routines. For the ‘SECTION’ formulation, three specific field-by-elements are required for operational reasons: ‘CHARACTERISTICS’ for the material and geometrical characteristics, and ‘CONTRAINTE’ (‘stress’) and ‘VARIABLES INTERNES’ (‘internal variables’) for the non-linear sections. The corresponding modifications/additions have been performed: • Material characteristics: the usual elastic material characteristics (‘YOUN’= E , ‘NU’= ν , ‘RHO’= ρ ) are specified by IDMATR, whereas the non-linear material characteristics of the two available ‘plastic’ models (‘ACIER_UNI’ and ‘BETON_UNI’) are specified as usual in the satellite of IDMATR called IDPLAS. The name of the non linear components for the ‘BETON_UNI fibres are (see Section 3.1): ‘STFC’ (maximum compression stress σ co ), ‘EZER’ (strain at the maximum compression stress ε co ), ‘STFT’ (maximum tensile stress σ t ), ‘ALF1’ ( α factor), ‘OME1’ ( ω ω factor), ‘ZETA’ (modulus Z of the descending part of the compression curve), ‘ST85’ (stress value σ pt defining the compression plateau), ‘TRAF’ (descending traction curve factor r = ε tm ⁄ ε t ). The name of the non linear components for the ‘ACIER_UNI’ fibres are (see Section 3.2): ‘STSY’ (yield stress σ sy ), ‘STSU’ (ultimate stress σ su ), ‘EPSH’ (hardening strain ε sh ), ‘EPSU’ (ultimate strain ε su ), ‘ROFA’ ( R o factor), ‘BFAC’ ( b factor), ‘A1FA’ ( a 1 factor), ‘A2FA’ ( a 2 factor), ‘FALD’ ( L ⁄ D factor), ‘A6FA’ ( a 6 factor), ‘CFAC’ ( c factor), ‘AFAC’ ( a factor). y

z

• 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


The manual of the operator ‘MATE’ has been changed according to the new possible characteristics.

6.4 Complex operators working for the ‘beam’ (‘POUTRE’) formulation The ‘simple’ operators working for the ‘beam’ (‘POUTRE’) formulation and the ‘TIMO’ element with fibres are those ones which do not use the model of the section. These operators does not have to be modified when fibre modelling is introduced. This is for example the case for the ‘EPSI’ operator (which computes the field-by-element of strain) and the operator ‘BSIG’ (which computes the ‘field-by-node’ of internal nodal force). Of course, all the operators using the information at the ‘SECTION’ level should be modified: this is the case of the operators ‘HOOK’, ‘SIGI’, ‘RIGI’ and ‘MASS’ used for performing elastic and/or dynamic analysis and of the operators ‘ZERO’, ‘ECOU’ and ‘KTAN’ used for performing non-linear analysis. All these ‘complex’ operators involve the particular field-by-element defined at the section level, introduced for performing the non-linear ‘SECTION’ computations. 6.4.1 Hook matrix calculation (‘HOOK’ operator) The ‘HOOK’ operator computes the field-by-element of (elastic) Hook matrix at the Gauss points of the mesh associated with the sub-zones of a given model with associated material characteristics.The Hook matrix which is constant for each element having the same material and geometric characteristic, is, for the TIMO element with fibres, the 6x6 section constitutive matrix K whose expression is given by (39). Due to the particular role played by this operator for building-up efficient elastic calculations (see Section 5.4.2) it is assumed that the component ‘MAHO’ of the ‘field-by-element’ of material characteristics is not (yet) existing. As a consequence the ‘HOOK’ operator works at both ‘beam’ and ‘SECTION’ level. As usual in this type of computation, a driver routine (here HOOK2P) performs a loop on each subzone of the input model, collecting and checking all the useful informations for each sub-zone and then, inside the loop, calls the routine (here HOOK2D) which achieves the work at the element level. HOOK2P has been modified in order to allow the ‘TIMO’ element to have ‘complex’ material characteristics (IF-THEN-ELSE structure on the value of CMATE). The subroutine HOOK2D has been modified using the same structure, and, when required, calls the new routine FRIGIE which returns two arrays of 12 integrals allowing to compute the Hook matrix and the section mass matrix. In fact, FRIGIE is a driver routine at the ‘SECTION’ level, which performs a loop on each sub-zone of the section model and which calls a new routine FRIGI2.FRIGI2 loops on the elements/Gauss points of the sub-zone, and computes the numerical contribution to the integrals which appear in the expressions of K (39) and M (46). When available at the ‘beam’ level, the array relative to the Hook matrix is elaborated (still in HOOK2D) by a new routine: DOHTIF. 6.4.2 Element elastic stress calculation (‘SIGM’ operators) Three possibilities exist in CASTEM 2000 for evaluating the field-by-element of (local generalized) stress at the Gauss point of a given mesh. These possibilities have to be compatible. • Compute the field-by-element of the (elastic) Hook matrix (‘HOOK’ operator) at the Gauss points of the mesh, compute the deformation using the ‘EPSI’ operator and make the product of the two results. • Use the operator ‘SIGM’ with the field-by-element of Hook matrix and the field-by-node of (generalized) displacement.


35 of 55


• 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


6.4.5 Initialisation (‘ZERO’ operator) The ‘ZERO’ operator is used, between other things, for making the initialization of the field-by-element of internal variables at the beginning of any non linear calculation. In particular, this is the case of the two components (‘CONS’ and ‘VAIS’) associated to the ‘beam’ formulation. As already discussed, these components are field-by-elements over the section. Instead of performing the initialization at the two levels, ‘beam’ and ‘SECTION’, which would require the knowledge of the ‘SECTION’ model defined as ‘beam’ material characteristics, only the ‘beam’ level is initialized according to the following trick: the two component fields are created, are correctly typed, but are associated to a “void” object. The further modelling routines, used for the non-linear ‘TIMO’ element with fibres, are allowed to recognize a “void” object as being a “null” object and to perform, if required, the initialization at the ‘SECTION’ level (calling usual slave routines of the operator ‘ZERO’). The trick used for the initialization at the ‘beam’ level is introduced in the routine ZEROP. 6.4.6 Non linear flow evaluation (‘ECOU’ operator) The ‘ECOU’ operator is used for performing all the incremental evaluation of the non-linear constitutive laws in CASTEM 2000 written under the form: σ

new

= σ

old

+ ∆σ ( σ

old

,η

old

e

, ∆σ )

new

(47) old

where the stress σ , at the end of the increment, is a function of the stress σ and of the internal old e variables η at the beginning of the increment, and of the increment of elastic stress ∆σ . The two levels of modelling are treated in a separate way. • The ‘beam’ (‘POUTRE’) level is implemented as usual [2]: introduction of a new model number (NOMATE), modification of ECOUL1 for accounting to the special value type of the non linear components of the material characteristics and of the internal variables (the default being real number value), addition of a slot in ECOUL2 and call to the new routine which is devoted to the fibre modelling: BIFLEX. • The subroutine BIFLEX, called for each integration point of each ‘TIMO’ element with fibres, is the interface between the two levels. In particular, the data coming from the ‘beam’ level are rearold ranged in such a way the model, the characteristics, the initial stress σ and the initial internal old variable η are accessible as usual. The behaviour of each fibre being specified kinematically (Eq. (4)), the increment of elastic generalized stress is inverted using the section stiffness. In consequence, the ‘SECTION’ models are generally formulated as σ

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


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.

6.5 Operators working for the ‘SECTION’ formulation Letting aside the operators ‘MODE’, ‘MATE’ and ‘CARA’ which work directly at the section level after having modified the model-related routines (see Section 6.2) and the characteristics-related identification routines (see Section 6.3), a new operator named ‘MOCU’ have been specifically developed for dealing with the section modelling. 6.5.1 Biaxial bending of a section (‘MOCU’ operator) The ‘MOCU’ operator (see chapter 5.4.1) allows to derive the bi-axial moment/curvature curves of any section described by a ‘SECTION’ model and the associated characteristics. Normal forces can A FiBRE/TIMOSHENKO BEAM ELEMENT IN CASTEM 2000

38 of 55


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.

6.6 How to implement a new fibre model The method to follow in order to implement a new ‘SECTION’ non linear model is more or less the one indicated in [2], with some slight differences due to the two-level implementation. A short guideline is here indicated. • Model: in the subroutine MODPLA add 1 to the number NMOD of possible models and add a new position to the array MOMODL by setting the name of the model. The index position of the model is the relative plastic model index INPLAS. • Model index: an absolute model index (INATU) is associated to the new model in the subroutine NOMATE. The correspondence between this absolute number and the relative ‘SECTION’ model index INFIBR should be established in the subroutine TEMANF. • Material characteristics: The names of the required non linear material characteristics are specified in the subroutine IDMATR by adding a new slot (ELSEIF structure for the local index IPLAC) corresponding to the new relative plastic model index. The names are stored in the array MOTS. • Internal variables: The names of the required internal variables are specified in the subroutine IDVARI by adding a new slot (ELSEIF structure for the local index MATPLA) corresponding to the new relative plastic model index. The names are stored in the array LESOBL. A position for tang the component ‘TANG’ (containing E ) should be presented. A position for the component ‘EPSO’ (containing the total elongation ε x of each fibre) is strongly advised for performing the post processing. • Interface with FCOUL2 (‘ECOU’): In FCOUL2 a new slot (ELSEIF structure for the relative fibre index INFIBR) should be added. A call to a subroutine, say FIXXXX, is then required. This routine is used as an interface between FCOUL2 and, say MOXXXX, which implements effectively tang the new model. It is also the role of the interface to collect the value of E and possibly the one of ε x . • Tangent modulus in FRIGT2 (‘KTAN’): in FRIGT2 a new slot (ELSEIF structure for the local index INFIBR) should be added. The local variable YOUNGT is set according to the value of the internal variable associated to the component ‘TANG’. This value is found in the array XVAR at the position decided by the developer.


39 of 55



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 LNEC’s 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.

7.2 Monotonic tests 7.2.1 Doubly reinforced rectangular cross sections The first numerical test is performed on a R/C rectangular section. Different volumetric steel ratios are considered at both sides of the beam. The loading, an increasing curvature path, is imposed to the structure and moments are computed by the fibre algorithm. The response curves (see Fig. 11) are to be compared with the theoretical results obtained by another algorithm but supposing the same behaviour (stress-strain laws) for the materials. The monotonic curve for the concrete considers the compression strength ( σ co = – 27.6MPa ) , for ( ε co = – 0.2%). The softening behaviour is described by ( Z = – 100.0 ). No compression plateau nor tensile stresses are considered. For the steel, an elastic perfectly plastic law is adopted. This two parameter law is defined by a yield stress ( σ sy = 276.0MPa ) and a Young modulus ( E = 207.5GPa ). Fig. 11 shows the fibre discretization considered in the section as well as the moment-curvature results for the different reinforcement ratios of Table 1. A very good agreement with the theoretical results presented by Park and Paulay [16] was achieved.


41 of 55


Table 1: Steel reinforcement ratios (see Fig. 11) Beam

1

2

3

4

5

6

7

ρ

0.0375

0.0375

0.0375

0.0250

0.0250

0.0125

0.0125

ρ’

0.0250

0.0125

0

0.0125

0

0.0125

0

Moment/(b*d2) (MN/m2) 11.0 A’s

d/10 1 2

8.8 3

d

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

A’s ρ’ = ------bd

2.80

X1.E-2

Figure 11- Fibre results according to Park and Paulay tests.

7.2.2 Influence of the concrete tensile strength The analysis of the influence of the concrete tensile strength in the global response of R/C structural elements under flexural moments is presented. The structure under analysis, a square doubly symmetrical section, is illustrated in Fig. 12 together with the fibre discretization. In this test, a compressive strength ( σ co = – 30.0MPa ) , for ( ε co = – 0.2 %) and a tensile strength ( σ t = 3.0MPa ) , were adopted for the concrete. For the steel, a Young modulus ( E o = 203.0GPa ) and a yield stress ( σ sy = 460.0MPa ) were considered. The concrete cover, ( c = 0.19m ), was modeled by two external fibres, one on the top and the other on the bottom of the section. For all the other fibres representing the concrete core, a confinement factor ( β = 1.331 ) was adopted. Since the aim of this test is the analysis of the response curve just before and after the cracking process, all the other parameters defining both the steel and the concrete constitutive laws are irrelevant to the response and so are not presented. Four numerical tests were performed and the associated results are presented in Fig. 12. In each one a different cracking model was adopted. For the first test a sudden rupture was chosen, ( ε tm = ε t ). In the second, third and fourth tests the ratio ( r = ε tm ⁄ ε t ) takes values equal to 1.5 , 2.0 , and 22.7 , respectively. The dashed line illustrates the response of the section when zero tensile strength is considered. The fourth result follows the suggestion of Barzegar [10].


42 of 55


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

3φ16mm

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.

7.3 Cyclic tests 7.3.1 Uniaxial Loading on a Cantilever T Beam A series of tests on R/C cantilever T beams were performed at LNEC within Task D of the Cooperative Research Program on Seismic Response of R/C Structures (2nd Phase). The main aim of these experimental tests was to provide information on the nonlinear behaviour of the beams of a four storey full scale building tested in ELSA using the PSeudo Dynamic (PSD) test method. A FiBRE/TIMOSHENKO BEAM ELEMENT IN CASTEM 2000

43 of 55


4φ12mm

4φ6mm

0.10m 3φ12mm

φ6mm // 7mm

0.20m

0.20m

0.40m

Figure 13- Section of the T beam.

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

Figure 14- Response curves of the cantilever T beam.


44 of 55


Z

31.0mm 0.25m

8∅16mm

Y 19.0mm

15.0mm

∅8mm // 70mm 0.25m

Figure 15- Column section tested by Bousias.

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.

6

70.

1 .0

50. 7

30.

5 Disp. UX (mm)

Disp. UY (mm)

3

10.

-140. -2.0

.0

2.0

4.0

-140.

.0

140.

Figure 16- Loading path: axial force and horizontal displacements.


45 of 55


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

Figure 17- Free top force-displacement response curves on Y direction.


46 of 55


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

Figure 18- Free top force-displacement response curves on Y direction.


47 of 55



48 of 55

SUMMARY AND CONCLUSION

8. SUMMARY AND CONCLUSION A fibre Timoshenko beam element has been implemented in CASTEM 2000. The Timoshenko beam theory 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 potentially taken into account. Although a suitable constitutive model is still to be derived, the most 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 followed by a straight line in the softening zone. Confinement is taken into account by a 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 been also considered. • Steel behaviour is represented by a modified Menegoto-Pinto model with a three stage monotonic curve (linear, plateau and hardening). The Bauschinger effect is taken into account and buckling can be simulated. The set of numerical tests included in the report shows the potential of this uniaxial modelling for the study of complex structural elements and loading paths. However some aspects need further improvements; in particular, the modelling of the concrete crack closing. Moreover, a suitable model able to represent the non-linear behaviour of structural elements dominated by shear (e.g. short bridge piers) should be developed. A series of tests programmed for the ELSA laboratory on bridges (bridge piers) will constitute the basis for this development. The implementation of the beam fibre-modelling in CASTEM 2000 has been performed following the standard procedures for model implementation. However, in order to be general enough and to benefit at most from the existing capabilities of the code, a two-level approach has been adopted. • The structural beam level allows to discretize complex boundary value problems, the Gauss integration points of each beam element being described at a lowest level: the section level. • The section level allows not only to specify the geometry and the material characteristics of any section, but also to easily analyse the distribution of stress, strain and internal variables resulting from the straining of this section. As a result of this two-level approach, it is very easy to implement other fibre constitutive laws.

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 d’utilisation, 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., Oñate, 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 d’utilisation, 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

ANNEX 1- A SAMPLE GIBIANE INPUT FOR FIBRE APPLICATION

ANNEX 1- A SAMPLE GIBIANE INPUT FOR FIBRE APPLICATION *-------------------------------------------* * EXAMPLE OF A "T" R.C. SECTION WITH STEEL ON * THE TOP AND BOTTOM OF THE BEAM * (L.N.E.C. TESTS) * *-----------------------------------------* Quadrangular Elements *-----------------------------------------opti dime 2 elem qua4; *-----------------------------------------* Triangular Elements (commented out) *-----------------------------------------* opti dime 2 elem tri3; *-----------------------------------------* DEFINITION OF THE STEEL GEOMETRY (SECTION) *-----------------------------------------ps1 = -1.41372d-2 -17.6d-2; ps2 = -1.41372d-2 -16.4d-2; ps3 = 2.82743d-2 0.0 ; cars1 = ps1 d 1 ps2 tran 1 ps3; * ps4 = -1.88496d-2 7.6d-2; ps5 = -1.88496d-2 6.4d-2; ps6 = 3.76991d-2 0.0 ; cars2 = ps4 d 1 ps5 tran 1 ps6; * ps7 = -30.9425d-2 7.0d-2; ps8 = -30.9425d-2 6.4d-2; ps9 = 1.88496d-2 0.0 ; cars3 = ps7 d 1 ps8 tran 1 ps9; * cars4 = cars3 plus (60.0d-2 0.0); * mesf = coul (cars1 et cars2 et cars3 et cars4) roug; *-----------------------------------------* DEFINITION OF THE UNCONFINED CONCRETE GEOMETRY (SECTION) *-----------------------------------------pc1 = -50.0d-2 10.0d-2; pc2 = -50.0d-2 0.0d-2; pc3 = 40.0d-2 0.0d-2; carc1 = pc1 d 4 pc2 tran 1 pc3; carc2 = carc1 plus (60.0d-2 0.0); * pc7 = -10.0d-2 10.0d-2; pc8 = 10.0d-2 10.0d-2; pc9 = 0.0d-2 -2.1d-2; carc3 = pc7 d 1 pc8 tran 1 pc9; * carc4 = carc3 plus (0.0 -27.9d-2); * mesu = coul (carc1 et carc2 et carc3 et carc4) turq; *-----------------------------------------* DEFINITION OF THE CONFINED CONCRETE GEOMETRY (SECTION) *-----------------------------------------pc10 = -10.0d-2 7.9d-2; pc11 = -10.0d-2 -17.9d-2;


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


2.00000000000000E-03 -6.07153216591883E-18 -2.00000000000001E-03); dtopref=dtopref et (prog -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 -2.20000000000000E-02 -2.40000000000000E-02 -2.60000000000000E-02 -2.80000000000000E-02 -3.00000000000000E-02 -3.20000000000000E-02 -3.40000000000000E-02 -3.60000000000000E-02 -3.80000000000000E-02); dtopref=dtopref et (prog -4.00000000000000E-02); ftopref=prog .00000000000000E+00 -1.14481311676523E-03 -1.22901793761370E-02 -2.32933890701765E-02 -3.41409060368847E-02 -4.48163634737750E-02 -5.52993760011572E-02 -6.25736672822648E-02 -6.34682582020739E-02 -6.37997368461639E-02 -6.39516852406453E-02 -6.40648825527967E-02 -6.41051047602156E-02 -6.40881950555459E-02 -6.40683885171247E-02 -6.40457149947723E-02 -6.40198575442582E-02 -5.31916838295068E-02 -4.24307708413968E-02 -3.18940241565422E-02 -2.15933027921714E-02; ftopref=ftopref et (prog -1.14642689013442E-02 -1.95490782397935E-03 4.34527534900837E-03 1.05039345010342E-02 1.64943367735540E-02 2.22738649284235E-02 2.77888455743349E-02 3.11954965063240E-02 3.11996112204543E-02 3.12036823442664E-02 3.12076823925978E-02 3.12115835624900E-02 3.12155450724221E-02 3.12193553958296E-02 3.12232037374745E-02 3.12269455176290E-02 3.13243163279959E-02 3.17866134984408E-02 3.22108792714725E-02 3.26000050114313E-02 3.29579270685599E-02); ftopref=ftopref et (prog 3.32842841187497E-02 3.35810040717798E-02 3.38496367881653E-02 3.40918486154093E-02 3.43092715455886E-02 2.76024463612365E-02 2.10585068426046E-02 1.46524661635717E-02 8.49910700480509E-03 2.69258722007945E-03 -2.62926487360908E-03 -8.01317688346363E-03 -1.25038980801999E-02 -1.60322003928289E-02 -1.87262001823851E-02 -2.07658260777817E-02 -2.23232577557867E-02 -2.35408842921363E-02 -2.45206665833014E-02 -2.53333845014858E-02 -2.60273227929928E-02); ftopref=ftopref et (prog -2.66354492379779E-02 -2.71804925526341E-02 -2.76783464385084E-02 -2.81403048598080E-02 -2.85745259272785E-02 -2.89869974898484E-02 -2.93821816783630E-02 -2.97634513185043E-02 -3.01333896091121E-02 -3.42706858331460E-02 -4.38448863358707E-02 -5.31853566475303E-02 -6.18463597099409E-02 -6.72796758934449E-02 -6.94003595649139E-02 -7.03019196590381E-02 -7.08358037573039E-02 -7.12371916595880E-02); ftopref=ftopref et (prog -7.16045522313159E-02); fdtopr=evol rouge manu ’DISPL. (m)’ dtopref ’FORCE (MN)’ ftopref; titre ’Computed curve (white) and reference curve (red)’; dess (fdtop et fdtopr) tt;


55 of 55

Fibre Beam Element

Recommend Documents