- Email: [email protected]
Efficient approaches to finite element analysis in earthquake engineering
Computers and Structures 81 (2003) 1223–1239 www.elsevier.com/locate/compstruc
Efficient approaches to finite element analysis in earthquake engineering L. Davenne a
a,*
, F. Ragueneau a, J. Mazars b, A. Ibrahimbegovic
a
Laboratoire de M ecanique et Technologie, Ecole Normale Sup erieure de Cachan, 61, Avenue du Pr esident Wilson, 94235 Cachan Cedex, France b Laboratoire Sols, Solides, Structures, INPG––BP 53, 38041 Grenoble Cedex 9, France
Abstract This paper deals with the modeling of reinforced concrete structures subjected to earthquake ground motion. Due to the complex behavior of both materials and structures, efficient numerical tools are developed herein in order to keep accuracy and robustness for large scale computations. We focus our attention on the use of simplified multifiber beam element describing the response of structural components and on macro-element accounting for soil–structure interaction. 2003 Elsevier Science Ltd. All rights reserved.
1. Introduction Earthquake engineering design of a civil engineering structure concerns in general several levels of strong motion earthquakes, each one with its specific design requirements. In particular, one should distinguish between: – A damage limitation criterion which allows to keep the operational features of the building intact for a seismic loading which has a large occurrence probability. – A collapse criterion which insures that for the nominal design seismic action the structure remains stable without significant local collapse. – A life safe objective requiring that for the strongest possible earthquake feasible for the site, the global collapse is prevented and the life is preserved. In the traditional approach to earthquake engineering design, the computations are carried out on the basis
* Corresponding author. Tel.: +33-1-4740-2243; fax: +33-14740-2240. E-mail address: [email protected] (L. Davenne).
of linear elastic static analysis. The nonlinear behavior and energy dissipation, can be accounted for in a trivial manner by a force-based approach, where the level of seismic loading is computed by analyzing the elastic response spectra on a single degree of freedom system and the so-called R reduction factor method is then used to introduce the ductility of materials (steel and concrete). A more rational concept developed during this last decade, named displacement-based procedure [19], turns towards the design based on the critical limit states of the structural elements. In other words, defining the ultimate behavior of materials, one has to deal with maximum strains or curvatures, thus taking into account the structure ductility demand in a more physical way. Whereas the force-based design uses the elastic properties of the structure, the displacement-based design characterizes it through the secant stiffness at maximum displacement. The evaluation of vibration frequencies during the loading is more appropriate than for the purely elastic analysis. These classical procedures for design (force or displacement formulation) are based on the study of an equivalent single degree of freedom system or in fact a simplified substitute structure which is not capable to account for the load redistribution inside the structures due to local nonlinearities. This is one of the major drawbacks preventing a realistic description of the
0045-7949/03/$ - see front matter 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0045-7949(03)00038-5
1224
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
global and local behavior of a structure up to failure. One remedy consists in using nonlinear finite element method to perform push-over analysis, i.e. a nonlinear static analysis under monotonic increasing lateral load. It allows to determine the maximum carrying capacity of the structure in terms of forces, displacements, ductility (deformation), crack pattern and failure mode. The main assumption of such a computation is that the response is related to an equivalent single degree of freedom and thus, controlled by a single mode which justifies that the position and amplitude of horizontal lateral loads are typically defined by an elastic analysis. An alternative choice to perform the earthquake resistant design is by making use of the nonlinear timehistory analysis, assuming a physical description of materials and applying transient loadings on the structure in terms of natural or simulated ground motion. The evolution of eigenmodes concomitant to the stiffness degradation governed by local yield criterion provides currently the most refined method of analysis for ultimate behavior of concrete structure [16,23]. For a simple reason of excessive computational costs, such approaches to structural dynamics in civil engineering are not commonly used. Several strategies have been developed in order to lower the CPU time in transient analysis by solving the equations of motion by a mode superposition techniques even in the nonlinear range, either by dynamic substructuring for limited spread of nonlinear zones (e.g. see [10,11]) or using the tangent stiffness to modify the basis [14]. Nonlinear dynamic analysis of complex civil engineering structures based on a detailed finite element model requires large scale computations and handles delicate solution techniques. The necessity to perform parametric studies due to the stochastic characteristic of the input accelerations imposes simplified numerical modeling which will reduce the computation cost. In this work the latter is achieved by selecting the classical Euler–Bernoulli beam model for representing the global behavior of the structural components of a complex civil engineering structure. We note that no simplification is made other than this choice of the beam model. In particular, the constitutive bahavior of the chosen model remains sufficiently general to take into account all the different inelastic phenomena (cracking by damage, irreversible deformation by plasticity and crack-closing by unilateral frictional contact condition). Another effect which is also taken into account concerns the interaction of structure with the foundation, which provides both a more realistic boundary condition as well as the damping model due to radiation effect. The outline of the paper is as follows. In Section 2 we briefly describe the chosen model of multilayer beam used for representing the structural components. The development pertaining to the modeling of structure– foundation interaction are described in Section 3. In
Section 4 we present several numerical simulations which illustrate the performance of the proposed models with respect to the experimental results. Concluding remarks are given in Section 5.
2. Multifibers beam elements for nonlinear transient analysis of reinforced concrete structures In this work we seek to develop the models capable of representing in a reliable manner nonlinear behavior of a structure damaged by a strong earthquake, with the main goal of providing a physically based description of damping (as opposed to a more convenient, but very adhoc Rayleigh damping model). This kind of approach is especially suitable for a novel concept in earthquake resistant design of shear-wall type reinforced concrete structures; referred nowadays as ÔFrench wallsÕ, where the choice of reinforcement is made so that the damaged and cracking spreads as much as possible through the structure which allows that the main part of the energy generated by an earthquake be absorbed through inelastic material behavior. In short we seek to develop the structural models capable of integrating the advanced damage models for concrete (e.g. [17]), in order to describe the nonlinear behavior of a reinforced concrete structures. With respect to the large spreading of the zone with nonlinear behavior we further seek to limit the model complexity (and resulting computational costs) by limiting the diversity of possible deformation global patterns which is achieved by choosing a multifiber beam model, with all the fibers restricted to beam kinematics and with each fiber employing its own constitutive model (see Fig. 1). The main advantage of using a multifiber type finite element concerns the simple uniaxial behavior which allows a very efficient implementation of inelastic constitutive. This is no longer possible for thick beams where shear strains play a major role [7]. 2.1. Cross-section behavior The multifiber beam element developed herein employs the standard Hermite polynomial shape functions to describe the variation of the displacement field along the beam. The difference with ‘‘classical’’ beam elements [21] concerns the cross-section behavior, that is the relation between the generalized strains e and the generalized stresses s. In the general 3D case the latter includes: s ¼ ðN
Mx
My
Mz ÞT
and e ¼ ð e hx
vy
vz ÞT ð1Þ
where N is the normal force, Mx the torque, My and Mz are the bending moments, e the axial strain, hx the twist,
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
1225
Fig. 1. Multifibers beam formulation for reinforced concrete structures.
vy and vz the curvatures. The cross-section behavior is expressed with the constitutive matrix: 2 3 0 C13 C14 C11 6 C22 0 0 7 7 C¼6 ð2Þ 4 C33 C34 5 sym C44 where the coefficients are obtained by integrating over the cross-section (y and z axes): Z Z E dS; C13 ¼ Ez dS; C11 ¼ S S Z Z C14 ¼ Ey dS; C22 ¼ Gðy 2 þ z2 Þ dS S S Z Z C33 ¼ Ez2 dS C34 ¼ Eyz dS S ZS 2 C44 ¼ Ey dS ð3Þ S
where E and G are YoungÕs and shear moduli which vary in y and z. The chosen moduli can be either initial, secant or tangent, depending upon the iterative algorithm used to solve the global equilibrium equations. The component of the constitutive matrix are computed by numerical integrations in (3), with one Gauss point per fiber. For the Euler–Bernoulli element, the shear forces are computed at the element level through the equilibrium equations (included in the Hermite polynomial shape functions). Reinforcement bars are introduced within special fibers, whose behavior is obtained as a combination of those of concrete and steel according to: rlayer ¼ ð1 aÞrconcrete þ arsteel
ð4Þ
where a is the relative area of the reinforcement in the layer. When dealing with structures which slenderness ratio is far from the classical beam such as shear wall, a more reliable representation of shear deformations and shear stresses has to be provided. One possibility in that respect is to use the classical Timoshenko beam model, which can describe the constant shear strain. The main
difficulty of developing the finite element implementation of the Timoshenko beam model concerns the socalled shear locking phenomena (e.g. [26]), or inability of the standard finite element approximations to represent pure bending vanishing shear modes. A number of different remedies to shear locking problem has been proposed (e.g. see [9]), ranging from selective or reduced integration, assumed shear strain, enhanced shear strain or hierarchical displacement interpolations. A very recent work of Kotronis [12] extends these ideas in order to construct shear locking remedies for a multifiber Timoshenko beam. Remark: In order to retain the simplicity of the constitutive model for the Timoshenko formulation, the normal and shear stress resultants are uncoupled and nonlinearities are only introduced in the stress–strain relationship in the direction of the normal to the crosssection, whereas the tangential behavior is assumed to remain elastic. A more complex constitutive relationships can be used for each fiber, expressed within the framework of irreversible processes thermodynamics which allows physical coupling of all stress components. Such a development for 2D Timoshenko multilayer beam was carried out by Dube [7], where the shear stress as well as the normal stress in each layer was used to compute the damage of the layer with the 2D version of the damage model presented hereafter. 2.2. Enhanced strain formulation For a multifiber beam model, as indicated in (2), the nonlinear behavior implies coupling between axial force N and bending moment M. It is thus necessary that the axial force and bending moment have the same variation along the element. The latter is not satisfied for an Euler–Bernoulli multifiber beam element, where axial strain remains constant along the element whereas the bending strain is linear. This problem may become important in the nonlinear range where a linear curvature along the beam implies a linear variation of the axial strain due to the shift of the neutral axis.
1226
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
x x u2 u1 u1 þ u2 ) eðxÞ ¼ ¼ cst uðxÞ ¼ 1 L L L
x2 x3 x3 x2 2 þ x h1 vðxÞ ¼ 1 3 2 þ 2 3 v1 þ L L L2 L
3
3 2 x x x x2 h2 ) jðxÞ þ 2 3 þ 3 2 v2 þ L L L2 L
6 12x 6x 4 ¼ 2 þ 3 v1 þ h1 L L L2 L
12x 6 6x 4 h2 þ 3 þ 2 v2 þ L L L2 L
u_ tþDt ¼
ð5Þ One way to remove this incompatibility of variations for axial strain and curvature is by using enriched description of the axial strain field with: eðxÞ ) eðxÞ þ ~eðxÞ;
~eðxÞ ¼ GðxÞa;
4 8x GðxÞ ¼ 2 L L ð6Þ
where GðxÞ is an enhanced axial strain which can be derived from the displacement bubble function. The variational basis of such an enhancement is provided by Hu–Washizu principle (e.g. [26]), which can be presented in the same manner as the modified method of incompatible modes of Ibrahimbegovic and Wilson [10,11] Z Z 0 ¼ ðe ðxÞN ðxÞ þ j ðxÞMðxÞÞ dx ¼ u ðxÞf ðxÞ dx L L Z ~e ðxÞN ðxÞ dx ð7Þ 0¼ L
2.3. Global solution procedure Dynamic analysis for earthquake ground motion for this kind of structural model reduces to solving the set of nonlinear equations which can be written as: M€ uðtÞ þ C u_ ðtÞ þ f int ðu; tÞ ¼ Mr€ug ðtÞ
ð8Þ
where M and C are, respectively, mass and damping matrix, €uðtÞ and u_ ðtÞ are nodal accelerations and velocities, f int ðu; tÞ is the internal force vector obtained by discretizing (6) and r€ug ðtÞ is the ground acceleration effect applied on the structure. By using a time-integration scheme, the differential equation of motion in (8) is reduced to an algebraic equation. In particular for low frequency response in earthquake engineering analysis one uses an implicit scheme such as Newmark one-step scheme (see [18]) which allows to express the velocity and acceleration vectors at time t þ Dt as the functions of their corresponding values at time t and incremental displacement vector: € utþDt ¼
1 1 1 2b €uðtÞ u_ ðtÞ Du b Dt2 b Dt 2c
ð9Þ
c c 1 2b Du þ 1 u_ ðtÞ u_ ðtÞ b Dt b 2c
c € þ Dt 1 uðtÞ 2b
ð10Þ
where c ¼ 1=2 and b ¼ 1=4 are typically chosen for optimal result accuracy. The discrete set of equations, obtained by introducing (9) and (10) into (8), is further solved by an iterative solution procedure, where at each iteration the problem reduces to
1 c C þ K Dui ¼ Mr€ M þ ug ðt þ tÞ b Dt2 b Dt M€ ui ðt þ DtÞ C u_ i ðt þ DtÞ f int ðui ðt þ DtÞÞ uiþ1 ðt þ DtÞ ¼ ui ðt þ DtÞ þ Dui
ð11Þ
Different iterative strategies can be used to compute the tangent stiffness matrix K in (11) above; if we seek to reduce the computational cost one would not update K at each iteration as for NewtonÕs method, but either keep it constant as for modified NewtonÕs method or use a secant stiffness for quasi-Newton method. The details of computing K and f int pertain to the constitutive model which is specified for each component as explained subsequently. 2.4. Materials constitutive relations for dynamics Both steel and concrete materials are described within the thermodynamic framework for irreversible processes (e.g. see [15]). In describing the nonlinear behavior of reinforcement bars, we choose the classical plasticity model and take into account the nonlinear kinematic hardening of Armstrong and Frederick [2] in order to be able to better describe the observed hysteresis loops. The free energy for this model can be written as: qw ¼ 12ðe ep Þ : H : ðe ep Þ þ 12ba : a
ð12Þ
in which H is the HookeÕs elasticity tensor, ep is the plastic strain and a the hardening internal variable which describes hardening. The constitutive equations for this kind of model can be written as: r¼
oðqwÞ ¼ C : ðe ep Þ; oe
X ¼
oðqwÞ ¼ ba oa
ð13Þ
where X is the stress-like hardening variable. The latter is used to describe a modified form of the plasticity criterion: f ¼ J2 ðr X Þ þ 34aX : X ry 6 0
ð14Þ
where a, b and ry are material parameters. Due to the particular geometric characteristics of steel bars, only a 1D implementation of the model is carried out. A typical
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
1227
Fig. 2. Hysteresis loop for the steel modeling.
stress–strain hysteresis loops predicted by this model is given in Fig. 2. The constitutive model for concrete for earthquake engineering ought to take into account some observed phenomena, such as decrease in material stiffness due to cracking, stiffness recovery which occurs at crack closure and inelastic strains concomitant to damage. To simulate this behavior we employ the damage model of La Borderie [13], which uses two scalar damage variables, D1 for damage in tension and D2 for damage in compression. Unilateral effect and stiffness recovery (damage deactivation) are also included. Inelastic strains are taken into account thanks to an isotropic tensor. The Gibbs free energy of this model can be expressed as: v¼
hriþ : hriþ hri : hri t þ þ ðr : r Trðr2 ÞÞ 2Eð1 D1 Þ 2Eð1 D2 Þ E b 1 D1 b2 D2 f ðrÞ þ TrðrÞ þ Eð1 D1 Þ Eð1 D2 Þ
ð15Þ
where b1 and b2 are material coefficients, whereas h iþ and h i denotes the positive or negative values of the given variable. f ðrÞ and rf are the crack closure function and the crack closure stress respectively. E is the initial YoungÕs modulus and m the Poisson ratio. For the 1D expression, the total strain is: e ¼ ee þ ein rþ r þ Eð1 D1 Þ Eð1 D2 Þ b 1 D1 b2 D2 f ðrÞ þ ein ¼ Eð1 D1 Þ Eð1 D2 Þ ee ¼
ð16Þ
where ee the elastic strain and ein the inelastic strain. The partition of the stress is obtained as:
r > 0 ! rþ ¼ r; r < 0 ! rþ ¼ 0;
r ¼ 0 r ¼ r
and for the crack-closure function:
ð17Þ
Fig. 3. Unidimensional constitutive relations for concrete.
8 < r > 0 ! f ðrÞ ¼ 1 rf < r < 0 ! f ðrÞ ¼ 1 rrf : r < rf ! f ðrÞ ¼ 0
ð18Þ
where rf is the stress where the cracks are completely reclosed. The role of the crack closure function is: – To restore the material stiffness, passing from tensile to compressive behavior. – To null the tensile inelastic strains creating after tension when excessive (over rf ) compression is reached. Damage criteria are expressed as fi ¼ Yi Zi (i ¼ 1 for tension or 2 for compression, Yi is the associated force to the damage variable Di and Zi a threshold dependent on the hardening variables). The evolution laws for the damage variables Di are written as Di ¼ 1
1 1 þ ½Ai ðYi Y0i Þ Bi
ð19Þ
where Y0i is the initial elastic threshold and Ai , Bi are material constants. A typical stress–strain response produced by this model for a uniaxial cyclic loading which passes through traction, compression and again traction is given in Fig. 3. The material coefficients Y0i , Ai ; Bi , bi (i ¼ 1, 2) and rf can be identified from uniaxial tests. They are only function of intrinsic material characteristics (aggregate size, concrete mix, . . .).
3. A soil–structure interaction macro-element It has already been noted that, after severe earthquakes, some structures with different properties from
1228
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
surrounding ground have been preserved from any damage. The interaction of the structure with the soil can lead to a less severe response of the structure either by modifying the support condition with respect to a clamped boundary or accounting for modification of the properties of the interacting system. In particular the low stiffness of the soil reduces the system frequency, and may, depending on the frequency content of the input motion, drive outside of the resonance zone of the input spectrum. Moreover, the nonlinearities in the surrounding soil and interface (yielding of the soil and uplift of the structure) lead to a strong energy dissipation and further reduce the forces in the structure. To model these phenomena we developed a simplified method for soil–structure interaction. The overall linear and nonlinear behavior in the soil and at the interface is reduced through a macro-element expressed in three degrees of freedom on the strip foundation. The constitutive laws governing the behavior of this macroelement are developed within the framework of the plasticity theory. The quasi-static cyclic loading model is developed first (see [4,5] followed by an extension to dynamic loading case [6]. The only variables of the macro-element are the forces at the base of the foundation (vertical force V , horizontal force H and moment M) and the corresponding displacements (vertical displacement z, horizontal displacement x and rotation h) measured at the centre of the foundation (assumed as a rigid body). These variables normalized by characteristic values to become dimensionless.
0 1 1 V V0 1 B C B C F @ H0 A @ H A and Bqmax 0 M=B M 0 01 0 1 z z B C 1B C u @ x0 A @ x A B h0 Bh 0
ð20Þ
where B is the width and Vmax ¼ Bqmax the bearing capacity of the foundation. Since it is a strip foundation, all terms are expressed for a unit transverse length (unit depth of the footing). Three different phenomena are represented by the forces acting on the foundation: the plasticity of the soil surrounding the foundation or near field, the uplift of the foundation and the elasticity of the soil in far field. They are modeled through three interacting models (Fig. 4). Plasticity of the soil and uplift are strongly nonlinear and coupled since they both concern the near field while elasticity concerns the far field. The total displacement increment is partitioned into elastic, plastic and uplift components, which can be written choosing the corresponding rate decomposition: u_ tot ¼ u_ el þ u_ pl þ u_ up
ð21Þ
The global model of constitutive behavior can be written: F_ ¼ ð½Kelpl 1 þ ½Kup 1 Þ1 : u_ tot
ð22Þ
where elastoplastic tangent stiffness matrix Kelpl is obtained from the plasticity model and the tangent stiffness matrix Kup from the chosen ‘‘uplift model’’.
Fig. 4. Idealization of the complete model.
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
1229
3.1. Elasticity and damping
3.2. Plasticity model
The energy dissipation mechanism in the near field arises exclusively from the nonlinearities included in the plastic and uplift cyclic models. The energy dissipation in the far field arises from the radiation boundary conditions in semi-infinite medium. The latter is modeled through the imaginary part of the elastic impedances [8], represented by a dashpot, and the real part is the elastic stiffness matrix. 0 0 1 Kzz 0 0 0 0 A Kel ¼ @ 0 Kxx 0 0 0 Khh 1 0 K zz 0 0 C B qmax C B Kxx B 0 C ð23Þ ¼B 0 C C B qmax A @ Khh 0 0 B2 qmax
Finite element analysis with an efficient model for the soil [4] or analytical results, when available, guided us in constructing the macro-model presented in the following. The elastoplastic model provide the elastic and plastic parts of the displacement increments:
0:73 2 2 G0 ð1 þ 2aÞ; Kxx ¼ G0 1 þ a ; 1m 2m 3
2
p B 1 Khh ¼ G0 ð24Þ 1þ a 2ð1 mÞ 2 3 Kzz ¼
where a is the gradient of the variation of the shear modulus G ¼ G0 ð1 þ afÞ (G0 the shear modulus at z ¼ 0) with respect to the depth f ¼ 2z=B and m the PoissonÕs ratio. In the same way, the coefficient of the damping matrix are [8]: 3:4 Ab ; pð1 mÞ 3:4 Ibh Chh ¼ qVs0 pð1 mÞ
Czz ¼ qVs0
F_ ¼ Kelpl : ðu_ el þ u_ pl Þ with
Q¼
ofi ; oF
ð27Þ
where fi is the yield surface and g the plastic potential for nonassociated flow rule employed herein. One also defines a failure surface from analytical [20] and numerical results [25], based on limit analysis for structure overturning with uplift mechanisms (the uplift influence is also included in the plasticity model). !2 !2 H0 M0 fi ¼ þ 1¼0 aV 0c ð1 V 0 Þd bV 0e ð1 V 0 Þf ð28Þ For a soil with constant gradient of cohesion with depth, the coefficients in the last expression are given as: a ¼ 0:32=x; e ¼ 0:8;
where q is the mass density of the soil, Ab the areapof theffi ffiffiffiffiffiffiffiffiffiffi foundation, Ibh its quadratic moment and Vs0 ¼ G0 =q the shear velocity at z ¼ 0. Remark: Step by step analysis prevents from including the frequency variation of the dynamic impedances. Constant values are thus chosen (static value for stiffness, which are already given for the model, and the value near resonance of the soil structure system for damping, which gives the better average damping). This choice is acceptable for a semi-infinite soil medium with a slow and continuous properties variation (nearly homogeneous soil). In this case the variation of the coefficients is low over the frequency range of interest (0–10 Hz). For a foundation lying on a thin soil layer, this approximation is no longer valid. In this case one could take the values corresponding to the fundamental frequency of the soil building system.
1 ðKel : PÞ ðQ : Kel Þ; h0 þ h0 og P¼ oF
Kelpl ¼ Kel
Cxx ¼ qVs0 Ab ; ð25Þ
ð26Þ
b ¼ 0:37=x0:2 ;
c ¼ 0:25; qmax f ¼ 0:8 with x ¼ q0 max
d ¼ 0:55;
where qmax is the actual ultimate bearing capacity of the foundation (soil with cohesion gradient c ¼ c0 þ gz) and q0 max the ultimate bearing capacity for an homogeneous soil with a constant cohesion c0 . The model is particularly well suited for the evolution of the yield surface towards the failure surface with a vertical force being almost constant, for reproducing the behavior of the soil and the foundation, initially submitted to the weight of the structure, and then loaded mainly in the (H ; M) plane (seismic loading).
fi ¼
H 0 aCcH qCcH
2
þ
M 0 bCcM qCcM
2 1¼0
ð29Þ
with CcH ¼ aV 0c ðc V 0 Þd ; CcM ¼ aV 0e ðc V 0 Þf ; c ¼ v þ ð1 vÞðq þ sÞ; v ¼ N =Vmax and N is the weight of the structure (initialisation under dead load). In (29) above s ð 0 a b ÞTqisffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the kinematic hardening vector with s ¼ ðs : sÞ1=2 ¼
parameter.
a2 þ b2 and is the isotropic hardening
1230
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
The consistency rule, the criterion of noninterpenetration of the failure surface with the loading surface and a relationship obtained from the observation of the foundation behavior during uplift allow us to define the evolution of the isotropic and kinematic hardening parameters. The kinematic hardening increment is defined by its amplitude s_ and its direction l: s_ ¼ s_
l jlj
with s_ ¼ hk_i ofi os
h0 l
i : jlj þ of oq
and
1 k_ ¼ 0 Q : F_ h
0
oz0up B oV 0 B F_ ¼ Kup : u_ up ¼ ðCup Þ1 : u_ up ; Cup ¼ B 0 @ 0up oh oV 0
The displacements induced by uplift are given as:
z0up 1 0 1 0 C cos w H B H C B C c C and w ¼ a tan C lB H B 1 C 0 @ CM sin w M A
ð35Þ C1H M 0 bC1M
C1M H 0 aC1H
!
CcM
ð31Þ The isotropic hardening increment is set equal to the kinematic one: q_ ¼ s_ . The plastic modulus h0 is fitted on the M 0 –h0 curves since it is the most important behavior of the foundation
K el ðdM 0 =dh0 Þ2M 0 of dV 0 of dH 0 of h0 ¼ elhh þ þ Khh ðdM 0 =dh0 Þn2 oV 0 dM 0 oH 0 dM 0 oM 0 ð32Þ K 0 h0 hh0 M1
dM 0 0 ¼ Khh exp dh0
!
The external normal to the yield surface is not convenient for defining the direction of the plastic displacements and a nonassociated flow rule is required. The plastic potential g has been chosen as:
g¼
H j
0 2
þ
M n
0 2
þ V 02 1 ¼ 0
ðd=dmax Þ2 and 1 d=dmax
jh0 j d=dmax þ lnð1 d=dmax Þ ¼ ð1 V 0 Þ 2 1 d=dmax
h0up ¼ ð1 V 0 Þh0
ð30Þ where 0
1 oz0up oM 0 C C 0 0 C 0up A oh 0 oM 0 ð34Þ 0
ð33Þ
3.3. Foundation uplift model The foundation uplift behavior on an elasto-plastic soil is a very complex phenomenon which depends strongly on the level of soil yielding. A cyclic model which takes into account the irreversibility of the uplift behavior due to plasticity is developed as a direct modification of the uplift model of Cremer [4] developed for elastic soils. The uplift model operates only on the uplift part of the total displacement introducing coupling between the vertical displacement zup and the rotation hup of the foundation:
0ð0Þ M0 j
0 is the uplift reached where dmax ¼ ð4=V 0 ÞjMmax for the maximum moment at failure for the given V 0 and 0 for H 0 ¼ 0, that is to say for Mmax ¼ bV 0e ð1 V 0 Þf . 0 For a given V , the uplift ratio d is obtained as follows: 0ðiÞ
• if jM 0 j 6 jM0 j, d ¼ 0 0ð0Þ 0ð0Þ • if jM 0 j P jMp0ðiÞ j, d ¼ ð4=V 0 ÞjM 0 M0 j with M0 ¼ 0 ðV 0 =4Þ expAV 0ðiÞ 0ðiÞ • if jM0 j < jM 0 j < jMp0ðiÞ j, d ¼ ð4=gV 0 ÞjM 0 M0 j 0ðiÞ
0ðiÞ
0 0ðiÞ where M0 is such that dðiÞ p ¼ ð4=gV ÞjMp M0 j ¼ 0ð0Þ ð4=V 0 ÞjMp0ðiÞ M0 j with Mp0ðiÞ the maximum moment reached during the loading history and dðiÞ p is the corresponding uplift. The right and left uplift are treated independently supposing they remain uncoupled. Note: The macro-element treats uplift in a continuous global model. No coefficient of restitution to describe the energy dissipation upon impact while the uplifted foundation will contact the soil has been introduced.
3.4. Use of the model This SSI macro-element which was implemented in computer code FEAP [24] is well suited for the dynamic simulation of slender structures with potentially large uplift and moderate yielding in the soil (and no sliding of the footing on the soil). When using this model it is recommended to have a safety factor SF ¼ Vmax =V > 3 (where Vmax is the bearing capacity of the foundation, V the vertical load on the foundation), and a slenderness ratio h=B > 1 (where h is the height of the dynamic centre of gravity and B the width of the foundation).
4. Numerical simulations and experimental results (CAMUS program) In this section we briefly present the main results of the experimental program carried out within the re-
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
search network CAMUS along with the comparative numerical simulations. The main purpose of this experimental program consists in demonstrating the ability of lightly reinforced concrete bearing walls to sustain very severe earthquake. The design is based on a concept favoring damage spreading over several storeys of a lightly reinforced concrete wall. This kind of design leads to a lower percentage of reinforcements with an optimized distribution generating a wide cracks pattern allowing for dissipation of great amounts of energy. As a consequence, one obtains a vertical rising of the masses resulting in energy transformation (from kinematic to potential instead of storing strain energy in the structure) and an increasing ductility thanks to this particular way of dissipating the earthquake energy. In order to test this concept experimentally, a one-third scaled model (see Fig. 5) has been installed on the shaking table of CEA, composed of two parallel braced walls linked by six square slabs. A massive reinforced concrete footing allows the anchorage to the shaking table. The similarity laws to a full-scale structure and the mock-up require additional masses of 6.55 t at each storey. The problem is assumed to be two-dimensional since the mock-up is loaded through horizontal acceleration parallel to the walls and the presence of steel bracing systems at each storey placed perpendicularly to the loading direction prevent occurrence of any torsional modes. pffiffiffiThe accelerograms are scaled in time with a ratio of 1= 3 in order to account for the similarity rules. Two types of accelerograms are applied: a synthetic one generated from the french code [1] Nice S1 as the far field type earthquake and the San Francisco earthquake of 1957 (Fig. 6) as the near field one. The response spectra show clearly the main difference of these two earthquakes: on one hand Nice has a broad frequency content and on the other hand, San Francisco has a
1231
narrower effective band width. The complete experimental sequence was: Nice 0.25 g, San Francisco 1.13 g, Nice 0.4 g and Nice 0.71 g. 4.1. Model calibration (modal analysis and damping calibration) Eigenfrequencies of the reduced scale numerical model anchored to the shaking table show a higher global stiffness. Namely the frequency of the first eigenmode measured before the test is 7.3 Hz, as opposed to 10.3 Hz predicted by the numerical model. This indicates the necessity to take into account the shaking table compliance as well as the anchorage system which can never be considered as perfectly clamped. The second eigenmode corresponds to pumping with an eigenfrequency of 40 Hz as predicted by the finite element model. During the test, the vertical displacement of the mass induced by the openings of cracks excited this second mode. The resulting vibration measured through the vertical acceleration of the shaking table, with a frequency of 20 Hz leads to variation of the vertical dynamic forces which can be very important for reinforced concrete structure failure. This effect must clearly be taken into account through the vertical stiffness of the connecting rods. A simplified model for the basement footing is assumed consisting in an horizontal beam whose properties are selected in accordance with the shaking table stiffness with Kv ¼ 48EI=l3 and Kh ¼ 12EI=l. This kind of footing element adds two eigenmodes, as required for the CAMUS structure, since as already discussed, the second vertical mode plays a very important role. Finally, YoungÕs modulus of the footing beam was set to 10 000 MPa in order to introduce anchorage and contact defects. Table 1 summarizes the correction of fundamental frequencies of the system due to these new boundary conditions.
Fig. 5. CAMUS structure.
1232
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Fig. 6. Nice and San-Francisco accelerograms and response spectrum.
4.2. Comparisons between experimental and numerical results
Table 1 Eigenmode adjustment
First mode Second mode
Initial modeling (Hz)
Adjusted modeling (Hz)
Measured frequencies (Hz)
10.3 40.0
7.4 19.0
7.3 20.0
Damping is another important concern. In particular, the classical viscous damping is assumed to account for various sources of energy dissipation not related to concrete cracking [22]. Despite the lack of clear physical interpretation, the damping is often introduced in the analysis through the Rayleigh damping matrix which is expressed as a linear combination of the mass matrix and the initial stiffness matrix: C ¼ aM þ bK The particular values of damping coefficients are set to ensure a relative damping value of 1% in the first mode and 2% in the second mode. Great care is taken to keep these damping values as stable as possible during all the analysis, despite the fact that cracking induces reduction of the stiffness and shift of the fundamental frequency.
The complete experimental sequence of loading (four earthquakes) was computed using the proposed numerical model, with the material parameters following: YoungÕs modulus of concrete Ec ¼ 30 000 MPa for a maximum compressive strength of fc ¼ 35 MPa and tensile strength ft ¼ 3 MPa, YoungÕs modulus of steel Es ¼ 200 000 MPa with yield stress ry ¼ 414 MPa accounting for hardening with ultimate stress ru ¼ 480 MPa. The numerical and experimental results are compared in terms of time-history displacements, ductility and local damage criterion (Figs. 7–9). For the lowest level of loading (Nice 0.24 g), when the concrete cracks in tension and the reinforcement bars reach the yield stress, the global behavior of the structure is well represented by the numerical model. The loss in stiffness during the later stages of loading and the decrease of the fundamental frequency is also predicted in a fairly acceptable manner, with a maximum displacement always slightly underestimated before and overestimated after the maximum of the load. Table 2 summarizes the comparisons of experimental and computed results for reaction forces for different
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Fig. 7. Nice 0.24 g: experiment/computation comparisons.
Fig. 8. San-Francisco 1.13 g: experiment/computation comparisons.
1233
accelerograms which indicates globally a very good agreement. Another feature of material modeling during cyclic loading, which is barely taken into account, is the cracks closing there described by the crack closure function of the constitutive relation. When re-closing the cracks, the point where the compressive stiffness is totally recovered depends on the crack closure stress of the model rf (see Fig. 10). In the CAMUS experiment, the shocks induced by cracks closure which activates the vertical eigenmode leads to variation of vertical dynamic forces. Experimentally, the phenomenon can be quantified by the measurement of an induced vertical acceleration at the shaking table. The discontinuity of the stiffness when rf is reached reproduces the shock, two computations performed with different forms of the crack closure function (Figs. 11 and 12) allow to reproduce this experimental observation. Modeling this kind of feature become very important for reinforced concrete structures where the design takes into account the coupling between flexural and axial loading. Through this short analysis we can easily remark that the local behavior of concrete plays a very important role on the global structural response and that the material parameters of the crack closure functions can be identified using different levels of loading (see Table 3). Due to cracking of the concrete, the global stiffness of the structure decreases. This variation does not remain
Fig. 10. Effect of the crack closure stress on the slope of the stiffness recovery during cyclic loading.
Fig. 9. Nice 0.71 g: experiment/computation comparisons.
Table 2 Global response comparisons
Nice 0.24 g SF 1.1 g Nice 0.4 g Nice 0.7 g
Displacement (cm)
Shear load (kN)
Moment (kN m)
Vertical l (kN)
Exp.
Comp.
Exp.
Comp.
Exp.
Comp.
Exp.
Comp.
0.72 1.2 1.35 4.4
0.61 1.1 1.1 3.9
65.9 106 86.6 111
65 90 75 120
200 280 280 345
200 240 240 380
202 271 217 312
190 270 225 310
1234
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Fig. 11. Initial model: vertical load variation. Fig. 13. Nice 0.24 g: displacement spectrum.
1.5
Fig. 12. Modified model: vertical load variation.
constant during the loading but depends on the cracks opening situation (opened or closed crack). This stiffness decrease shifts greatly the natural frequencies. The analysis of the fundamental frequency of the structure subject to a white noise after each sequence is not satisfactory because it is applied on the cracked structure but for closed cracks (this is due to the low level of the white noise which does not allow to open cracks). Such experimental measures cannot capture this kind (changes in stiffness reduction during the earthquake) of structural behavior and provide the first eigenfrequency of 6 Hz instead of the 3 Hz., during the strongest part of
Fig. 14. Nice 0.71 g: displacement spectrum.
the earthquake motion. One way to estimate the dominant frequency during the loading is to plot the response spectrum of horizontal displacement on the top of the structure. Considering that this is the spectrum of a nonlinear response signal, the results should only have a qualitative sense. These are shown on Figs. 13 and 14 for the Nice 0.24 g and Nice 0.71 g earthquakes. We observe that such a frequency shift, although far from being negligible, is realistic. It indicates that the
Table 3 Vertical load: crack closure function identification (rf ) (kN) Nice 0.25 g SF 1.13 g Nice 0.40 g Nice 0.71 g
3.5 MPa
1.75 MPa
1.3 MPa
1.0 MPa
Experiment
115 150 119 140
119 160 132 190
120 200 155 240
120 218 165 265
138 198 146 248
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
local degradations is satisfactorily accounted for and permits to follows the global stiffness reduction of the structure. The final crucial point in modeling pertains to computing corresponding local level of degradation for both concrete and steel. Fig. 15 shows a typical result of this kind with the level of damage in concrete and the irreversible strains in steel at the end of the complete loading sequence. One can observe that the global trend observed experimentally is recovered in the computations at this local level. The location of the critical region are positioned on the upper level. This particular behavior is mainly due to the effect of the second pumping eigenmode, whose main effect results in shifting the failure region. Fig. 15 clearly demonstrates that the computed strains always underestimate the experimentally obtained values. The lack of information at local scale (for example real strain in a steel bar at the location of a crack) is the major drawback preventing the designer from expressing physical criteria describing rupture. Due to the small reinforcement ratio, failure phenomenon can only happen by rupture of the steel bars under tension, thus post-peak behavior cannot be represented. 4.3. Simulation of Camus IV test We present in this paper a comparison of the numerical results obtained with the SSI model with those measured in the test Camus IV [3]. 4.3.1. Test characteristics The five storey mock-up (with the same dimensions as one used in Camus I test) is mounted on the shaking table simply supported on a sand container (Fig. 16). The main goal is to study the boundary conditions influence due to the presence of soil and to the possible uplift and sliding. The two shallow foundations (Fig. 17) where dimensioned with a safety factor equal to 3 to avoid too early failure of the soil when uplift will occur. The sand
1235
Fig. 16. Camus IV mock-up.
container dimensions limited by the shacking table capacity where chosen 4 m 4 m, but computations have shown a low influence of the lateral boundaries. As it was impossible to have a realistic depth of sand capable of reproducing the real failure mechanisms of the soil, the goal was limited to find a more or less realistic overturning stiffness of the foundation. The depth was thus chosen equal to 0.4 m. The mock-up was submitted to the Nice synthetic accelerogramm applied horizontally at different levels and followed by the natural signal of Melandy Ranch. The natural frequency of the interacting system tablesand-mock-up was 3.5 Hz which is much lower than for the clamped mock-up of Camus I with frequency of about 7.0 Hz. 4.3.2. Choice of the parameters Although the test conditions are far from the assumptions used for elaboration of the model, by
Fig. 15. Degradation of the structure at the end of the analysis. Location of the cracks on the instrumented wall at the end of the loading sequence. Measured and computed strains (maximum values).
1236
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Fig. 18. Modeling of Camus IV mock-up.
Fig. 17. Details of the foundations.
adapting the parameters (stiffness, failure criterion and damping) to the particular conditions due to the low depth of soil, it is possible to use the macro-element since the yielding of the soil was limited during the test (SF 4; h=B 1:5). The soil acted as an isolator and during the test no damage was observed in the structure. The structure is thus modeled with elastic beams (Fig. 18). Concentrated mass and inertial moments are placed at each storey. The test is modeled in 2D, which is realistic for the structure (since only one wall is modeled) but brings the strong assumption on a shallow foundation. The following fitting procedure for the parameters has been used to take into account the particular shape of the soil in the test (the chosen values of parameters are given in Table 4): – The horizontal and rotational stiffness have been obtained from the very low level response with quasi-linear behavior. – The vertical stiffness has been obtained by calibration of the vertical natural frequency of the soil–structure system. – The radiation damping coefficients should not exist here since the soil is not semi-infinite, nevertheless a small value is kept to take into account the damping in the whole system (since no damping is put in the structure).
Table 4 Chosen parameter values B ¼ 2:1 m, l ¼ 0:8 m, Vmax ¼ 0:8 MN qmax ¼ 0:476 MPa el Kzzel ¼ 320 MN/m, Kxx ¼ 50 MN/m el Khh ¼ 150 MN m/rad el Czzel ¼ 2 MN s/m, Cxx ¼ 2 MN s/m el Chh ¼ 2 MN m s/rd a ¼ 0:7, b ¼ 0:57 c ¼ 1, d ¼ 1, e ¼ 1, f ¼ 1 j ¼ 0:25, n ¼ 0:17
– The bearing capacity was calculated with a safety factor equal to 3 for a rough foundation on a sand layer with a limited depth. Since the sand container is also laterally limited, the safety factor is estimated to 4. We can thus estimate the bearing capacity Vmax . – The failure criterion coefficients have been adapted from the maximum moment measured for high excitation levels of the test and assuming a friction ratio equal to 0.7 at the interface. – Due to the low depth of soil and to its high density, the uplift in the test was greater than in a foundation lying on an infinite soil. An amplification coefficient of the uplift displacements equal to 3 has been used. 4.3.3. Results In Fig. 19a, the S shape hysteresis loops of the moment––rotation response due to the strong uplift is
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
1237
Fig. 19. Comparison between experients and simulation results for Camus IV test.
clearly observed indicating energy dissipation, despite a low depth of soil. The residual settlement (positive) visible in Fig. 19b confirms that the soil yielding is also present.
The shear force and moment are reduced due to the uplift (see Table 5). We see on the other hand a large dynamical variation of the vertical load which can reach as much as half of the static load (Fig. 19c). Namely, a
1238
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Table 5 Comparison between Camus I et Camus IV [3] Storey
Camus I 0.71 g
Camus IV 0.52 g
Camus IV 0.90 g
Camus IV 1.11 g
Fifth floor M (kN m) T (kN)
36.5 40.0
14.3 15.9
25.2 28
26.4 29.4
First floor M (kN m) T (kN)
345 111
136.7 51.4
137.2 56.5
133 58.6
variation of 80 kN is measured while the static load is 200 kN. This is due to the uplift and the impact of the foundation on the soil. Despite it was not developed for this kind of conditions, the model gives good results after some adaptations of the parameters. 5. Conclusions Efficient approaches to finite element analysis in earthquake engineering are presented in this paper. They have been developed for the nonlinear transient analysis of reinforced concrete structures under seismic loading with various boundary conditions including soil–structure interaction. Accurate constitutive relations are used with dedicated simplified structural finite elements (multifiber beam element or soil–structure interaction macro-element). This brings robustness and efficiency while keeping the accuracy. Numerical simulations of two tests of the European research network program CAMUS are presented. The first one focuses on the reinforced concrete behavior while the second concerns the soil–structure interaction influence. The simulation results are in good agreement with the experimental ones. These simplified methods permit to perform parametric studies since they are not overly computer time consuming. We can for instance study the influence of the input motion (type, frequency content, . . .), the influence of the type of design concept of the structure (shear wall with ‘‘plastic hinge’’ versus ‘‘French shear wall’’) or the influence of the soil characteristics.
References [1] AFPS. Guide AFPS 92 Pour la protection des ponts. Presse de lÕENPC. France. 1995. [2] Armstrong PJ, Frederick CO. A mathematical representation of the multiaxial Bauschinger effect. GEGB Report RD/B/N 731. 1966. [3] Combescure D, Chaudat T. Icons European program seismic tests on R/C walls with uplift; Camus IV specimen. ICONS project, CEA/DRN/DMT report, SEMT/EMSI/ RT/00-27/4. 2000.
[4] Cremer C. Modelisation du comportement non lineaire des fondations superficielles sous seisme, PhD thesis, Civil Engng, LMT Cachan, ENS Cachan, France. 2001. [5] Cremer C, Pecker A, Davenne L. Cyclic macro-element of soil–structure interaction: material and geometrical nonlinearities. Int J Num Anal Meth Geomech 2001;25:1257– 84. [6] Cremer C, Pecker A, Davenne L. Modelling of non linear dynamic behavior of a shallow strip foundation with macro-element. Earthquake Eng Struct Dyn, in press. [7] Dube J.F. Modelisation simplifiee et comportement viscoendommageable des structures en beton. PhD thesis, ENS, Cachan. 1994. [8] Gazetas G. Foundations vibrations. In: Fang H-Y, editor. Foundation engineering handbook. NY: van Nostrand Reinhold; 1991 [Chapter 15]. [9] Ibrahimbegovic A, Frey F. Finite element analysis of linear and non-linear planar deformations of elastic initially curved beams. Int J Num Meth Engng 1993;36:3239–58. [10] Ibrahimbegovic A, Wilson EL. A methodology for dynamic analysis of linear structure–foundation systems with local nonlinearities. Earthquake Eng Struct Dyn 1991; 19:1197–208. [11] Ibrahimbegovic A, Wilson EL. A modified method of incompatible modes. Commun Appl Num Meth 1991; 7:187–94. [12] Kotronis P. Cisaillement Dynamique de Murs en Beton Arme. Modeles Simplifies 2D et 3D. PhD thesis, ENS, Cachan. 2000. [13] La Borderie Ch. Phenomenes unilateraux dans un materiau endommageable: modelisation et application a lÕanalyse de structures en beton. PhD thesis, University of Paris VI. 1991. [14] Leger P, Dussault S. Nonlinear seismic response analysis using vector superposition methods. Earthquake Eng Struct Dyn 1991;20. [15] Lemaitre J, Chaboche JL. Mechanics of solids material. Cambridge University Press; 1990. [16] Lowes LN, Moehle JP. Evaluation and retrofit of beamcolumn T-joints in older RC bridge structures. ACI Struct J 1999;96(4):519–32. [17] Mazars J. A description of micro- and macroscale damage of concrete structures. J Eng Fract Mech 1986;25(5/6):729– 37. [18] Newmark NM. A method of computation for structural dynamics. ASCE J Eng Mech Div 1959;85:67–94. [19] Paulay T, Priestley MJN. Seismic design of concrete and masonry buildings. New-York: John Wiley & Sons Inc.; 1992. 744 p. [20] Pecker A. Analytical formulae for the seismic bearing capacity of shallow strip foundations. In: Seco e Pinto PS, editor. Seismic behavior of ground and geotechnical structures. Rotterdam: Balkema; 1997. p. 261–8. [21] Przemieniecky JS. Theory of matrix structural analysis. New-York: MCGraw-Hill; 1985. [22] Ragueneau F, La Borderie Ch, Mazars J. Damage model for concrete like materials coupling cracking and friction, contribution towards structural damping: first uniaxial application. Mech Cohes Frict Mater 2000;5:607–25. [23] Rashid JYR, Dameron RA, Dunham RS. Finite element analysis of reinforced concrete in bridge seismic design
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239 practice ASCE Commitee Report In: Shing PB, Tanabe T, editors. Modeling of inelastic behavior of RC structures under seismic loads. Structural Engineering Institute; 2001. p. 217–33. ISBN 0-7844-0553-0. [24] Taylor RL. FEAP: A finite element analysis program, version 5.01 manual. University of California, Berkeley. 1996.
1239
[25] Ukritchon B, Whittle AJ, Sloan SW. Undrained limit analysis for combined loading of strip footings on clay. J Geotech Geoenv Engng March 1998:265–76. [26] Zienkiewicz OC, Talyor RL. 5th ed. The finite element method, Vols. I, II and III. Oxford: Butterworth Heinemann; 2000.
Efficient approaches to finite element analysis in earthquake engineering L. Davenne a
a,*
, F. Ragueneau a, J. Mazars b, A. Ibrahimbegovic
a
Laboratoire de M ecanique et Technologie, Ecole Normale Sup erieure de Cachan, 61, Avenue du Pr esident Wilson, 94235 Cachan Cedex, France b Laboratoire Sols, Solides, Structures, INPG––BP 53, 38041 Grenoble Cedex 9, France
Abstract This paper deals with the modeling of reinforced concrete structures subjected to earthquake ground motion. Due to the complex behavior of both materials and structures, efficient numerical tools are developed herein in order to keep accuracy and robustness for large scale computations. We focus our attention on the use of simplified multifiber beam element describing the response of structural components and on macro-element accounting for soil–structure interaction. 2003 Elsevier Science Ltd. All rights reserved.
1. Introduction Earthquake engineering design of a civil engineering structure concerns in general several levels of strong motion earthquakes, each one with its specific design requirements. In particular, one should distinguish between: – A damage limitation criterion which allows to keep the operational features of the building intact for a seismic loading which has a large occurrence probability. – A collapse criterion which insures that for the nominal design seismic action the structure remains stable without significant local collapse. – A life safe objective requiring that for the strongest possible earthquake feasible for the site, the global collapse is prevented and the life is preserved. In the traditional approach to earthquake engineering design, the computations are carried out on the basis
* Corresponding author. Tel.: +33-1-4740-2243; fax: +33-14740-2240. E-mail address: [email protected] (L. Davenne).
of linear elastic static analysis. The nonlinear behavior and energy dissipation, can be accounted for in a trivial manner by a force-based approach, where the level of seismic loading is computed by analyzing the elastic response spectra on a single degree of freedom system and the so-called R reduction factor method is then used to introduce the ductility of materials (steel and concrete). A more rational concept developed during this last decade, named displacement-based procedure [19], turns towards the design based on the critical limit states of the structural elements. In other words, defining the ultimate behavior of materials, one has to deal with maximum strains or curvatures, thus taking into account the structure ductility demand in a more physical way. Whereas the force-based design uses the elastic properties of the structure, the displacement-based design characterizes it through the secant stiffness at maximum displacement. The evaluation of vibration frequencies during the loading is more appropriate than for the purely elastic analysis. These classical procedures for design (force or displacement formulation) are based on the study of an equivalent single degree of freedom system or in fact a simplified substitute structure which is not capable to account for the load redistribution inside the structures due to local nonlinearities. This is one of the major drawbacks preventing a realistic description of the
0045-7949/03/$ - see front matter 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0045-7949(03)00038-5
1224
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
global and local behavior of a structure up to failure. One remedy consists in using nonlinear finite element method to perform push-over analysis, i.e. a nonlinear static analysis under monotonic increasing lateral load. It allows to determine the maximum carrying capacity of the structure in terms of forces, displacements, ductility (deformation), crack pattern and failure mode. The main assumption of such a computation is that the response is related to an equivalent single degree of freedom and thus, controlled by a single mode which justifies that the position and amplitude of horizontal lateral loads are typically defined by an elastic analysis. An alternative choice to perform the earthquake resistant design is by making use of the nonlinear timehistory analysis, assuming a physical description of materials and applying transient loadings on the structure in terms of natural or simulated ground motion. The evolution of eigenmodes concomitant to the stiffness degradation governed by local yield criterion provides currently the most refined method of analysis for ultimate behavior of concrete structure [16,23]. For a simple reason of excessive computational costs, such approaches to structural dynamics in civil engineering are not commonly used. Several strategies have been developed in order to lower the CPU time in transient analysis by solving the equations of motion by a mode superposition techniques even in the nonlinear range, either by dynamic substructuring for limited spread of nonlinear zones (e.g. see [10,11]) or using the tangent stiffness to modify the basis [14]. Nonlinear dynamic analysis of complex civil engineering structures based on a detailed finite element model requires large scale computations and handles delicate solution techniques. The necessity to perform parametric studies due to the stochastic characteristic of the input accelerations imposes simplified numerical modeling which will reduce the computation cost. In this work the latter is achieved by selecting the classical Euler–Bernoulli beam model for representing the global behavior of the structural components of a complex civil engineering structure. We note that no simplification is made other than this choice of the beam model. In particular, the constitutive bahavior of the chosen model remains sufficiently general to take into account all the different inelastic phenomena (cracking by damage, irreversible deformation by plasticity and crack-closing by unilateral frictional contact condition). Another effect which is also taken into account concerns the interaction of structure with the foundation, which provides both a more realistic boundary condition as well as the damping model due to radiation effect. The outline of the paper is as follows. In Section 2 we briefly describe the chosen model of multilayer beam used for representing the structural components. The development pertaining to the modeling of structure– foundation interaction are described in Section 3. In
Section 4 we present several numerical simulations which illustrate the performance of the proposed models with respect to the experimental results. Concluding remarks are given in Section 5.
2. Multifibers beam elements for nonlinear transient analysis of reinforced concrete structures In this work we seek to develop the models capable of representing in a reliable manner nonlinear behavior of a structure damaged by a strong earthquake, with the main goal of providing a physically based description of damping (as opposed to a more convenient, but very adhoc Rayleigh damping model). This kind of approach is especially suitable for a novel concept in earthquake resistant design of shear-wall type reinforced concrete structures; referred nowadays as ÔFrench wallsÕ, where the choice of reinforcement is made so that the damaged and cracking spreads as much as possible through the structure which allows that the main part of the energy generated by an earthquake be absorbed through inelastic material behavior. In short we seek to develop the structural models capable of integrating the advanced damage models for concrete (e.g. [17]), in order to describe the nonlinear behavior of a reinforced concrete structures. With respect to the large spreading of the zone with nonlinear behavior we further seek to limit the model complexity (and resulting computational costs) by limiting the diversity of possible deformation global patterns which is achieved by choosing a multifiber beam model, with all the fibers restricted to beam kinematics and with each fiber employing its own constitutive model (see Fig. 1). The main advantage of using a multifiber type finite element concerns the simple uniaxial behavior which allows a very efficient implementation of inelastic constitutive. This is no longer possible for thick beams where shear strains play a major role [7]. 2.1. Cross-section behavior The multifiber beam element developed herein employs the standard Hermite polynomial shape functions to describe the variation of the displacement field along the beam. The difference with ‘‘classical’’ beam elements [21] concerns the cross-section behavior, that is the relation between the generalized strains e and the generalized stresses s. In the general 3D case the latter includes: s ¼ ðN
Mx
My
Mz ÞT
and e ¼ ð e hx
vy
vz ÞT ð1Þ
where N is the normal force, Mx the torque, My and Mz are the bending moments, e the axial strain, hx the twist,
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
1225
Fig. 1. Multifibers beam formulation for reinforced concrete structures.
vy and vz the curvatures. The cross-section behavior is expressed with the constitutive matrix: 2 3 0 C13 C14 C11 6 C22 0 0 7 7 C¼6 ð2Þ 4 C33 C34 5 sym C44 where the coefficients are obtained by integrating over the cross-section (y and z axes): Z Z E dS; C13 ¼ Ez dS; C11 ¼ S S Z Z C14 ¼ Ey dS; C22 ¼ Gðy 2 þ z2 Þ dS S S Z Z C33 ¼ Ez2 dS C34 ¼ Eyz dS S ZS 2 C44 ¼ Ey dS ð3Þ S
where E and G are YoungÕs and shear moduli which vary in y and z. The chosen moduli can be either initial, secant or tangent, depending upon the iterative algorithm used to solve the global equilibrium equations. The component of the constitutive matrix are computed by numerical integrations in (3), with one Gauss point per fiber. For the Euler–Bernoulli element, the shear forces are computed at the element level through the equilibrium equations (included in the Hermite polynomial shape functions). Reinforcement bars are introduced within special fibers, whose behavior is obtained as a combination of those of concrete and steel according to: rlayer ¼ ð1 aÞrconcrete þ arsteel
ð4Þ
where a is the relative area of the reinforcement in the layer. When dealing with structures which slenderness ratio is far from the classical beam such as shear wall, a more reliable representation of shear deformations and shear stresses has to be provided. One possibility in that respect is to use the classical Timoshenko beam model, which can describe the constant shear strain. The main
difficulty of developing the finite element implementation of the Timoshenko beam model concerns the socalled shear locking phenomena (e.g. [26]), or inability of the standard finite element approximations to represent pure bending vanishing shear modes. A number of different remedies to shear locking problem has been proposed (e.g. see [9]), ranging from selective or reduced integration, assumed shear strain, enhanced shear strain or hierarchical displacement interpolations. A very recent work of Kotronis [12] extends these ideas in order to construct shear locking remedies for a multifiber Timoshenko beam. Remark: In order to retain the simplicity of the constitutive model for the Timoshenko formulation, the normal and shear stress resultants are uncoupled and nonlinearities are only introduced in the stress–strain relationship in the direction of the normal to the crosssection, whereas the tangential behavior is assumed to remain elastic. A more complex constitutive relationships can be used for each fiber, expressed within the framework of irreversible processes thermodynamics which allows physical coupling of all stress components. Such a development for 2D Timoshenko multilayer beam was carried out by Dube [7], where the shear stress as well as the normal stress in each layer was used to compute the damage of the layer with the 2D version of the damage model presented hereafter. 2.2. Enhanced strain formulation For a multifiber beam model, as indicated in (2), the nonlinear behavior implies coupling between axial force N and bending moment M. It is thus necessary that the axial force and bending moment have the same variation along the element. The latter is not satisfied for an Euler–Bernoulli multifiber beam element, where axial strain remains constant along the element whereas the bending strain is linear. This problem may become important in the nonlinear range where a linear curvature along the beam implies a linear variation of the axial strain due to the shift of the neutral axis.
1226
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
x x u2 u1 u1 þ u2 ) eðxÞ ¼ ¼ cst uðxÞ ¼ 1 L L L
x2 x3 x3 x2 2 þ x h1 vðxÞ ¼ 1 3 2 þ 2 3 v1 þ L L L2 L
3
3 2 x x x x2 h2 ) jðxÞ þ 2 3 þ 3 2 v2 þ L L L2 L
6 12x 6x 4 ¼ 2 þ 3 v1 þ h1 L L L2 L
12x 6 6x 4 h2 þ 3 þ 2 v2 þ L L L2 L
u_ tþDt ¼
ð5Þ One way to remove this incompatibility of variations for axial strain and curvature is by using enriched description of the axial strain field with: eðxÞ ) eðxÞ þ ~eðxÞ;
~eðxÞ ¼ GðxÞa;
4 8x GðxÞ ¼ 2 L L ð6Þ
where GðxÞ is an enhanced axial strain which can be derived from the displacement bubble function. The variational basis of such an enhancement is provided by Hu–Washizu principle (e.g. [26]), which can be presented in the same manner as the modified method of incompatible modes of Ibrahimbegovic and Wilson [10,11] Z Z 0 ¼ ðe ðxÞN ðxÞ þ j ðxÞMðxÞÞ dx ¼ u ðxÞf ðxÞ dx L L Z ~e ðxÞN ðxÞ dx ð7Þ 0¼ L
2.3. Global solution procedure Dynamic analysis for earthquake ground motion for this kind of structural model reduces to solving the set of nonlinear equations which can be written as: M€ uðtÞ þ C u_ ðtÞ þ f int ðu; tÞ ¼ Mr€ug ðtÞ
ð8Þ
where M and C are, respectively, mass and damping matrix, €uðtÞ and u_ ðtÞ are nodal accelerations and velocities, f int ðu; tÞ is the internal force vector obtained by discretizing (6) and r€ug ðtÞ is the ground acceleration effect applied on the structure. By using a time-integration scheme, the differential equation of motion in (8) is reduced to an algebraic equation. In particular for low frequency response in earthquake engineering analysis one uses an implicit scheme such as Newmark one-step scheme (see [18]) which allows to express the velocity and acceleration vectors at time t þ Dt as the functions of their corresponding values at time t and incremental displacement vector: € utþDt ¼
1 1 1 2b €uðtÞ u_ ðtÞ Du b Dt2 b Dt 2c
ð9Þ
c c 1 2b Du þ 1 u_ ðtÞ u_ ðtÞ b Dt b 2c
c € þ Dt 1 uðtÞ 2b
ð10Þ
where c ¼ 1=2 and b ¼ 1=4 are typically chosen for optimal result accuracy. The discrete set of equations, obtained by introducing (9) and (10) into (8), is further solved by an iterative solution procedure, where at each iteration the problem reduces to
1 c C þ K Dui ¼ Mr€ M þ ug ðt þ tÞ b Dt2 b Dt M€ ui ðt þ DtÞ C u_ i ðt þ DtÞ f int ðui ðt þ DtÞÞ uiþ1 ðt þ DtÞ ¼ ui ðt þ DtÞ þ Dui
ð11Þ
Different iterative strategies can be used to compute the tangent stiffness matrix K in (11) above; if we seek to reduce the computational cost one would not update K at each iteration as for NewtonÕs method, but either keep it constant as for modified NewtonÕs method or use a secant stiffness for quasi-Newton method. The details of computing K and f int pertain to the constitutive model which is specified for each component as explained subsequently. 2.4. Materials constitutive relations for dynamics Both steel and concrete materials are described within the thermodynamic framework for irreversible processes (e.g. see [15]). In describing the nonlinear behavior of reinforcement bars, we choose the classical plasticity model and take into account the nonlinear kinematic hardening of Armstrong and Frederick [2] in order to be able to better describe the observed hysteresis loops. The free energy for this model can be written as: qw ¼ 12ðe ep Þ : H : ðe ep Þ þ 12ba : a
ð12Þ
in which H is the HookeÕs elasticity tensor, ep is the plastic strain and a the hardening internal variable which describes hardening. The constitutive equations for this kind of model can be written as: r¼
oðqwÞ ¼ C : ðe ep Þ; oe
X ¼
oðqwÞ ¼ ba oa
ð13Þ
where X is the stress-like hardening variable. The latter is used to describe a modified form of the plasticity criterion: f ¼ J2 ðr X Þ þ 34aX : X ry 6 0
ð14Þ
where a, b and ry are material parameters. Due to the particular geometric characteristics of steel bars, only a 1D implementation of the model is carried out. A typical
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
1227
Fig. 2. Hysteresis loop for the steel modeling.
stress–strain hysteresis loops predicted by this model is given in Fig. 2. The constitutive model for concrete for earthquake engineering ought to take into account some observed phenomena, such as decrease in material stiffness due to cracking, stiffness recovery which occurs at crack closure and inelastic strains concomitant to damage. To simulate this behavior we employ the damage model of La Borderie [13], which uses two scalar damage variables, D1 for damage in tension and D2 for damage in compression. Unilateral effect and stiffness recovery (damage deactivation) are also included. Inelastic strains are taken into account thanks to an isotropic tensor. The Gibbs free energy of this model can be expressed as: v¼
hriþ : hriþ hri : hri t þ þ ðr : r Trðr2 ÞÞ 2Eð1 D1 Þ 2Eð1 D2 Þ E b 1 D1 b2 D2 f ðrÞ þ TrðrÞ þ Eð1 D1 Þ Eð1 D2 Þ
ð15Þ
where b1 and b2 are material coefficients, whereas h iþ and h i denotes the positive or negative values of the given variable. f ðrÞ and rf are the crack closure function and the crack closure stress respectively. E is the initial YoungÕs modulus and m the Poisson ratio. For the 1D expression, the total strain is: e ¼ ee þ ein rþ r þ Eð1 D1 Þ Eð1 D2 Þ b 1 D1 b2 D2 f ðrÞ þ ein ¼ Eð1 D1 Þ Eð1 D2 Þ ee ¼
ð16Þ
where ee the elastic strain and ein the inelastic strain. The partition of the stress is obtained as:
r > 0 ! rþ ¼ r; r < 0 ! rþ ¼ 0;
r ¼ 0 r ¼ r
and for the crack-closure function:
ð17Þ
Fig. 3. Unidimensional constitutive relations for concrete.
8 < r > 0 ! f ðrÞ ¼ 1 rf < r < 0 ! f ðrÞ ¼ 1 rrf : r < rf ! f ðrÞ ¼ 0
ð18Þ
where rf is the stress where the cracks are completely reclosed. The role of the crack closure function is: – To restore the material stiffness, passing from tensile to compressive behavior. – To null the tensile inelastic strains creating after tension when excessive (over rf ) compression is reached. Damage criteria are expressed as fi ¼ Yi Zi (i ¼ 1 for tension or 2 for compression, Yi is the associated force to the damage variable Di and Zi a threshold dependent on the hardening variables). The evolution laws for the damage variables Di are written as Di ¼ 1
1 1 þ ½Ai ðYi Y0i Þ Bi
ð19Þ
where Y0i is the initial elastic threshold and Ai , Bi are material constants. A typical stress–strain response produced by this model for a uniaxial cyclic loading which passes through traction, compression and again traction is given in Fig. 3. The material coefficients Y0i , Ai ; Bi , bi (i ¼ 1, 2) and rf can be identified from uniaxial tests. They are only function of intrinsic material characteristics (aggregate size, concrete mix, . . .).
3. A soil–structure interaction macro-element It has already been noted that, after severe earthquakes, some structures with different properties from
1228
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
surrounding ground have been preserved from any damage. The interaction of the structure with the soil can lead to a less severe response of the structure either by modifying the support condition with respect to a clamped boundary or accounting for modification of the properties of the interacting system. In particular the low stiffness of the soil reduces the system frequency, and may, depending on the frequency content of the input motion, drive outside of the resonance zone of the input spectrum. Moreover, the nonlinearities in the surrounding soil and interface (yielding of the soil and uplift of the structure) lead to a strong energy dissipation and further reduce the forces in the structure. To model these phenomena we developed a simplified method for soil–structure interaction. The overall linear and nonlinear behavior in the soil and at the interface is reduced through a macro-element expressed in three degrees of freedom on the strip foundation. The constitutive laws governing the behavior of this macroelement are developed within the framework of the plasticity theory. The quasi-static cyclic loading model is developed first (see [4,5] followed by an extension to dynamic loading case [6]. The only variables of the macro-element are the forces at the base of the foundation (vertical force V , horizontal force H and moment M) and the corresponding displacements (vertical displacement z, horizontal displacement x and rotation h) measured at the centre of the foundation (assumed as a rigid body). These variables normalized by characteristic values to become dimensionless.
0 1 1 V V0 1 B C B C F @ H0 A @ H A and Bqmax 0 M=B M 0 01 0 1 z z B C 1B C u @ x0 A @ x A B h0 Bh 0
ð20Þ
where B is the width and Vmax ¼ Bqmax the bearing capacity of the foundation. Since it is a strip foundation, all terms are expressed for a unit transverse length (unit depth of the footing). Three different phenomena are represented by the forces acting on the foundation: the plasticity of the soil surrounding the foundation or near field, the uplift of the foundation and the elasticity of the soil in far field. They are modeled through three interacting models (Fig. 4). Plasticity of the soil and uplift are strongly nonlinear and coupled since they both concern the near field while elasticity concerns the far field. The total displacement increment is partitioned into elastic, plastic and uplift components, which can be written choosing the corresponding rate decomposition: u_ tot ¼ u_ el þ u_ pl þ u_ up
ð21Þ
The global model of constitutive behavior can be written: F_ ¼ ð½Kelpl 1 þ ½Kup 1 Þ1 : u_ tot
ð22Þ
where elastoplastic tangent stiffness matrix Kelpl is obtained from the plasticity model and the tangent stiffness matrix Kup from the chosen ‘‘uplift model’’.
Fig. 4. Idealization of the complete model.
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
1229
3.1. Elasticity and damping
3.2. Plasticity model
The energy dissipation mechanism in the near field arises exclusively from the nonlinearities included in the plastic and uplift cyclic models. The energy dissipation in the far field arises from the radiation boundary conditions in semi-infinite medium. The latter is modeled through the imaginary part of the elastic impedances [8], represented by a dashpot, and the real part is the elastic stiffness matrix. 0 0 1 Kzz 0 0 0 0 A Kel ¼ @ 0 Kxx 0 0 0 Khh 1 0 K zz 0 0 C B qmax C B Kxx B 0 C ð23Þ ¼B 0 C C B qmax A @ Khh 0 0 B2 qmax
Finite element analysis with an efficient model for the soil [4] or analytical results, when available, guided us in constructing the macro-model presented in the following. The elastoplastic model provide the elastic and plastic parts of the displacement increments:
0:73 2 2 G0 ð1 þ 2aÞ; Kxx ¼ G0 1 þ a ; 1m 2m 3
2
p B 1 Khh ¼ G0 ð24Þ 1þ a 2ð1 mÞ 2 3 Kzz ¼
where a is the gradient of the variation of the shear modulus G ¼ G0 ð1 þ afÞ (G0 the shear modulus at z ¼ 0) with respect to the depth f ¼ 2z=B and m the PoissonÕs ratio. In the same way, the coefficient of the damping matrix are [8]: 3:4 Ab ; pð1 mÞ 3:4 Ibh Chh ¼ qVs0 pð1 mÞ
Czz ¼ qVs0
F_ ¼ Kelpl : ðu_ el þ u_ pl Þ with
Q¼
ofi ; oF
ð27Þ
where fi is the yield surface and g the plastic potential for nonassociated flow rule employed herein. One also defines a failure surface from analytical [20] and numerical results [25], based on limit analysis for structure overturning with uplift mechanisms (the uplift influence is also included in the plasticity model). !2 !2 H0 M0 fi ¼ þ 1¼0 aV 0c ð1 V 0 Þd bV 0e ð1 V 0 Þf ð28Þ For a soil with constant gradient of cohesion with depth, the coefficients in the last expression are given as: a ¼ 0:32=x; e ¼ 0:8;
where q is the mass density of the soil, Ab the areapof theffi ffiffiffiffiffiffiffiffiffiffi foundation, Ibh its quadratic moment and Vs0 ¼ G0 =q the shear velocity at z ¼ 0. Remark: Step by step analysis prevents from including the frequency variation of the dynamic impedances. Constant values are thus chosen (static value for stiffness, which are already given for the model, and the value near resonance of the soil structure system for damping, which gives the better average damping). This choice is acceptable for a semi-infinite soil medium with a slow and continuous properties variation (nearly homogeneous soil). In this case the variation of the coefficients is low over the frequency range of interest (0–10 Hz). For a foundation lying on a thin soil layer, this approximation is no longer valid. In this case one could take the values corresponding to the fundamental frequency of the soil building system.
1 ðKel : PÞ ðQ : Kel Þ; h0 þ h0 og P¼ oF
Kelpl ¼ Kel
Cxx ¼ qVs0 Ab ; ð25Þ
ð26Þ
b ¼ 0:37=x0:2 ;
c ¼ 0:25; qmax f ¼ 0:8 with x ¼ q0 max
d ¼ 0:55;
where qmax is the actual ultimate bearing capacity of the foundation (soil with cohesion gradient c ¼ c0 þ gz) and q0 max the ultimate bearing capacity for an homogeneous soil with a constant cohesion c0 . The model is particularly well suited for the evolution of the yield surface towards the failure surface with a vertical force being almost constant, for reproducing the behavior of the soil and the foundation, initially submitted to the weight of the structure, and then loaded mainly in the (H ; M) plane (seismic loading).
fi ¼
H 0 aCcH qCcH
2
þ
M 0 bCcM qCcM
2 1¼0
ð29Þ
with CcH ¼ aV 0c ðc V 0 Þd ; CcM ¼ aV 0e ðc V 0 Þf ; c ¼ v þ ð1 vÞðq þ sÞ; v ¼ N =Vmax and N is the weight of the structure (initialisation under dead load). In (29) above s ð 0 a b ÞTqisffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the kinematic hardening vector with s ¼ ðs : sÞ1=2 ¼
parameter.
a2 þ b2 and is the isotropic hardening
1230
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
The consistency rule, the criterion of noninterpenetration of the failure surface with the loading surface and a relationship obtained from the observation of the foundation behavior during uplift allow us to define the evolution of the isotropic and kinematic hardening parameters. The kinematic hardening increment is defined by its amplitude s_ and its direction l: s_ ¼ s_
l jlj
with s_ ¼ hk_i ofi os
h0 l
i : jlj þ of oq
and
1 k_ ¼ 0 Q : F_ h
0
oz0up B oV 0 B F_ ¼ Kup : u_ up ¼ ðCup Þ1 : u_ up ; Cup ¼ B 0 @ 0up oh oV 0
The displacements induced by uplift are given as:
z0up 1 0 1 0 C cos w H B H C B C c C and w ¼ a tan C lB H B 1 C 0 @ CM sin w M A
ð35Þ C1H M 0 bC1M
C1M H 0 aC1H
!
CcM
ð31Þ The isotropic hardening increment is set equal to the kinematic one: q_ ¼ s_ . The plastic modulus h0 is fitted on the M 0 –h0 curves since it is the most important behavior of the foundation
K el ðdM 0 =dh0 Þ2M 0 of dV 0 of dH 0 of h0 ¼ elhh þ þ Khh ðdM 0 =dh0 Þn2 oV 0 dM 0 oH 0 dM 0 oM 0 ð32Þ K 0 h0 hh0 M1
dM 0 0 ¼ Khh exp dh0
!
The external normal to the yield surface is not convenient for defining the direction of the plastic displacements and a nonassociated flow rule is required. The plastic potential g has been chosen as:
g¼
H j
0 2
þ
M n
0 2
þ V 02 1 ¼ 0
ðd=dmax Þ2 and 1 d=dmax
jh0 j d=dmax þ lnð1 d=dmax Þ ¼ ð1 V 0 Þ 2 1 d=dmax
h0up ¼ ð1 V 0 Þh0
ð30Þ where 0
1 oz0up oM 0 C C 0 0 C 0up A oh 0 oM 0 ð34Þ 0
ð33Þ
3.3. Foundation uplift model The foundation uplift behavior on an elasto-plastic soil is a very complex phenomenon which depends strongly on the level of soil yielding. A cyclic model which takes into account the irreversibility of the uplift behavior due to plasticity is developed as a direct modification of the uplift model of Cremer [4] developed for elastic soils. The uplift model operates only on the uplift part of the total displacement introducing coupling between the vertical displacement zup and the rotation hup of the foundation:
0ð0Þ M0 j
0 is the uplift reached where dmax ¼ ð4=V 0 ÞjMmax for the maximum moment at failure for the given V 0 and 0 for H 0 ¼ 0, that is to say for Mmax ¼ bV 0e ð1 V 0 Þf . 0 For a given V , the uplift ratio d is obtained as follows: 0ðiÞ
• if jM 0 j 6 jM0 j, d ¼ 0 0ð0Þ 0ð0Þ • if jM 0 j P jMp0ðiÞ j, d ¼ ð4=V 0 ÞjM 0 M0 j with M0 ¼ 0 ðV 0 =4Þ expAV 0ðiÞ 0ðiÞ • if jM0 j < jM 0 j < jMp0ðiÞ j, d ¼ ð4=gV 0 ÞjM 0 M0 j 0ðiÞ
0ðiÞ
0 0ðiÞ where M0 is such that dðiÞ p ¼ ð4=gV ÞjMp M0 j ¼ 0ð0Þ ð4=V 0 ÞjMp0ðiÞ M0 j with Mp0ðiÞ the maximum moment reached during the loading history and dðiÞ p is the corresponding uplift. The right and left uplift are treated independently supposing they remain uncoupled. Note: The macro-element treats uplift in a continuous global model. No coefficient of restitution to describe the energy dissipation upon impact while the uplifted foundation will contact the soil has been introduced.
3.4. Use of the model This SSI macro-element which was implemented in computer code FEAP [24] is well suited for the dynamic simulation of slender structures with potentially large uplift and moderate yielding in the soil (and no sliding of the footing on the soil). When using this model it is recommended to have a safety factor SF ¼ Vmax =V > 3 (where Vmax is the bearing capacity of the foundation, V the vertical load on the foundation), and a slenderness ratio h=B > 1 (where h is the height of the dynamic centre of gravity and B the width of the foundation).
4. Numerical simulations and experimental results (CAMUS program) In this section we briefly present the main results of the experimental program carried out within the re-
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
search network CAMUS along with the comparative numerical simulations. The main purpose of this experimental program consists in demonstrating the ability of lightly reinforced concrete bearing walls to sustain very severe earthquake. The design is based on a concept favoring damage spreading over several storeys of a lightly reinforced concrete wall. This kind of design leads to a lower percentage of reinforcements with an optimized distribution generating a wide cracks pattern allowing for dissipation of great amounts of energy. As a consequence, one obtains a vertical rising of the masses resulting in energy transformation (from kinematic to potential instead of storing strain energy in the structure) and an increasing ductility thanks to this particular way of dissipating the earthquake energy. In order to test this concept experimentally, a one-third scaled model (see Fig. 5) has been installed on the shaking table of CEA, composed of two parallel braced walls linked by six square slabs. A massive reinforced concrete footing allows the anchorage to the shaking table. The similarity laws to a full-scale structure and the mock-up require additional masses of 6.55 t at each storey. The problem is assumed to be two-dimensional since the mock-up is loaded through horizontal acceleration parallel to the walls and the presence of steel bracing systems at each storey placed perpendicularly to the loading direction prevent occurrence of any torsional modes. pffiffiffiThe accelerograms are scaled in time with a ratio of 1= 3 in order to account for the similarity rules. Two types of accelerograms are applied: a synthetic one generated from the french code [1] Nice S1 as the far field type earthquake and the San Francisco earthquake of 1957 (Fig. 6) as the near field one. The response spectra show clearly the main difference of these two earthquakes: on one hand Nice has a broad frequency content and on the other hand, San Francisco has a
1231
narrower effective band width. The complete experimental sequence was: Nice 0.25 g, San Francisco 1.13 g, Nice 0.4 g and Nice 0.71 g. 4.1. Model calibration (modal analysis and damping calibration) Eigenfrequencies of the reduced scale numerical model anchored to the shaking table show a higher global stiffness. Namely the frequency of the first eigenmode measured before the test is 7.3 Hz, as opposed to 10.3 Hz predicted by the numerical model. This indicates the necessity to take into account the shaking table compliance as well as the anchorage system which can never be considered as perfectly clamped. The second eigenmode corresponds to pumping with an eigenfrequency of 40 Hz as predicted by the finite element model. During the test, the vertical displacement of the mass induced by the openings of cracks excited this second mode. The resulting vibration measured through the vertical acceleration of the shaking table, with a frequency of 20 Hz leads to variation of the vertical dynamic forces which can be very important for reinforced concrete structure failure. This effect must clearly be taken into account through the vertical stiffness of the connecting rods. A simplified model for the basement footing is assumed consisting in an horizontal beam whose properties are selected in accordance with the shaking table stiffness with Kv ¼ 48EI=l3 and Kh ¼ 12EI=l. This kind of footing element adds two eigenmodes, as required for the CAMUS structure, since as already discussed, the second vertical mode plays a very important role. Finally, YoungÕs modulus of the footing beam was set to 10 000 MPa in order to introduce anchorage and contact defects. Table 1 summarizes the correction of fundamental frequencies of the system due to these new boundary conditions.
Fig. 5. CAMUS structure.
1232
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Fig. 6. Nice and San-Francisco accelerograms and response spectrum.
4.2. Comparisons between experimental and numerical results
Table 1 Eigenmode adjustment
First mode Second mode
Initial modeling (Hz)
Adjusted modeling (Hz)
Measured frequencies (Hz)
10.3 40.0
7.4 19.0
7.3 20.0
Damping is another important concern. In particular, the classical viscous damping is assumed to account for various sources of energy dissipation not related to concrete cracking [22]. Despite the lack of clear physical interpretation, the damping is often introduced in the analysis through the Rayleigh damping matrix which is expressed as a linear combination of the mass matrix and the initial stiffness matrix: C ¼ aM þ bK The particular values of damping coefficients are set to ensure a relative damping value of 1% in the first mode and 2% in the second mode. Great care is taken to keep these damping values as stable as possible during all the analysis, despite the fact that cracking induces reduction of the stiffness and shift of the fundamental frequency.
The complete experimental sequence of loading (four earthquakes) was computed using the proposed numerical model, with the material parameters following: YoungÕs modulus of concrete Ec ¼ 30 000 MPa for a maximum compressive strength of fc ¼ 35 MPa and tensile strength ft ¼ 3 MPa, YoungÕs modulus of steel Es ¼ 200 000 MPa with yield stress ry ¼ 414 MPa accounting for hardening with ultimate stress ru ¼ 480 MPa. The numerical and experimental results are compared in terms of time-history displacements, ductility and local damage criterion (Figs. 7–9). For the lowest level of loading (Nice 0.24 g), when the concrete cracks in tension and the reinforcement bars reach the yield stress, the global behavior of the structure is well represented by the numerical model. The loss in stiffness during the later stages of loading and the decrease of the fundamental frequency is also predicted in a fairly acceptable manner, with a maximum displacement always slightly underestimated before and overestimated after the maximum of the load. Table 2 summarizes the comparisons of experimental and computed results for reaction forces for different
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Fig. 7. Nice 0.24 g: experiment/computation comparisons.
Fig. 8. San-Francisco 1.13 g: experiment/computation comparisons.
1233
accelerograms which indicates globally a very good agreement. Another feature of material modeling during cyclic loading, which is barely taken into account, is the cracks closing there described by the crack closure function of the constitutive relation. When re-closing the cracks, the point where the compressive stiffness is totally recovered depends on the crack closure stress of the model rf (see Fig. 10). In the CAMUS experiment, the shocks induced by cracks closure which activates the vertical eigenmode leads to variation of vertical dynamic forces. Experimentally, the phenomenon can be quantified by the measurement of an induced vertical acceleration at the shaking table. The discontinuity of the stiffness when rf is reached reproduces the shock, two computations performed with different forms of the crack closure function (Figs. 11 and 12) allow to reproduce this experimental observation. Modeling this kind of feature become very important for reinforced concrete structures where the design takes into account the coupling between flexural and axial loading. Through this short analysis we can easily remark that the local behavior of concrete plays a very important role on the global structural response and that the material parameters of the crack closure functions can be identified using different levels of loading (see Table 3). Due to cracking of the concrete, the global stiffness of the structure decreases. This variation does not remain
Fig. 10. Effect of the crack closure stress on the slope of the stiffness recovery during cyclic loading.
Fig. 9. Nice 0.71 g: experiment/computation comparisons.
Table 2 Global response comparisons
Nice 0.24 g SF 1.1 g Nice 0.4 g Nice 0.7 g
Displacement (cm)
Shear load (kN)
Moment (kN m)
Vertical l (kN)
Exp.
Comp.
Exp.
Comp.
Exp.
Comp.
Exp.
Comp.
0.72 1.2 1.35 4.4
0.61 1.1 1.1 3.9
65.9 106 86.6 111
65 90 75 120
200 280 280 345
200 240 240 380
202 271 217 312
190 270 225 310
1234
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Fig. 11. Initial model: vertical load variation. Fig. 13. Nice 0.24 g: displacement spectrum.
1.5
Fig. 12. Modified model: vertical load variation.
constant during the loading but depends on the cracks opening situation (opened or closed crack). This stiffness decrease shifts greatly the natural frequencies. The analysis of the fundamental frequency of the structure subject to a white noise after each sequence is not satisfactory because it is applied on the cracked structure but for closed cracks (this is due to the low level of the white noise which does not allow to open cracks). Such experimental measures cannot capture this kind (changes in stiffness reduction during the earthquake) of structural behavior and provide the first eigenfrequency of 6 Hz instead of the 3 Hz., during the strongest part of
Fig. 14. Nice 0.71 g: displacement spectrum.
the earthquake motion. One way to estimate the dominant frequency during the loading is to plot the response spectrum of horizontal displacement on the top of the structure. Considering that this is the spectrum of a nonlinear response signal, the results should only have a qualitative sense. These are shown on Figs. 13 and 14 for the Nice 0.24 g and Nice 0.71 g earthquakes. We observe that such a frequency shift, although far from being negligible, is realistic. It indicates that the
Table 3 Vertical load: crack closure function identification (rf ) (kN) Nice 0.25 g SF 1.13 g Nice 0.40 g Nice 0.71 g
3.5 MPa
1.75 MPa
1.3 MPa
1.0 MPa
Experiment
115 150 119 140
119 160 132 190
120 200 155 240
120 218 165 265
138 198 146 248
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
local degradations is satisfactorily accounted for and permits to follows the global stiffness reduction of the structure. The final crucial point in modeling pertains to computing corresponding local level of degradation for both concrete and steel. Fig. 15 shows a typical result of this kind with the level of damage in concrete and the irreversible strains in steel at the end of the complete loading sequence. One can observe that the global trend observed experimentally is recovered in the computations at this local level. The location of the critical region are positioned on the upper level. This particular behavior is mainly due to the effect of the second pumping eigenmode, whose main effect results in shifting the failure region. Fig. 15 clearly demonstrates that the computed strains always underestimate the experimentally obtained values. The lack of information at local scale (for example real strain in a steel bar at the location of a crack) is the major drawback preventing the designer from expressing physical criteria describing rupture. Due to the small reinforcement ratio, failure phenomenon can only happen by rupture of the steel bars under tension, thus post-peak behavior cannot be represented. 4.3. Simulation of Camus IV test We present in this paper a comparison of the numerical results obtained with the SSI model with those measured in the test Camus IV [3]. 4.3.1. Test characteristics The five storey mock-up (with the same dimensions as one used in Camus I test) is mounted on the shaking table simply supported on a sand container (Fig. 16). The main goal is to study the boundary conditions influence due to the presence of soil and to the possible uplift and sliding. The two shallow foundations (Fig. 17) where dimensioned with a safety factor equal to 3 to avoid too early failure of the soil when uplift will occur. The sand
1235
Fig. 16. Camus IV mock-up.
container dimensions limited by the shacking table capacity where chosen 4 m 4 m, but computations have shown a low influence of the lateral boundaries. As it was impossible to have a realistic depth of sand capable of reproducing the real failure mechanisms of the soil, the goal was limited to find a more or less realistic overturning stiffness of the foundation. The depth was thus chosen equal to 0.4 m. The mock-up was submitted to the Nice synthetic accelerogramm applied horizontally at different levels and followed by the natural signal of Melandy Ranch. The natural frequency of the interacting system tablesand-mock-up was 3.5 Hz which is much lower than for the clamped mock-up of Camus I with frequency of about 7.0 Hz. 4.3.2. Choice of the parameters Although the test conditions are far from the assumptions used for elaboration of the model, by
Fig. 15. Degradation of the structure at the end of the analysis. Location of the cracks on the instrumented wall at the end of the loading sequence. Measured and computed strains (maximum values).
1236
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Fig. 18. Modeling of Camus IV mock-up.
Fig. 17. Details of the foundations.
adapting the parameters (stiffness, failure criterion and damping) to the particular conditions due to the low depth of soil, it is possible to use the macro-element since the yielding of the soil was limited during the test (SF 4; h=B 1:5). The soil acted as an isolator and during the test no damage was observed in the structure. The structure is thus modeled with elastic beams (Fig. 18). Concentrated mass and inertial moments are placed at each storey. The test is modeled in 2D, which is realistic for the structure (since only one wall is modeled) but brings the strong assumption on a shallow foundation. The following fitting procedure for the parameters has been used to take into account the particular shape of the soil in the test (the chosen values of parameters are given in Table 4): – The horizontal and rotational stiffness have been obtained from the very low level response with quasi-linear behavior. – The vertical stiffness has been obtained by calibration of the vertical natural frequency of the soil–structure system. – The radiation damping coefficients should not exist here since the soil is not semi-infinite, nevertheless a small value is kept to take into account the damping in the whole system (since no damping is put in the structure).
Table 4 Chosen parameter values B ¼ 2:1 m, l ¼ 0:8 m, Vmax ¼ 0:8 MN qmax ¼ 0:476 MPa el Kzzel ¼ 320 MN/m, Kxx ¼ 50 MN/m el Khh ¼ 150 MN m/rad el Czzel ¼ 2 MN s/m, Cxx ¼ 2 MN s/m el Chh ¼ 2 MN m s/rd a ¼ 0:7, b ¼ 0:57 c ¼ 1, d ¼ 1, e ¼ 1, f ¼ 1 j ¼ 0:25, n ¼ 0:17
– The bearing capacity was calculated with a safety factor equal to 3 for a rough foundation on a sand layer with a limited depth. Since the sand container is also laterally limited, the safety factor is estimated to 4. We can thus estimate the bearing capacity Vmax . – The failure criterion coefficients have been adapted from the maximum moment measured for high excitation levels of the test and assuming a friction ratio equal to 0.7 at the interface. – Due to the low depth of soil and to its high density, the uplift in the test was greater than in a foundation lying on an infinite soil. An amplification coefficient of the uplift displacements equal to 3 has been used. 4.3.3. Results In Fig. 19a, the S shape hysteresis loops of the moment––rotation response due to the strong uplift is
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
1237
Fig. 19. Comparison between experients and simulation results for Camus IV test.
clearly observed indicating energy dissipation, despite a low depth of soil. The residual settlement (positive) visible in Fig. 19b confirms that the soil yielding is also present.
The shear force and moment are reduced due to the uplift (see Table 5). We see on the other hand a large dynamical variation of the vertical load which can reach as much as half of the static load (Fig. 19c). Namely, a
1238
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239
Table 5 Comparison between Camus I et Camus IV [3] Storey
Camus I 0.71 g
Camus IV 0.52 g
Camus IV 0.90 g
Camus IV 1.11 g
Fifth floor M (kN m) T (kN)
36.5 40.0
14.3 15.9
25.2 28
26.4 29.4
First floor M (kN m) T (kN)
345 111
136.7 51.4
137.2 56.5
133 58.6
variation of 80 kN is measured while the static load is 200 kN. This is due to the uplift and the impact of the foundation on the soil. Despite it was not developed for this kind of conditions, the model gives good results after some adaptations of the parameters. 5. Conclusions Efficient approaches to finite element analysis in earthquake engineering are presented in this paper. They have been developed for the nonlinear transient analysis of reinforced concrete structures under seismic loading with various boundary conditions including soil–structure interaction. Accurate constitutive relations are used with dedicated simplified structural finite elements (multifiber beam element or soil–structure interaction macro-element). This brings robustness and efficiency while keeping the accuracy. Numerical simulations of two tests of the European research network program CAMUS are presented. The first one focuses on the reinforced concrete behavior while the second concerns the soil–structure interaction influence. The simulation results are in good agreement with the experimental ones. These simplified methods permit to perform parametric studies since they are not overly computer time consuming. We can for instance study the influence of the input motion (type, frequency content, . . .), the influence of the type of design concept of the structure (shear wall with ‘‘plastic hinge’’ versus ‘‘French shear wall’’) or the influence of the soil characteristics.
References [1] AFPS. Guide AFPS 92 Pour la protection des ponts. Presse de lÕENPC. France. 1995. [2] Armstrong PJ, Frederick CO. A mathematical representation of the multiaxial Bauschinger effect. GEGB Report RD/B/N 731. 1966. [3] Combescure D, Chaudat T. Icons European program seismic tests on R/C walls with uplift; Camus IV specimen. ICONS project, CEA/DRN/DMT report, SEMT/EMSI/ RT/00-27/4. 2000.
[4] Cremer C. Modelisation du comportement non lineaire des fondations superficielles sous seisme, PhD thesis, Civil Engng, LMT Cachan, ENS Cachan, France. 2001. [5] Cremer C, Pecker A, Davenne L. Cyclic macro-element of soil–structure interaction: material and geometrical nonlinearities. Int J Num Anal Meth Geomech 2001;25:1257– 84. [6] Cremer C, Pecker A, Davenne L. Modelling of non linear dynamic behavior of a shallow strip foundation with macro-element. Earthquake Eng Struct Dyn, in press. [7] Dube J.F. Modelisation simplifiee et comportement viscoendommageable des structures en beton. PhD thesis, ENS, Cachan. 1994. [8] Gazetas G. Foundations vibrations. In: Fang H-Y, editor. Foundation engineering handbook. NY: van Nostrand Reinhold; 1991 [Chapter 15]. [9] Ibrahimbegovic A, Frey F. Finite element analysis of linear and non-linear planar deformations of elastic initially curved beams. Int J Num Meth Engng 1993;36:3239–58. [10] Ibrahimbegovic A, Wilson EL. A methodology for dynamic analysis of linear structure–foundation systems with local nonlinearities. Earthquake Eng Struct Dyn 1991; 19:1197–208. [11] Ibrahimbegovic A, Wilson EL. A modified method of incompatible modes. Commun Appl Num Meth 1991; 7:187–94. [12] Kotronis P. Cisaillement Dynamique de Murs en Beton Arme. Modeles Simplifies 2D et 3D. PhD thesis, ENS, Cachan. 2000. [13] La Borderie Ch. Phenomenes unilateraux dans un materiau endommageable: modelisation et application a lÕanalyse de structures en beton. PhD thesis, University of Paris VI. 1991. [14] Leger P, Dussault S. Nonlinear seismic response analysis using vector superposition methods. Earthquake Eng Struct Dyn 1991;20. [15] Lemaitre J, Chaboche JL. Mechanics of solids material. Cambridge University Press; 1990. [16] Lowes LN, Moehle JP. Evaluation and retrofit of beamcolumn T-joints in older RC bridge structures. ACI Struct J 1999;96(4):519–32. [17] Mazars J. A description of micro- and macroscale damage of concrete structures. J Eng Fract Mech 1986;25(5/6):729– 37. [18] Newmark NM. A method of computation for structural dynamics. ASCE J Eng Mech Div 1959;85:67–94. [19] Paulay T, Priestley MJN. Seismic design of concrete and masonry buildings. New-York: John Wiley & Sons Inc.; 1992. 744 p. [20] Pecker A. Analytical formulae for the seismic bearing capacity of shallow strip foundations. In: Seco e Pinto PS, editor. Seismic behavior of ground and geotechnical structures. Rotterdam: Balkema; 1997. p. 261–8. [21] Przemieniecky JS. Theory of matrix structural analysis. New-York: MCGraw-Hill; 1985. [22] Ragueneau F, La Borderie Ch, Mazars J. Damage model for concrete like materials coupling cracking and friction, contribution towards structural damping: first uniaxial application. Mech Cohes Frict Mater 2000;5:607–25. [23] Rashid JYR, Dameron RA, Dunham RS. Finite element analysis of reinforced concrete in bridge seismic design
L. Davenne et al. / Computers and Structures 81 (2003) 1223–1239 practice ASCE Commitee Report In: Shing PB, Tanabe T, editors. Modeling of inelastic behavior of RC structures under seismic loads. Structural Engineering Institute; 2001. p. 217–33. ISBN 0-7844-0553-0. [24] Taylor RL. FEAP: A finite element analysis program, version 5.01 manual. University of California, Berkeley. 1996.
1239
[25] Ukritchon B, Whittle AJ, Sloan SW. Undrained limit analysis for combined loading of strip footings on clay. J Geotech Geoenv Engng March 1998:265–76. [26] Zienkiewicz OC, Talyor RL. 5th ed. The finite element method, Vols. I, II and III. Oxford: Butterworth Heinemann; 2000.