TOWARDS REALISTIC FE MODELS FOR REINFORCED CONCRETE SHEAR WALL DOMINANT BUILDINGS SUBJECTED TO LATERAL LOADING by John E KBolander Jr. Jaime'K. Wight A Report on Research Sponsored by National Science Foundation Grant No. ECE-8603624 Report No. UMCE 89-2 Department of Civil Engineering The University of Michigan Ann Arbor, Michigan January 1989

0 ->) I S"^b <$ 3. `3) I'VI'l;"~ jt

ABSTRACT TOWARDS REALISTIC FE MODELS FOR REINFORCED CONCRETE SHEAR WALL DOMINANT BUILDINGS SUBJECTED TO LATERAL LOADING A finite element based program is developed for the planar analysis of reinforced concrete shear wall dominant buildings subjected to quasi-static loadings. Twodimensional continuum elements are used to represent the shear walls which are idealized to be in a state of plane stress. Basic established material models are utilized to represent the most essential features of reinforced concrete nonlinear behavior, namely, tensile cracking of the concrete, yielding of the reinforcing steel and nonlinear behavior of the concrete in biaxial compression. Background information related to the material models employed is provided and proper operation of the program is substantiated through comparisons with test results. Loading is applied incrementally with iterations to reestablish structural equilibrium within each increment. The program traces the history of load redistribution among structural members during inelastic deformations leading up to failure conditions. The solution process also provides information useful in characterizing the i

ultimate failure mode (e.g. yielding of the flexural reinforcing steel, diagonal tensile cracking of the concrete, crushing of the concrete in compression). The program formulation could be extended to perform a more general dynamic three-dimensional loading analysis which seems to be mandatory for a more complete assessment of building performance during strong earthquake type motions. The program is used to investigate the inelastic trends leading to failure in a ten story reinforced concrete shear wall building situated in the city of Vinia del Mar along the Chilean coastline, approximately due west of Santiago. The building under consideration is typical of other Chilean buildings which, for the most part, performed exceptionally well during the March 3, 1985 earthquake experienced in that region. ii

ACKNOWLEDGEMENTS This report was submitted by John E. Bolander, Jr. in partial fulfillment of the requirement for the degree of Doctor of Philosophy (Civil Engineering) in the Horace H. Rackham School of Graduate Studies at The University of Michigan. Dr. Bolander wishes to express his sincere appreciation to Professor James K. Wight, chairman of the doctoral committee, for his guidance, encouragement, and continuous support throughout this research. The authors would like to thank Professors Subhash C. Goel, Antoine E. Naaman, and Noboru Kikuchi (members of the doctdral committee) for reviewing the report and offering helpful suggestions. The investigation was sponsored by the National Science Foundation through Grant Number ECE-8603624. The continuing support of Dr. S. C. Liu of the National Science Foundation is gratefully acknowledged The conclusions and opinions expressed in this report are solely those of the authors and do not necessarily represent the views of the sponsors. iii

TABLE OF CONTENTS ABSTRACT................................. i ACKNOWLEDGEMENTS..........................iii LIST OF FIGURES......................... vii LIST OF TABLES................................ xi CHAPTER I. INTRODUCTION........................... 1 1.1 Problem Definition and Discussion............... 1 1.2 Overview of Building Analysis................. 5 1.3 Local Modeling vs. Global Modeling.............. 6 1.4 Aspects of Material Nonlinearity................ 7 1.5 Literature Review........................ 10 II. CONSTITUTIVE RELATIONS FOR PLAIN CONCRETE. 12 2.1 Experimental Observations................... 12 2.1.1 Uniaxial stress states............. 16 2.1.2 Multiaxial stress states................ 18 2.2 Failure Theories......................... 22 2.2.1 General forms..................... 23 2.2.2 Modeling........................25 2.3 Mathematical Modeling of Constitutive Equations...... 28 2.3.1 Incremental (tangent) elasticity-based approach.. 31 2.3.2 Classical plasticity theory based approach......35 III. TENSION CRACKING OF REINFORCED CONCRETE. 47 3.1 Crack Initiation...................... 47 3.2 Crack Representation...................... 48 3.2.1 Discrete cracking models...............49 3.2.2 Smeared cracking models............... 50 iv

3.3 Crack Propagation........................ 3.3.1 Models based on strength theory......... 3.3.2 Models based on fracture mechanics concepts.... 53 55 57 IV. SOLUTION TECHNIQUES FOR NONLINEAR PROBLEMS 67 4.1 4.2 4.3 Fundamental Concepts................ Solution Techniques....................... Linear Equation Solver..................... 4.3.1 Basic methods..................... 4.3.2 Coping with limited storage availabilty....... 67 70 75 75 78 V. THE SNAC PROGRAM................. 81 5.1 Program Description.......................... 5.1.1 Element library.............. 5.1.2 Load application................ 5.1.3 Material only nonlinearity.............. 5.1.4 Convergence criteria.................. 5.1.5 Linear equation formation/solution......... 5.1.6 Constraint equations................. 5.1.7 Graphics post-processing............... 5.2 QUADC - Two-dimensional Continuum Element....... 82 82 85 87 88 89 90 90 91 92 106 116 118 118 122 127 5.2.1 5.2.2 5.3 TRUSS5.4 BEAM5.4.1 5.4.2 Non-orthogonal cracking model w/tensile softening. Plasticity based model w/abrupt tension cutoff... Two Node Truss Element............... Two Node Beam-Column Element..... Two surface plasticity model in 1-D action space.. Hinge plastic stiffness determination.... VI. CHILEAN BUILDING ANALYSIS........... 6.1 Background............................ 127 6.2 Building Description..................... 129 6.3 Modeling of the Structure................... 130 6.3.1 Upper story idealization............... 135 6.3.2 Rigid diaphram assumption.............. 138 6.3.3 Equivalent column modeling of shear walls.....139 6.3.4 Coupling beams................... 139 6.3.5 Orthogonal walls................... 140 6.3.6 Fixed base assumption................143 6.3.7 Reinforcing mesh representation...........144 6.3.8 Nonstructural elements................144 6.4 Modeling of the Loading................... 145 6.4.1 Gravity loading................... 145 v

6.4.2 6.4.3 Lateral loading................... Uniform Building Code lateral load distribution. 147.. 148 6.5 Analysis Results......................... 150 6.5.1 Global ductility measures............... 6.5.2 Nonlinear trends.................... 6.5.3 Interstory drift indices................ 6.5.4 Neutral axis migration................ 6.5.5 Failure state characterization............. 6.6 Additional Considerations.................... 6.6.1 Concrete tensile strength............... 6.6.2 Floor slab flexibility.................. 6.6.3 Lateral load distribution profile........... 152 155 163 165 167 169 169 173 174 VII. SUMMARY.............................. 179 7.1 SNAC Program Development.................. 179 7.2 Chilean Building Analysis Observations............ 180 7.3 Future Considerations...................... 182 BIBLIOGRAPHY............................... 184 vi

LIST OF FIGURES Figure 1.1 Finite element modeling of beam under two-point loading......11 2.1 Evolution of microcracks during loading.............. 14 2.2 Concrete compressive stress-strain relations.............. 17 2.3 Typical stress-strain curve for concrete in uniaxial tension...... 18 2.4 Biaxial strength envelope........................ 19 2.5 General form of concrete failure surface in 3-D principal stress space 21 2.6 Cylendrical coordinate system for describing concrete failure surface 21 2.7 Failure models........................ 26 2.8 Equivalent uniaxial strain for linear material............. 32 2.9 Loading surfaces of concrete in biaxial stress plane..........38 2.10 Graphical representation of three simple hardening rules...... 39 3.1 Nodal separation using two and four coincident nodes........ 49 3.2 Crack pattern and stress distribution................. 50 3.3 The three modes of cracking..................... 54 3.4 Three dimensional crack tip coordinate and stress systems..... 55 3.5 Dugdale and Hillerborg cracking models............... 60 3.6 General variation of stress with crack width............ 61 3.7 Size effect for plain concrete structures................ 64 vii

3.8 4.1 4.2 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 5.19 Size effects for reinforced concrete................... 65 Cracking during iterative process................... 71 Iteration techniques represented in one dimension.......... 73 Main and global driving routines.................... 83 Nonlinear solution sequence driving routine.............. 84 Decomposition of total strain...................... 93 Crack normal softening relations................... 97 Fracture specimen geometry and material properties......... 99 Fracture specimen deformation during loading............ 101 Fracture specimen load-deflection response.............. 102 Shear critical beam under two-point loading (symmetric portion). 103 Numeric cracking compared to experimental results.........105 Yield function in biaxial principal stress space............ 108 Yield function augmented with tension cut-off surface........ 110 Concrete model response in biaxial compression.......... 111 Shear wall geometry and undeformed finite element mesh......113 Shear wall load-deflection response comparison........... 114 Shear wall deformation at selected load levels............. 115 Shear wall crack pattern during testing (6top 4.9 in.)........ 116 TRUSS element bilinear material idealization............. 117 BEAM element idealization..................... 120 Hinge moment-rotation response during cyclic loading........ 121 viii

5.20 Uniaxial material properties for sectional analysis.......... 124 5.21 Moment vs. rotation at element end node............... 125 5.22 Tip load vs. tip deflection in equivalent cantilever.......... 126 6.1 Villa Real building plan views..................... 131 6.2 Villa Real cross-sections (E-W direction)............... 132 6.3 Villa Real cross-sections (N-S direction)................ 133 6.4 Building finite element model (E-W direction)............ 136 6.5 TRUSS element modeling of flexural reinforcement.......... 137 6.6 Uniform Building Code lateral load distribution........... 151 6.7 Roof lateral displacement vs. normalized base shear......... 152 6.8 Building displacement ductility ratios................. 154 6.9 Nonlinear actions during loading history............... 156 6.10 Wall nonlinear trends at t = 0.50................... 158 6.11 Wall nonlinear trends at t = 1.00................... 159 6.12 Wall nonlinear trends at t = 2.00.................. 160 6.13 Wall nonlinear trends at t = 2.79................... 161 6.14 Wall nonlinear trends at t = 2.79 (undeformed plot)......... 162 6.15 Wall nonlinear trends at t = 3.10 (westward loading)........ 164 6.16 Lateral drift profiles at several load levels............... 165 6.17 Interstory drift indices at several load levels............ 166 6.18 Migration of wall 4-4 neutral axis during loading.......... 168 6.19 Post-converged state for Iteration Cycle 2............... 170 6.20 Post-converged state for Iteration Cycle 4............... 171 ix

6.21 Post-converged state for Iteration Cycle 6............... 172 6.22 Sensitivity of building response to concrete tensile strength..... 173 6.23 Flexible diaphram idealization..................... 175 6.24 Building response with quasi-flexible diaphram........... 176 6.25 Load distribution profile response comparison.............178 x

LIST OF TABLES Table 5.1 Material parameters for shear critical beam.............. 104 5.2 Material input parameters for shear wall analysis.......... 112 6.1 Villa Real dimensions and specified material properties....... 134 6.2 Villa Real building gravity load summary............... 145 xi

CHAPTER I INTRODUCTION 1.1 Problem Definition and Discussion Within current earthquake resistant design practice for buildings, each type of building component (e.g. beams, columns, shear walls) must satisfy its own set of building code safety requirements for a specified amount of gravity and lateral loading. While the magnitude of the lateral loading is related to the overall structure type, such practice still tends to emphasize design of the individual members within the building and neglects the design of the building as a composite structure. Relatively little information is provided by the design process regarding the overall strength reserves possessed by a building due to redistribution of member loadings during inelastic deformations. In addition, potential collapse modes of the building are by no means readily identifiable. During the past two decades, considerable efforts have been directed towards the development of analysis techniques both able to account for the composite behaviour of the building components during inelastic deformations and simple enough in nature to be useful in a design setting. A popular finite element based approach is to idealize the major structural components of the building as beam column elements which have nonlinear capabilities lumped at their end nodes. These nonlinear capa 1

2 bilities are tuned to provide a degrading hysteretic behavior which agrees with data obtained from either test results or from a more refined analysis procedure. This technique, which will hereafter be called the 'equivalent column approach' is attractive in that it often reduces the computational degrees of freedom of a structure to a near minimum. The equivalent column approach has provided good approximations for the load-deflection response, story shears, overturning moments, and ductility demands for frame type systems. This same approach is often extended to cover the analysis of shear wall buildings where the shear wall is represented as an equivalent column. The distance between coupling beam end nodes and the column centerline is usually represented as a rigid zone which is effectively modeled using displacement constraints. While the equivalent column approach to modeling shear wall type structures has undoubtably contributed to our understanding of building behavior during earthquake type excitations, the approach suffers from several fundamental drawbacks. The line element approach to modeling shear walls is unable to directly account for such multi-dimensional phenomena as nonlinear strain distribution over the wall length and neutral axis migration during loading. This deficiency becomes more critical as the wall aspect ratio decreases. An alternate approach to the inplane analysis of the shear wall dominant structures would be to use two dimensional continuum finite elements for modeling the behavior of the shear walls. Such an approach circumvents many of the ambiguities and inconsistencies which arise when idealizing a two-dimensional wall as a onedimensional line element. Continuum element models applied to reinforced concrete structures have proven accurate for a variety of smaller scale problems including the analysis of beams, plates and isolated walls. One major difficulty is the increased

3 expense associated with this type of model when extending analyses to complete building systems. This may be seen as only a temporary disadvantage, however, as the benifits from current trends in supercomputing (e.g. reduced cycle times, vector and parallel processing) should soon become more widely available to researching engineers. Apart from the structural idealization, an equally important matter is the type of loading history the structure will be subjected to. Traditionally, a primary concern of designers in seismic regions has been to provide a structure with enough base shear capacity to withstand earthquake loadings prescribed according to a certain design spectra. An equivalent lateral force procedure is allowed by many building codes for simple building types [4,91]. In their current form, these equivalent lateral force methods significantly underestimate the actual forcasted earthquake forces and, in doing so, implicitly rely on a good deal of inelastic deformation, energy absorption and dissipation during the forecasted earthquake. Due to such considerations it would be useful to have some indication of building strength reserves beyond the 'code defined' load levels. One estimate of building strength could be determined by assuming a lateral load distribution profile and then proportionally increasing the profile magnitude up until the structure reaches some predefined failure state. The load profile could be chosen to approximate the fundamental mode during dynamic excitations. It must be emphasized that the lateral strength of a building is just one indication of building performance during a forecasted earthquake. The emphasis placed upon lateral strength in the building codes over the years has evolved, in part, due to a preoccupation with peak acceleration values recorded from various earthquakes. Studies have strongly indicated that high frequency spikes of acceleration are not

4 a dominant influence on the normal response of structural systems, but rather the "repetetive strong shaking characterized by strong energy content is that characteristic of the time history which leads to structural deformation and damage [44]". In light of such long duration ground motions as those witnessed in the 1985 Mexico City earthquake, it is emerging that energy absorption capacity is a fundamentally more sound basis to design upon. Equivalent lateral force methods are unable to account for the tremendous amount of energy which may be supplied to a structure during an earthquake of long duration, particularly when a structure is in resonance with the ground motion. While collapse state analyses of ductile frame structures are numerous, relatively little work has been reported with regards to predicting the failure state of structures whose primary stiffness against lateral loads is wall systems. Building trends in the United States favor the ductile frame system over the shear wall system, however the shear wall system is popular abroad in such countries as Japan and Chile. As discussed in Chapter VI of this study, the reinforced concrete Chilean buildings performed quite well during the March 3, 1985 earthquake experienced in that region. Moreover, design practice for reinforced concrete shear wall type structures in Chile neglects the structural detailing deemed necessary to insure ductile behavior under extreme loading conditions. Looking over the analyses stemming from earthquake reconasance reports, there seems to be a preponderance of studies investigating heavily damaged buildings and the causes. Less attention has been devoted to analyses of buildings which performed relatively well. The building analysis performed in this study is an attempt to reduce this imbalance.

5 1.2 Overview of Building Analysis Within the present study a finite element based numerical technique has been developed for predicting the progressive inelastic actions leading up to the failure state of shear wall dominant buildings subjected to inplane quasistatic lateral loading. For lack of a better acronym, the program has been named SNAC (Simple Nonlinear Analysis of reinforced Concrete structures). Due account is given to the inelastic material behavior of reinforced concrete through application of basic material models taken from the wealth of literature on the subject. A failure state analysis will be performed for the Villa Real building, a typical building designed according to Chilean building code practice. This building is located in the coastal city of Vinia del Mar and is one of several Chilean buildings being studied in a United States National Science Foundation funded joint effort between the University of California at Berkeley, the University of Illinois (Urbana), and The University of Michigan. For the most part, these buildings performed well during the March 3, 1985 earthquake (magnitude M, = 7.8) experienced in that region [90,80,94]. The Chilean buildings under consideration in the joint project may all be considered as reinforced concrete shear wall dominant systems and vary in height from 10 to 23 stories. The Villa Real building structure will be subjected to lateral load distributions approximating the fundamental mode response. This load distribution will be monotonically increased up to a level at which overall structural failure occurs. Failure will be assumed to occur when either the structural displacements become unacceptably large or the structure begins to lose its load carrying capacity. The solution process developed within this study will trace the history of load distribution among

6 structural members during inelastic deformation leading up to the failure state. For a given building, the magnitude of the lateral load distribution at the failure state will be compared with the magnitude of the Uniform Building Code specified lateral load distribution, thus providing some indication of the strength reserves possessed by the overall structural system for the type of loading considered. In addition, the solution process also provides information useful in characterizing the corresponding ultimate failure state of the structure (e. g. crushing of the concrete in compression, tensile cracking of the concrete, yielding of the main reinforcing steel). Results will be compared to those obtained from an equivalent column modeling of the shear walls. Further details regarding the Villa Real building description and analysis are presented in Chapter VI of this study. 1.3 Local Modeling vs. Global Modeling Roughly speaking, most reinforced concrete finite element analyses may be classified into two types, namely: 1) those whose primary concerns are accuracy in modeling the localized behavior of the reinforced concrete constituents within the structure, and 2) those which are primarily concerned with accurately representing global behavior of the reinforced concrete structure. The former analysis type, based on what may be called a "micro-model" approach [43], through detailed modeling of the reinforced concrete constituents provides a means of realistically describing the localized physical phenomena such as crack propagation, crack width and spacing, transfer of stress across cracks, and bond stresses. A micro-model based study should also provide an accurate solution to the global structural behavior, albeit at the cost of considerable computational effort and program formulation complexity. A micro-model analysis is generally most useful for

7 research purposes where the fundamental nature of reinforced concrete under loading is of primary interest. Meaningful results have been attained for investigating such problems as crack propagation, shear transfer mechanisms, and reinforcing bar slip or pull-out. The latter type of analysis, based on what may be termed a "macro-model" approach [43], attempts to provide an expedient solution for the overall structural behavior often foresaking information regarding localized details in the process. The effects of the reinforcing steel or cracking in the concrete are usually assumed to be distributed or "smeared" out in a continuous fashion over the finite element(s) in question. Such models have been successful in describing the structural response for a wide variety of reinforced concrete problems, such as those involving concrete beams, shear walls, containment structures, and reactor vessels. The proposed study will approach the building analysis from a macro-model viewpoint, gaining simplicity in formulation and economy in program execution while foresaking detailed information regarding localized behaviors. The consequences of employing a macro-model based analysis, and in particular a smeared cracking approach, will be explored in more detail later in this report. 1.4 Aspects of Material Nonlinearity It seems most appropriate to present an introductory discussion on the nonlinear material modeling at this point, as it was undoubtably one of the most challenging and consequential aspects of this study. Indeed, the material modeling in an inelastic study is often recognized as the weakest link in the analysis chain. This is especially true for reinforced concrete where a number of complex and highly nonlinear material behaviors may exist. Considerable research has been performed over the past two

8 decades in an attempt to find material models which can adequately describe the nonlinear behavior of reinforced concrete [88,92]. However, the success of a material description has often been application dependent and as of this date there is still little agreement among researchers as to what material models should be adopted for general analysis purposes. When attempting a failure load analysis of a reinforced concrete structure one must identify and account for the essential nonlinearities arising in the reinforced concrete for the type of loading considered. Causes of composite nonlinear behavior include: 1. nonlinear behavior of concrete under multiaxial stress states, 2. tensile cracking of the concrete, 3. nonlinear behavior of the reinforcing steel, 4. bond slippage at the reinforcing steel-concrete interface, 5. doweling action of the reinforcing steel across cracks, 6. stress transfer across cracks due to aggregate interlock. One should recognize that some of the above causes of nonlinear behavior may play an important role in the overall load-deflection behavior of the structure, while the others might be of lesser importance. Given the amount of computational effort required for analyzing even a simple reinforced concrete structure, a collapse load analysis of an entire building system requires gross simplyfing assumptions in order to remain computationally feasible. With this in mind, the modeling of the more dominant nonlinear effects will be pursued vigorously, while other less significant effects might only be roughly accounted for, or neglected altogether. Considering

9 the type of loading proposed, it is the author's viewpoint that the tensile cracking of the concrete will dominate the nonlinear response of the structure during the early stages of the loading history with the nonlinear behavior of the concrete under high compressive stress becoming more significant as the ultimate failure state is approached. Careful consideration will be given to modeling these two aspects of nonlinear behaviour in Chapters II and III of this report. Further insights into these and the other factors listed above will arise when discussing the element material routines in Chapter V. It should be remembered that many of the analysis input parameters, the ground motion, the distribution of reactive live load masses, and to varying degrees the material properties, are probalistic rather than deterministic in nature. Many assumptions regarding the load application and structural idealization also involve gross simplifications. A high degree of refinement in the material models and the solution procedure relative to these other factors may not be warranted and in fact may lead to a false sense of security in the analysis results. Finally, it should be noted that the ultimate goal of current reinforced concrete modeling research is to develop analytical models capable of predicting the structural response to general loading histories. Within this study, frequent reference will be made to research activities related to problems involving three dimensional material behavior, post-failure response, and reversed loading. While such considerations may not appear to be within the scope of this study, however, a basic understanding of such topics is nonetheless essential for continued work in this field of study.

10 1.5 Literature Review The earliest published work concerning the finite element idealization of a reinforced concrete structure is that of Ngo and Scordelis in 1967 [67]. Their study investigated the behavior of simply supported concrete beams under primarily flexural loading. The actual three-dimensional reinforced concrete beam was modeled using a two-dimensional plane stress analytical model. The finite element idealization of the beam is shown in Figure 1.1. The steel and concrete were each modeled separately using simple constant strain triangular -finite elements. The steel-concrete bond characteristics were represented by special 'linkage' elements which were spaced periodically along the steel-concrete interface. Cracking was constrained to follow along predetermined inter-elemental boundaries. Using the above modeling assumptions in conjunction with finite element solution procedures, the above study was able to retrieve the displacement components of each nodal point, the concrete or steel stresses in each element, and the forces in the special bond linkage elements. Since that time, the amount of material published on the analysis of reinforced concrete by the finite element method is exhaustive and hence no attempt will be made here to make a complete summary. Instead, the interested reader is referred to two 'state-of-the-art' works published by The American Society of Civil Engineers [88,92]. The first of these works, published in 1982, is a collection of reports submitted by an ASCE Task Committee formed in 1978. This work provides a broad coverage of the material essential to reinforced concrete analysis by finite element methods. The more recent publication is a collection of 40 papers presented by researchers at the Joint U.S.-Japan Seminar on Finite Element Analysis of Reinforced

11 -Typical Triangular / Concrete Element + y/ f Steel Reinforcement..... 222 ZZ. 22 *: F K.'::x::':::...... -. -. -:.:,i ":'''::''::''v-':-~~============ = = =:: ' ':":::-:::::-::::::::'[.Discrete Cracksi LTypical Triangular Steel Element Figure 1.1: Finite element modeling of beam under two-point loading [67] Concrete Structures held in Tokyo, Japan (May 1985). Many of these papers are application orientated and more focused in subject matter relative to the previous ASCE publication. That the field of reinforced concrete analysis continues to attract a seemingly disproportionate amount of attention may be due, in part, to the fact that many important questions remain unresolved.

CHAPTER II CONSTITUTIVE RELATIONS FOR PLAIN CONCRETE 2.1 Experimental Observations Before proceeding to characterize concrete stress-strain relations it is instructive to examine the mechanisms of internal damage within concrete material during loading. Concrete stress-strain behavior is directly related to these mechanisms of internal damage (i.e. microcracking) and the propagation of these microcracks during the loading process can be characterized using fracture mechanics related concepts. Even before any external loading is applied numerous microcracks exist within the concrete material, mostly along the aggregate-matrix interfaces. The number, location, and extent of these preexisting microcracks depend mainly on the following factors [38]: * type of cement * mineralogical nature of aggregate * geometry of aggregate * water/cement ratio * curing conditions 12

13 For loading up to a proportional limit the lengths of these microcracks remain virtually unchanged (i.e. the material behaves in nearly a linear elastic manner). However, for loading beyond the proportional limit the microcracks will begin to propagate, first along paths of least energy demand and later continuing into more resistant media as the loading is increased. Often in concrete the aggregate-matrix bond strength is less than that of the matrix itself, which in turn is less than the strength of the aggregate. In such cases initial cracking would occur along the aggregate-matrix interface and only later proceed into the surrounding matrix, the material therefore showing some increase in resistance to fracture. Note that this increase in fracture resistance is a necessary criterion for slow crack growth to occur since there must be a continuous balance between released and consumed energy during slow stable crack propagation [52,26]. One can see that it is precisely the multiphase nature of concrete which enhances its resistance to fast crack propagation relative to more brittle materials such as hardened cement paste or mortar. Upon further loading microcracks propagating into the matrix eventually bridge (Fig. 2.1) forming larger continous cracks which eventually propagate in an unstable manner leading to failure of the material. That is, when cracks reach their critical length the available internal energy of the system becomes larger than the required crack release energy and the crack then propagates with increasing rate. It has been determined that the evolution of these microcracks during the loading process depends mainly on the following factors [38,49,68]. * aggregate/matrix stiffness ratio * strength of aggregate-matrix bond * percentage of voids in matrix

14 INCREASE OF THE STRESS FIELD O IIA.,IN / ORIGINAL v INCREASE CRACK DEBONOING \ IN OEBONOING BRANCHING If — ^-CRACK. ARREST 4 ^BRIDGING CRACKS IN MATRIX. LINKING. VO. Taken from DiTonmaso (38] Figure 2.1: Evolution of microcracks during loading With regards to experimental work, the majority of concrete material testing has been performed for the uniaxial loading case. Fewer results are available for biaxial loading situations and results for the triaxial loading case are even more scarce. The uniaxial results, and to a larger extent the multiaxial results, show significant scatter due to variance in both materials and testing apparatus. The degree of "confinement" unintentionally provided by the testing apparatus greatly influences the strength and failure mode of the specimen. In general, the more confinement the stronger and more ductile the specimen will behave. It should also be noted that testing programs have concentrated almost exclusively on proportional loading histories. Nonproportional loading (i.e. loading in which the directions of principal stress rotate) may markedly affect the concrete response [10].

15 Both the complexity of concrete material characterization and response on the microstructure level and the relative scarcity of test results for general loading histories should be kept in mind when developing mathematical models describing the constitutive behavior of concrete. To date there is little agreement as to a model which effectively captures the essential aspects for general loading. Historically, models have been phenomenological in nature [47,68]. Furthermore, nearly all existing models treat concrete as a homogeneous material in spite of the fact that the composite nature of concrete plays a prominent role in shaping its material behavior. Such models don't really represent failure theories, but are rather just sophisticated mathematical descriptions of macroscopically observed test results. It follows also that a large data base of test information is required in order to adequately simulate concrete or rock response over a wide variety of load or deformation histories. More recently, researchers have turned their attention towards employing a combination of micromechanical considerations and so-called "homogenization" (or averaging) techniques. The composite nature of concrete may be accounted for by treating the matrix and aggregate as separate materials and then predicting the overall response using the theory of interacting continua or the theory of mixtures [68]. These approaches build on an increased physical basis and therefore provide hope in reducing dependence on test results. However, the random and chaotic nature of the concrete microstructure, which depends upon many factors, essentially precludes the development of any practical material response models based on purely physical theory. It seems some degree of phenomenological modeling will always be necessary. Throughout the remainder of this report repeated reference will be made to the failure surface (or curve) of concrete under general stress states. Failure is an arbitrary definition and can based on yielding of the material, initiation of cracking, load

16 carrying capacity, extent of deformation, and so forth. Each failure definition has its own merits depending upon the application. Henceforth, the failure surface of the concrete will be regarded as an ultimate strength envelope based on either limiting stress or limiting strain criteria. 2.1.1 Uniaxial stress states The uniaxial stress strain relations are often taken as a starting point when describing the material behavior of concrete. This is probably due both to the conceptual simplicity of uniaxial relations and to the fact that these relations are well established through the vast number of test results which have been reported. Typical stress-strain and stress-volumetric strain relationships for concrete under uniaxial compressive loading are shown in Figure 2.2. The concrete behaves almost as a linear elastic material for loading up to proportional limit of about 30 percent of its ultimate load carrying capacity, fi. As the loading is increased the stressstrain curve shows first a gradual increase in curvature related to microcracking along the aggregate-matrix interface (debonding). During subsequent loading, a more dramatic increase in curvature is seen as the microcracks propagate into the mortar and eventually bridge in the form of continuous mortar cracks. The load at which continuous mortar cracks begin to propagate in an unstable manner has been termed "critical" stress and coincides with the point of minimum volumetric strain. Beyond this point, the volumetric strain increases (dilitancy) and Poisson's ratio apparently increases. One must think in terms of an apparent Poisson ratio since the cracked material may no longer be idealized as a continuum. Beyond the peak load the stress-strain curve exhibits a descending, or "strain softening", branch until all strength is lost at some value of ultimate strain, E,.

17 Critical stress eVeV = el1 + e2 + e3 Taken from Chen [32] Figure 2.2: Concrete compressive stress-strain relations Plain concrete is relatively weak in tension, measured uniaxial tensile strengths being on the order of E/1000, where E represents Young's modulus. A typical stressstrain curve for a concrete element in uniaxial tension is shown in Figure 2.3 (note the similarity to the compressive stress-strain curve of Figure 2.2). The concrete behaves almost linear elastically up to roughly 60 percent of its maximum tensile load carrying capacity, ft. With increasing loading the stress-strain curve shows a dramatic increase in curvature as the tensile stress limit, f, is approached. The magnitude of stress at which tensile failure initiates and the precise shape of the post-peak region are the subject of much debate due, in part, to the variabilty of concrete material constituents and the experimental difficulties encountered. Theoretically, an ideal continuum material should show no strain softening branch as it implies instability. The heterogeneous nature of concrete material allows a strain softening branch through unstable strain localizations in the material [9]. It is interesting to note that experiments have been able to produce a strain softening branch

18 a Af Micro cracking/ Fracture Linear Elastic Figure 2.3: Typical stress-strain curve for concrete in uniaxial tension only when the size of the inhomogeneities (i.e. the size of the aggregates) is not too small relative to the overall specimen dimensions. 2.1.2 Multiaxial stress states With regards to biaxial stress states several studies deserve noteworthy attention [56,58,89]. In particular, the test results of Kupfer, et al. [56] have generally been regarded as an accurate representation of concrete biaxial stress-strain relations for proportional loading. In their study, steel brush platens were used to reduce the confining effects of the loading apparatus. These results are plotted in the form of a biaxial strength envelope where the rays directed outward from the origin indicate proportional loading paths (Figure 2.4). When comparing the biaxial response with that of the uniaxial case one sees that there is, in effect, a stress induced anisotropy. Studies have also indicated that the biaxial strength envelope is virtually independent of loading path [88].

19 Figure 2.4: Biaxial strength envelope [56] The biaxial concrete material response may be divided into three separate regions of interest: 1. Compression-compression: biaxial compression in effect impedes microcrack propagation leading to increases in ultimate strength. For a principal stress ratio of about al/a2 = 0.50 ultimate strength increases to roughly 1.25f,. For al/a2 = 1.00 an ultimate strength of roughly 1.16ff is obtained. 2. Compression-tension: the tension in one principal stress direction in effect fosters microcrack propagation leading to failure at drastically reduced strength levels. 3. Tension-tension: little stress induced anisotropy is witnessed as the material strengths in each principal, direction are roughly equivalent to their uniaxial

20 values. Test results for concrete under triaxial loading conditions indicate a failure surface of the form shown in Figure 2.5 in principal stress space. For reasons to be described later, it is convenient to describe to failure surface not in principal stress space, but rather in a cyllendrical coordinate system (see Fig. 2.6) whose generator is the hydrostatic axis (i. e. al = a2 = a3 where ai are the principal stresses). Planes containing the hydrostatic axis are called meridian planes. Planes normal to the hydrostatic axis are called deviatoric planes. Correspondingly, the state of stress at a point, represented by the Cauchy stress tensor, oij, can be decomposed into a hydrostatic (or mean normal) stress part, a,, and a deviatoric stress part, sij, as follows [60]: aij = s ij+ aij (2.1) where 6ij is the Kronecker delta. The deviatoric stress tensor by itself represents a state of pure shear and may be thought of as the " deviation" from a hydrostatic stress state. In an analogous manner the state of strain at a point, ij, may also be decomposed into hydrostatic, e,, and deviatoric, eij, strain components. From Figure 2.5 one sees that the failure surface is roughly in the shape of a paraboloid open in the direction of hydrostatic compression. Cross-sections of the failure surface in the deviatoric plane transition from a nearly triangular shape close to the tensile limit to a more circular appearance as one proceeds in the direction of increasing hydrostatic compression. It follows from an isotropic material assumption that the labeling of the principal stress directions is arbitrary and one therefore sees three-fold symmetry of the failure surface in the deviatoric planes. Six-fold symmetry is not realized in concrete materials since the failure strength in simple tension is generally not equal to the failure strength in simple compression. A practical result of

21 Elastic limit \ surface O // /N/ ^ -02 Taken from Chen [32] Figure 2.5: General form of concrete failure surface in 3-D principal stress space a1 l p F 6 a3 Figure 2.6: Cylendrical coordinate system for describing concrete failure surface

22 the symmetry property is that testing need only be performed within a 0~ < 8 < 60~ sector of the failure envelope, where 0 is defined as the angle of similarity. Most experiments have been performed for stress states along the bounding meridians of this 60 degree sector. The intersection of the failure surface with the meridian plane 0 = 00 corresponds to a hydrostatic stress state with a tensile stress superimposed in one direction and therefore is called the tensile meridian of the failure surface. Special cases found on this meridian include the uniaxial tensile strength, ft, and the equal biaxial compressive strength, fH. The intersection of the failure surface with the meridian plane 0 = 60~ corresponds to a hydrostatic stress state with a compressive stress superimposed in one direction and therefore is called the compressive meridian of the failure surface. Special cases found on this meridian include the uniaxial compressive strength, f, and the equal biaxial tensile strength, fb. 2.2 Failure Theories A mathematical description of the concrete failure surface is often necessary when performing analyses of reinforced concrete structures. For example, a failure surface may be used in conjuction with an elasticity based constitutive model for an ultimate load analysis of a structure. Alternatively, failure surfaces may be used to define the initial yield and subsequent loading surfaces employed in plasticity based constitutive formulations. These concepts will be discussed in more detail in Section 2.3.

23 2.2.1 General forms It is of primary importance that the mathematical description of the failure surface be invariant under changes in frame of reference (the principal of material frame indifference is one of the fundamental postulates of a purely mechanical theory [60]). For example, failure criteria based upon the state of stress must be invariant functions of the state of stress, that is, independent of the choice of frame of reference by which stress is defined. The general functional form of the failure criteria based upon the state of stress may therefore take suitable invariant stress measures as arguements (i. e. functions of invariant quantities are also invariant). The principal stresses are invariant functions of the stress state and the form of the failure surface may be expressed as: f(Ol,7a2,oT3)=0 (2.2) However, a failure function in terms of principal stresses is often difficult to formulate and also provides little physical or geometric insight into the failure conditions [32]. Instead, a failure function in terms of certain other invariants of the stress state appears advantageous. The following function is particulary suited to concrete and rock-like materials [88]: f(Il, J2, J3) =0 (2.3) where I, is the first invariant of the Cauchy stress tensor, aij, and J2, J3 are the second and third invariants of the deviatoric stress tensor, sij, respectively. Details regarding these and other stress (or strain) tensor invariants may be found in [60]. Consider a plane which passes through a point in a stressed body and makes equal angles with each of the principal stress directions (called octahedral plane). It can be shown that the normal stress on the octahedral plane, aoct, and the shear on

24 the octahedral plane, Tact, are related to the first two stress invariants as follows [32]: oct = I, = a, (2.4) Toct= J2 (2.5) Furthermore, it can be shown that the direction of the octahedral shear stress is related to the J3 invariant through the angle of similarity, 9: cos 30 = (2.6) Toct From the above discussion it is seen that the three invariants, I,, J2, and J3, possess some physical meaning. It also follows that the failure surface may be expressed by the function: f(Coct, oct, cos 30) = 0 (2.7) To obtain a geometric interpretation of these invariants consider the state of stress as being represented as a vector OP as shown in Figure 2.6. Resolving OP into a component ON of length 5 along the hydrostatic axis, and a component NP of length r within the deviatoric plane, it can be shown [32]: = I = = am = Vaot (2.8) r = VJ2 = VToct (2.9) Again, it can be shown that the J3 invariant is directly related to the angle of similarity, 0, which expresses the direction of NP in the deviatoric plane (see Eq. 2.6). It follows that the failure surface may also be expressed by the function: f(, r,) = 0 (2.10) While the above discussion on three dimensional failue surfaces is not essential to this study, it has been included here for the sake of completeness. It should be

25 noted that the above representations of the failure surfaces are essentially the same and one is chosen over the other merely as a matter of convenience for the problem at hand. The description of the failure surface is an important aspect in many engineering analyses and a basic understanding of the aforementioned models will also prove useful in modelling failure surfaces of other engineering materials, many of which take similar forms. 2.2.2 Modeling Earlier analytical studies of reinforced concrete expressed the failure surface in simple forms based on just one or two material parameters. However, these simple models invariably fail to describe one or more aspects considered essential to the concrete failure surface. For example, the von Mises formulation represents the failure surface as a circular cyllendar whose generator is the hydrostatic axis. It shows no hydrostatic pressure dependence and failure is therefore a function only of shear. In particular, the octahedral shear stress is the key variable considered to cause yield. The von Mises criterion may be expressed most conveniently in terms of the J2 invariant: f(J2) = J2- k2 = (2.11) where k represents the yield stress in pure shear. By itself the von Mises is a one parameter model. This type of formulation has traditionally been applied to ductile metals which show this type of yield based on shearing stress. The von Mises formulation has also been applied to reinforced concrete in tandem with a tensile stress cutoff (we now have a two parameter model) to account for the limited tensile capacity of concrete [86]. The tension cutoff reasonably models the behavior of concrete in the

26 02 Q3 (a) Von Mises 'r0cst 8,e ' a'i roct 9.6 ct 4*3 (b) Mohr-Coulomb roct (c) 3rucker - Prger (c) Orucker-Prager roct Ci Thc- octt C.2 C.3 (e) Chen- Chen ict 60oct 02 03 (f) Hsieh-Ting-Chen Toct 860 so iaoct - War (g) Wi Iliam-Warnke () se'P ) Oto (d) Bnelw - Piswr (h) Ottow Taken from Hegemier, et al. [47] Figure 2.7: Failure models so-called "brittle region" of the response and the failure surface does approach the shape of the von Mises cyllendar for stress states having a high hydrostatic stress component (the so-called "ductile region"). However, failure surface dependence on both the hydrostatic pressure, (, and the angle of similarity, 0, in the intermediate hydrostatic stress range is not accounted for. Cross-sections of the failure surface in both the meridian planes and the deviatoric plane are compared for several popular models in Figure 2.7. The Mohr-Coulomb model used extensively in soil mechanics has also been applied to concrete materials. This is a two parameter model based on the cohesion and angle of internal friction of the material. A tensile cutoff may also be utilized

27 making it a three parameter model. The Mohr-Coulomb model has the principal advantages of being both simple and able to provide a linear failure surface dependence on hydrostatic stress which is a good approximation for concrete for the intermediate range of hydrostatic stress levels. However, this model rapidly deviates from experimental results as hydrostatic pressure increases. The model also has numerous other shortcomings as outlined in [32]. Another two parameter model is the Drucker-Prager model which is basically an attempt to overcome the numeric difficulties associated with the hexagonal corners of the Mohr-Coulomb model. It becomes an extension of the von Mises formulation by simply adding a term to account for hydrostatic pressure dependence of the failure surface: f(1, J2) = SI1 + -k=O (2.12) where P3 is a constant. While the Drucker-Prager model produces a smooth failure surface (i. e. no numeric difficulties with corners), it too is based on a linear dependence between I, and /7 and therefore loses validity under high hydrostatic pressure states. In addition, the model shows no dependence upon the angle of similarity, 0. The Brestler-Pister model and the Chen-Chen model are both three parameter representations which have fitted a parabolic dependence of rToc on aoct (or VJ' on II) to match test results. Again there is no dependence on the angle of similarity. The remaining models, which range from four to five input parameters, incorporate both a parabolic roct-Coct relation plus an angle of similarity dependence. The William-Warnke model, for example, compares quite well with triaxial test results [88]. It should also be noted that~ failure models have been developed specifically for

28 the biaxial loading case. In [55] analytical expressions in terms of principal stresses have been fitted to Kupfer's [56] biaxial test results. An expression is provided for each region of the biaxial loading response: (- + f2)2 + + 3.65-2 = 0 compression-compression (2.13) fc fJ fc Jc 02 a0 = 1 + 0.8-1 compression-tension (2.14) ft fc r2 = ft tension-tension (2.15) All the above failure models are described in terms of stress variables. Test results have shown and it also seems reasonable that in certain situations failure may depend on the state of strain. It might therefore be beneficial to utilize "dual" failure criteria (i. e. two independent failure surfaces, one for the state of stress and one for the state of strain) [32]. An alternative is to base the failure surface function simultaneously on both stress and strain invariants [16]. 2.3 Mathematical Modeling of Constitutive Equations Numerous constitutive models exist for describing the nonlinear response of concrete or rock materials to short term loading. As of 1982 these models could, for the most part, be classified into one of four general catagories, namely [32,88]: 1. Elasticity theory and its generalizations 2. Perfect and work hardening plasticity theory 3. Plastic fracturing theory 4. Endochronic theory of plasticity Intense research in the area of concrete and rock constitutive modeling has introduced several new models in recent years which may not fit into one of the aforemen

29 tioned catagories. A representative cross-section of noteworthy three dimensional models is listed below [47]: * Nonlinear elasticity (Cedolin et al., 1977), * Classical plasticity (Chen and Chen, 1975), * Hypoelasticity (Coon and Evans, 1972), * Endochronic plasticity (Bazant and Bhat, 1976; Bazant and Shieh, 1978), * Bounding surface plasticity (Fardis et al., 1983; Yang et al., 1983), * Continuous damage theory (Krajcinovic and Selvaraj, 1983; Suaris and Shah, 1983), * Plastic fracturing theory (Bazant and Kim, 1979), * Hysteretic fracturing endochronic theory (Bazant and Shieh, 1980), * New endochronic theory with singular kernels (Valanis and Read, 1985). A detailed discussion on the relative merits of each of the above constitutive modeling schemes is well beyond the scope of this report. Most of the models accurately predict the concrete response for relatively simple loading histories which neither cause the principal stress directions to vary significantly nor subject the concrete to extreme stress states near or beyond ultimate strength [88]. More recent models, such as those based on endochronic theory, have been able to provide a more comprehensive modeling of concrete response for monotonic as well as cyclic loading situations. Endochronic theory based models have had success in predicting the correct volume changes of concrete for general stress states, an important consideration when investigating the effects of confining reinforcement. Furthermore, the strain softening branches and failure conditions are incorporated into the constitutive relations and,thus, these models are free of cumbersome yield criteria and hardening rules necessary in a plasticity based approach. However, current formulations based on endochronic theory are both relatively complex and computationally expensive. In addition, these models show a dependency on a large number of material input

30 parameters many of which no experimental data base exists for the problems in question. As previously mentioned, damage mechanics based models are attempting to relieve some of this dependence on test data by basing the concrete response on actions occuring at the microstructural level. While the damage mechanics approach shows great promise, the theory is still in its infancy and it too appears to demand a considerable degree of complexity in order to be of general nature. Taking such considerations to mind, the author has decided to base the program computations outlined later in this report on one of the simpler formulations taken from the extensive wealth of literature on this subject. Within a total (or secant) elasticity-based formulation the current state of stress, aij, is uniquely determined as a function of the state of strain, eij, or visa versa. The constitutive relations may be expressed in the form: aij = Fij(ekl) (2.16) where Fij is a tensorial material response function. The unique relation between stress and strain implies both irreversability and path independence which is certainly not a valid assumption for concrete behavior in general. To address this shortcoming, this type of formulation has been extended by introducing the concept of a loading function, f(aij), in terms of the state of stress. The quantity df = af daij may then be used as an indication of whether loading or unloading occurs. For loading (df > 0) the material responds according to Eq. 2.16. For unloading (df < 0) the material responds according to some other relation, often by an isotropic elastic model. Problems arise when either a neutral loading condition (df = 0) occurs or when unloading and loading occurs when the material response function, Fij, is not symmetric as pointed out in [88] and [32], respectively. For the type of loading considered in this study (monotonic proportional) this

31 secant type of formulation is attractive as no unloading would occur and the whole question of a loading function and its associated difficulties becomes academic. However, such a method is not representative of concrete behavior in general and is not in the spirit of extending to more general loading histories. 2.3.1 Incremental (tangent) elasticity-based approach Incremental elasticity-based formulations are useful in describing material response when the state of stress depends not only on the current state of strain, as in Eq. 2.16, but also on the stress (or strain) path to reach that state. Due to the path dependent nature of concrete materials this is a much more realistic than the total formulation when unloading or nonproportional loading occurs. Incremental stress daij and strain dEij tensors are assumed to be linearly related through material response moduli which depend on either the state of stress, the state of strain, or both '88,32]: daij = Dijkl(mn)dekl (2.17) dij = Dijkl(mn)dEkl (2.18) where Dijkl is a tangential material stiffness tensor of fourth order. The material response described by Eqs. 2.17 and 2.18 is reversible only in an infinitesimal (or incremental) sense which has led-some authors to describe this type of constitutive equation as "hypo-elastic" [60] or "minimum elastic" [32]. The above model has been well able to predict various test results but some researchers argue that the model loses validity under higher stress states as if fails to account for the coupling between hydrostatic and deviatoric components witnessed in testing. Researchers have also questioned the use of isotropic constitutive relations as concrete shows pronounced stress induced anisotropy under multiaxial loading.

32 -0I, Bioxial Loading cr2 ~ I-v cZ Unioxiol Stress-Strain Curve I a2 -I { r2 al ' a2 I T I i 7<u-T -( Taken from Darwin, et al. [36] Figure 2.8: Equivalent uniaxial strain for linear material The problem of stress induced anisotropy has been addressed for the biaxial loading case by Darwin and Pecknold [36,37] who proposed the concept of equivalent uniaxial strain. As an introduction to the equivalent uniaxial strain concept, consider the linear material response to biaxial loading (a = alla/2 1) shown in Figure 2.8. For a given stress level, the true strain in the major compressive direction may be found from Eq. 2.19 where the term containing Poisson's ratio effectively stiffens the response for increasing principal stress ratio a. = - a (2.19) e For the same stress level the equivalent uniaxial strain may be found using the same relation but with the Poisson effect removed. That is, for a biaxially (or more generally triaxially) loaded elastic material, the equivalent uniaxial strain in the ith

33 principal direction, eiu, is: = sE (2.20) Eo The concept is extended to cover nonlinear material response by defining the equivalent uniaxial strain in principal direction i as the summation of equivalent uniaxial strain increments in direction i over each step in the loading path. Ei = E a.: = E a (2.21) where Aoi is the incremental change in principal stress, Ei is a varying tangent stiffness for the corresponding load increment, and the summations are taken over all load increments. It must be emphasized that the equivalent uniaxial strains are not real strains, but rather fictitious strains which serve only as a basis from which to alter the material properties to account for the effects biaxial stresses on the internal damage state of the concrete. The material properties during any load increment are based on equivalent uniaxial stress-strain curves which take the general form [88]: ai =,i(ei,, a, loading history) (2.22) In [37] the equivalent uniaxial stress strain curves for biaxial compressive loading are based on an equation of Saenz [83] which is generally accepted for describing the ascending branch of one-dimensional stress-strain curves: a= (2.23) 1 = + [(Eo/E,) - 2](e/)+ (/)2 (2.23) where ec is the strain corresponding to the maximum compressive stress, fc and Es is the secant modulus E, = fc'c. Writing Eq. 2.23 in terms of equivalent uniaxial strain to approximate the biaxial

34 case, one gets the equivalent uniaxial stress-strain relation for principal direction i: = Eociu (2.24) i = 1 + [(Eo/Es) - 2](eiu/,c) + (iu/)2 (2.24) where aic is the maximum compressive stress, eic is the equivalent uniaxial strain corresponding to maximum compressive stress aic, and Es is the secant modulus (= aic/eic) The equivalent uniaxial stress-strain relations purposefully do not depend on Poisson effect. The effect of principal stress ratio is contained primarily in ai, and ec. For a given principal stress ratio a, the value for the current tangent modulus in principal direction i, Ei, may be found from the slope of the a, —i curve (Eq. 2.24) at the current value of ei, which has accumulated during the loading history according to Eq. 2.21. dai E,_ d= (i = 1,2) (2.25) dciu These tangent moduli values are substituted into incrementally orthotropic constitutive relations. fdal E1 i. /jE E 0 de da2 E = 0 de2 (2.26) dr12 sym '(Ei + E2 - 2pv/T ) d1 in which p = /vivi is an "equivalent" Poisson ratio. Poisson's ratio is therefore used in describing the overall material response in Eq. 2.26, but does not enter into the determination of the material moduli in Eqs. 2.25. The equivalent uniaxial strains are accumulated in each principal stress direction which, in general, change during the loading process. It follows that Ce, does not provide a strain history in a fixed direction, but rather in a continually changing direction corresponding to the principal stress ri. This is a desirable facility in that it

35 somewhat accounts for rotation of the principal stress directions during nonproportional loading. However, the concept has been subjected severe criticism. In [10,88] it is pointed out that damage accumulates within the concrete microstructure in the form of defects (microcracks) during the loading process. The defects will, in almost all cases, have a directional nature (e.g. microcracks forming perpendicular to the direction of principal tensile strain). In the equivalent uniaxial strain model, rotation of orthotropic axes during nonproportional loading implies rotating these defects against the material. This in turn implies the directional nature of these microstructural defects is soley a function of the current state of stress, a sort of path independence which is not true of concrete behavior in general. 2.3.2 Classical plasticity theory based approach Referring back to the experimental observations reported in Section 2.1, concrete under compressive loading was seen to behave in a nearly linear elastic manner up to a point at about 30 to 50 percent of the ultimate stress, fL, called either the proportional limit, discontinuous stress, or elastic limit. Further loading activates the mechanisms of internal damage and the material begins to deform in an irreversable manner. It appears, from a phenomenological viewpoint, that the material deformation is similar to the permanent plastic deformation in metals. Such observations lead researchers to employ classical plasticity theory based formulations to concrete. Two basic approaches to developing the constitutive relations for materials showing plastic behavior exist within the theory of plasticity [32,51], namely: 1) Total strain or deformation theory of plasticity, and 2) Incremental or flow theory of plasticity. The deformation theory of plasticity allows for the determination of finite relations

36 connecting stress and strain for material under simple loading conditions (i.e. when the components of the stress deviator, sij, increase in proportion to one parameter). Under some instances the use of this formulation can be extended to more general loading histories. For monotonically increasing loads it is immaterial whether work is dissipated in the form of heat or stored in the form of elastic potential energy. Along these lines, it can be shown that the deformation theory stress-strain relations are equations of a nonlinear elastic material. Flow theory based models are incremental (or differential) in nature as they relate the increments of plastic strain, deij, to both the increments of stress, daij, and the state of stress itself, aij. Unlike deformation theory models, flow theory plasticity models are applicable to the more general case where the stress at any state depends on the strain path. Flow theory plasticity models are therefore more appropriate for situations where unloading or nonproportional loading occur. Simple constitutive models based on the flow theory of plasticity are considered within this study. Concrete will be idealized as an elastic work-hardening plastic material, as opposed to an elastic perfectly plastic material. It will be seen that the case of perfect plasticity can be represented as a special case of the work-hardening model. Several basic ingredients are generally employed in developing a constitutive model based on the flow theory of plasticity. Each is listed below and will be described thereafter. * An initial yield criterion, * A hardening rule, * A loading/unloading criterion,

37 * A plastic flow rule. The initial yield criterion may be visualized as an "initial yield surface" in stress space. It may be expressed mathematically as: f(oij) = 0 (2.27) Stress points inside the yield surface represent elastic states and strain increments within the elastic range are assumed to produce stress increments in accordance with generalized Hooke's law. Stress points on the yield surface are said to be in a plastic state. Strain increments in the plastic state are assumed to be composed of a reversable elastic component and an irreversable plastic component. For concrete materials the initial yield (or elastic limit) surface will typically appear as shown in Figure 2.5 in three dimensional principal stress space or as shown in Figure 2.9 in two dimensional principal stress space. Often the initial yield surface is rpresented by a scaled down version of the failure surface. In particular, the failure surface may be constructed from the failure criteria in terms of stress invariants which were described in Section 2.2. Researchers have suggested scaling the failure surface to the range of 30 to 50 percent of its original size in direct correlation to the proportional limit in the uniaxial case (0.30 to 0.50 f/) [27,33]. For an elastic perfectly plastic material the initial yield surface and the failure surface coincide. For a work-hardening material for loading beyond initial yield, the initial yield surface evolves in stress space in the form of subsequent yield or loading surfaces. The current loading surface replaces the initial yield surface in that it defines the new limits of elastic behavior. If the material is unloaded from and reloaded within the current loading surface no plastic deformations will occur until the current loading surface is again reached.

38 Compression 0 Fracture ///// surface/ I definedby // / [ components Chn 32] yield \ j / surface c isotropic hardenloading- Hill 1950; * kinematic hardening - Prager 1955, Ziegler 1959; * mixed hardening - Hodge 1957. surfsurfaces can be expressed as:es Taken froa Chen [32] Figure 2.9: Loading surfaces of concrete in biaxial stress plane The evolution of the loading surfaces in stress space is governed by a hardening rule. Various hardening rules have been proposed and three simple cases are presented hereafter: * isotropic hardening - Hill 1950; * kinematic hardening - Prager 1955, Ziegler 1959; * mixed hardening - Hodge 1957. The isotropic hardening rule assumes that the loading surface expands uniformly during plastic flow as shown in Figure 2.10a. If the initial yield surface has the form F(aij) - k2 = 0 where k is a constant value, then the subsequent yield (loading) surfaces can be expressed as: f(aij, k) = F(aj) - k2(e) = (2.28)

39 Initial yield surface F k2 0 Subsequent yield (loading) surface,F> k2 (-A, a) Isotropic hardening 02 Initial yield surface F =k B Subsequent yield /C o w // )surface.F= k2 0 / (I 0 iu b) kinematic hardening - - \ Translation and m7nI~ \ expansion )Th ) F> k2 Translation only ranslation only F" k2 c) mixed hardening Taken from Chen [32] Figure 2.10: Graphical representation of three simple hardening rules

40 where the scalar k2 > 0 is now a function of effective strain, Ep, which is an integrated increasing function of plastic strain increments. The concepts of effective stress ae and effective strain ep are used in plasticity analyses in order to relate the hardening parameters in the loading function (i.e. k and JP in Eqs. 2.28 through 2.31) to the experimental results from the uniaxial stress-strain curve. That is, the concept of effective stress vs. effective strain relations enables one to extrapolate uniaxial test information into multidimensional situations. The effective stress as a function of effective strain commonly takes the same form as the function of stresses which governs yielding (i. e. F(aij) = ca(ep)). The effective strain is defined as the accumulation of plastic strain increments: p = dep (2.29) where the effective strain increment is often taken as dep = C e^i where C is a constant. As the isotropic hardening rule implies self-similar expansion of the loading surface during plastic deformation it can not account for stress induced anisotropies (e.g. Bauschinger effect). An isotropic hardening rule is therefore not valid for complex loading histories where significant changes in the direction of stress tensor occur. The kinematic hardening rule assumes that the loading surface translates as a rigid surface in stress space during plastic deformation as shown in Figure 2.10b. Mathematically, the kinematic hardening rule has the general form: f(aij, aJ,) = F(,ij - aij)- k2 = 0 (2.30) where cij is the hardening parameter denoting the coordinates of the "origin" of the loading surface which can change during plastic deformation (00' in Figure 2.10b). The hardening parameter acj may be defined in a number of ways, the simplest being

41 as a linear dependence of daij on deij (i. e. daij = cdeij) where c is a work-hardening constant dependent on the material. The kinematic hardening rule is useful in that it is able to account for the Bauschinger effect witnessed in many engineering materials during loading. The mixed hardening rule is simply a combination of the isotropic and kinematic hardening rules (Fig. 2.10c). The loading surface is now able to both translate rigidly and expand uniformly in stress space. It may be expressed in the general form: f(ij, j,, k) = F( - ai) -k2(e) = 0 (2.31) For obvious reasons, the mixed hardening rule provides the best description of concrete behavior for general loading histories. However, the isotropic hardening rule may still be appropriate for loadings which produce stress histories that monotonically increase in magnitude and don't vary significantly in principal direction. For any current plastic state (i.e. f = 0 in Eqs. 2.28 through 2.31) and any specified strain increment loading/unloading criteria are necessary in order to distinguish continuing plastic flow from elastic unloading. Two general procedures defining loading/unloading are listed below: 1. Assume loading occurs and calculate the magnitude of the plastic strain increment, dA (see Eq. 2.34). A positive magnitude indicates loading does occur, while a negative magnitude indicates unloading. 2. Assume that the material has unloaded an infinitesimal amount so that the current stress state resides just within the loading surface. Calculate the elastic stress increment, daij, corresponding to the specified strain increment and test to see whether daf violates the loading surface. If daft stays within the loading

42 surface the behavior of the material is truly elastic; otherwise loading (plastic deformation) occurs. For the case of neutral loading (i.e. stress increments daij lying in the plane tangential to the loading surface) only elastic deformations are assumed to occur. This continuity condition is necessary for a continuous transition from plastic to elastic deformation. Up until this point, the shape and evolution of the loading surface has been outlined while the stress-strain relations in the plastic state still need to be considered. During plastic flow, the total strain increments, deij, are assumed to be the sum of elastic and plastic strain components: dij = de, + dep (2.32) The elastic strain increment de!. can be related to the stress increments through Hooke's law or a suitable nonlinear elastic relation. Using Hooke's law one obtains: dai = Dklde, = DkL(deki - dI) (2.33) The plastic stress-strain relations (or equations of plastic flow) can be developed by introducing the concept of a potential plastic function g(aij, aJ, k), which may be visualized as a surface in stress space. The equations of plastic flow become: deP = dXA (2.34) where dA > 0 is a scalar function of the hardening parameters. The gradient of the plastic potential surface in stress space dX. indicates the direction of the plastic strain increment deP,, while the magnitude of dePj is defined by the scalar function dA.

43 For the. sake of simplicity, lack of experimental data, and certain other practical considerations, the plastic potential surface is often chosen to coincide with the loading surface: g(oj,, k)f = f(,ij, k, ) (2.35) The flow rule may then be stated as: de.3 = dAf (2.36) which is called the associated flow rule since the plastic flow relations are now associated with the loading surface. According to Drucker's Stability Postulate, a stable material is defined as one which: * Positive work is done by any loading acting through the changes in displacements it produces, * Nonnegative work is- done by any complete loading/unloading cycle acting through the changes in displacements it produces. Drucker's stability postulate was mentioned here in that it sets up the inequality relations: daijdeij > 0 for loading (2.37) daijdeij > 0 for complete cycle (2.38) which lead to important conclusions regarding the convexity of the loading surface and the necessity of the normality condition as expressed by the associated flow rule. The complete incremental stress-strain relations for a elastic plastic material based on the flow rule (Eq. 2.34 or Eq. 2.36) may be formulated after one more condition is noted. The consistency condition is necessary to insure that both the

44 initial yield and subsequent stress states satisfy the yield condition f(oi, j k) = 0. That is, both the initial yield function, f, and the subsequent yield function f + df must be equal to zero. It follows that the consistency condition may be expressed as: 09f af af df = dai + dP + Ad = 0 (2.39) This expression, together with Eq. 2.33 expressing Hooke's law and either Eq. 2.34 or Eq. 2.36 expressing the nonassociated and the associated flow rule, respectively, combine to form the constitutive relations for the elastoplastic material: dai = D^kIEkl (2.40) and for the nonassociated flow rule De (De..)( f0,) nD DP = - DC tjftu 8aorr \Oat,, r rskid Dijkl Dk (2.41) where h is a hardening parameter controlled by Of Og Of Ok ag h = _ (2.42) &,5 o,~j ok 04 o~,j The elastoplastic stiffness tensor DiPkl completely defines the incremental stressstrain relations. Note that if the associated flow rule is used in the formulation, then Dekl is symmetric (i.e. substitute f (uij,, ik, k) for g(ari, dki, k) in Eq. 2.41). This has important implications as solving the nonsymmetric problem not only may be several times more expensive but also will require additional data storage area relative to a symmetric problem having the same dimensions. The above brief discussion outlines the premises upon which the ow thereory of plasticity is based. The actual implementation of plasticity based formulations into a finite element code involves additional considerations which will not be expounded upon here. One should refer to [32,33,64,97] for further details.

45 A remarkable number of studies have used the theory of plasticity to model the constitutive behaviors of frictional materials such as concretes, soils, and rocks. The theory is attractive in that it is able to account for the stress history dependent nature of these materials in a relatively simple manner. On the other hand, it has been known for some time that geomaterials don't obey the associated flow rule and it appears the same conclusion applies to concrete as well [88]. Consequentially, plasticity theory based models using an associated flow rule have had trouble predicting certain aspects of the concrete response, such as volumetric strains (dilitancy). Such problems lead some researchers to adopt the more general nonassociated flow rule and the use of 'dilitancy factors' [45]. Another shortcoming of the classical plasticity approach is its inability to account for strain softening behavior (the descending branch of the concrete uniaxial stressstrain curve appears to violate Drucker's stability postulate). Various techniques have been employed to simulate the strain softening behavior within the plasticity theory approach [41,70,46]. Moreover, many researchers believe the plasticity theory is fundamentally unsuitable for modeling concrete materials. The plasticity theory is basically a phenomenological description of the yielding process in metalic materials. On the microstructure level, metalic yield may be described as plastic slip along various crystal planes. Such behavior is possible due to the homogeneity and relatively regular crystal structure of metals. On the other hand, concrete is strongly heterogeneous and has a randomly varying microstructure. While some of the inelastic deformation in concrete may be due to slip type mechanisms (especially under high hydrostatic pressure), a larger portion of the inelastic deformation is a result of microcracking or microfracturing. If the two modes of irreversable rearrangements of the microstructure (i. e. plastic

46 slip and microfracture) were similar, there might be some basis for using the plasticity theory for modeling concrete, but as pointed out in [54] they are profoundly different.

CHAPTER III TENSION CRACKING OF REINFORCED CONCRETE Tensile cracking of the concrete is often considered the major source of nonlinearity in reinforced concrete systems [32]. Excessive tensile stresses in plain concrete cause cracking to develop in a direction roughly perpendicular to the offending prin-. cipal stress direction. It should be noted that this observation pertains to plain concrete only. The presence of reinforcing makes the prediction of what stress levels cause cracking and which directions these cracks form considerably more uncertain. However, most researchers extend the plain concrete cracking criteria to include reinforced concrete as well [88]. When developing analytical models for describing the cracking process in reinforced concrete, three main ingredients are necessary [88]: 1. A criterion which indicates crack initiation, 2. A method for crack representation, 3. A criterion which indicates crack propagation. 3.1 Crack Initiation Crack initiation is usually assumed to occur when the stress at a certain location reaches a prescribed limiting value. The magnitude of this stress value, however, is 47

48 the subject of much debate. Several test procedures have been conducted to ascertain concrete tensile strength. This quantity is often correlated in terms of the concrete compressive strength, fc. One measure of concrete strength, the modulus of rupture fr, may be established by performing simple bending tests of plain concrete beams. The 1983 ACI Building Code [1] recommends the modulus of rupture to be taken as = 7.5gI, though variations in this quantity in the range of 7s to 13i/ have been observed. As previously mentioned, a type of direct tension testing on plain concrete was performed by Kupfer, et al. [56] where the influence of the loading surface was minimized by using steel brush platens. For the uniaxial case, tensile strength values fall inbetween 0.085 and 0.11 times the prism compressive strength. For biaxial stress states, the concrete strength varies as shown earlier in Figure 2.4. A strength criterion based on Kupfer's results, or results from similar ldrect tension testing methods, seems reasonable for situations where nearly uniform strain states are expected over the area in question. A modulus of rupture based strength criterion might be more applicable when considering regions with higher strain gradients as in flexurally dominant situations. 3.2 Crack Representation In general, two basic methods have been applied to represent cracking in reinforced concrete [88]: 1. Discrete cracking models, 2. Smeared (or densely distributed) cracking models. Each of these models has advantages and drawbacks. The choice of one model or the other depends largely upon the problem being investigated.

49 (a) (b) Taken from Chen [32] Figure 3.1: Nodal separation using two and four coincident nodes 3.2.1 Discrete cracking models In a discrete cracking model, cracking is usually represented as a separation of nodes along inter-elemental boundaries (Fig. 3.1). Such an approach requires either the definition of coincident nodes at the onset of program execution or the use of "adaptive" procedures to continually redefine the finite element mesh in the cracking regions. Discrete cracking models are often utilized in a "micro-model" based analysis where details of the localized material behaviors are of primary interest. The discrete cracking model realistically represents the strain discontinuity across dominant cracks in reinforced concrete. The effects of aggregate interlocking can be represented using special link elements that cross the crack and have stiffness both normal and tangent to the crack surfaces. The stiffness provided by these link elements can be decreased as the crack opens, thus providing a realistic representation of the reduction of interlocking forces with increasing crack width.

50 y(z) ay y\' I -. ~~~\ h^~~~~~. Kx(r) (a) (b) Taken from Chen [32] Figure 3.2: Crack pattern and stress distribution 3.2.2 Smeared cracking models Within smeared cracking models it is assumed that the cracked concrete remains a continuum [77]. That is, the effects of cracking are smeared out in a continuous fashion. Upon first cracking, the concrete constitutive models local to the crack region (within the element) become orthotropic, one of the material axes being orientated along the direction of cracking as shown in Figure 3.2. Smeared cracking models are appropriate in a "macro-model" based study where an expedient solution is being sought for the global behavioral characteristics of the system. Some analyses using the smeared cracking approach assume that after cracking occurs all concrete strength vanishes in the direction of principal tension and the cracked surface is also unable to transfer any shear stress. For plane stress situations: dal O O 0 del da72 = 0 E2 0 dE2 ) (3.1) d712 0 0 J d'712

51 where E2 represents the stiffness of concrete in a one-dimensional loading state. Shear strength reserves due to aggregate interlocking can be accounted for by specifying a positive shear modulus, fiG, where the reduction factor fi ranges from zero to one. daf 00 0 dfe dar2 0= 0 E2 0 de2 (3.2) dr72 0 0 /G dyl2 Retention of this shear factor may also be useful in suppressing the real singularity which arises when all elements surrounding a particular node have cracked in the same direction. Without the shear modulus factor no stiffness would be provided to such a node normal to the cracking direction. Incorporating a positive shear capacity also has the effect that secondary cracking need not appear perpendicular to the first direction of cracking. Reasonable results may be obtained using the two constitutive relations when considering the behavior of the structure in its most global sense. However, local to cracking the use of these orthotropic models may not be totally realistic. Actual cracks in concrete will possess rough surfaces and any sliding parallel to the crack direction will tend to cause local stresses and displacements, called "crack dilitancy", normal to the crack surfaces. In order to account for such behavior the off-diagonal terms in the constitutive matrices which relate shear strain to normal stress should also be nonzero with the relative magnitude of these terms decreasing as the crack opening increases [19]. In some ways the smeared cracking model appears more realistic than a discrete cracking method since the nonhomogeneity of concrete tends to spread cracking over a large zone near the crack tip. Fracture is preceded by a formation of microcracks

52 which gradually grow and join into a single dominant crack near the cracking front. Furthermore, the stresses and strains used in any practical reinforced concrete analysis are not the actual stress and strain at a point but rather "the stress and strains in an equivalent homogeneous continuum which smears the inhomogeneous microstructure and averages the properties of the actual material over a certain characteristic volume (called 'representative volume') [18]." In order to reasonably characterize the material, the representative volume must be taken sufficiently large with respect to the maximum size of the inhomogeneities present within the material. Use of the smeared cracking approach also has several practical advantages, including: * Cracking within a finite element code may be represented simply as an adjustment of elemental stiffness matrices. Dicrete cracking models require additional nodes in order to represent the crack not only increasing the number of equations to be solved, but also possibly upsetting any banded nature of the global stiffness matrix. * Cracking may occur in any direction within the finite element model by simply using inclined directions of orthotropy when computing the cracked element material stiffness matrix. Cracks represented by the discrete cracking approach are confined to proceed along inter-elemental boundaries and some adaptive procedure is therefore necessary to represent general cracking behavior. * Cracking criteria are based on the stresses (or strains) measured at interior points of the element, commonly at Gauss integration points. This allows for partial cracking representations (i. e. some Guass points have cracked within the element while others have not).

53 However, when trying to model the behavior of members in which a few large cracks dominate the response, the discrete cracking model may be preferable to the smeared cracking approach. The smeared cracking method has a tendency to diffuse the cracking system in a manner which effectively prohibits any one crack from dominating the response. Such observations should be taken into consideration when investigating shear failure development in some members [86]. 3.3 Crack Propagation Methods for determining crack propagation may generally be classified into one of two catagories, namely: 1) models based on strength criteria, and 2) models based on fracture mechanics concepts. Before proceeding to discuss each of these models, it is instructive to examine the stress conditions which exist in the vicinity of a crack tip. By investigating the stress distributions around eliptical holes which could be degenerated to simulate sharp cracks, Inglis (1913) was the first to show that the stresses near a sharp crack tip are much greater than elsewhere within the stressed body. Such stress fields due to cracking may be conveniently characterized by three independent modes (Irwin 1960) according to whether the displacements contribute to an opening mode (Mode I), an inplane sliding mode (Mode II), or an antiplane sliding mode (Mode III) of the crack surfaces in relation to one another, as shown in Figure 3.3. A general cracking case may be represented by superpositioning any combination of these three modes. Most researchers investigating cracking in concrete have focused on Mode I type cracking, although recently interest has grown in modeling mixed mode behavior. The most direct approach for determining the elastic stress and displacement

54 1 modeI model mode M opening mode sliding mode tearing mode Taken from Broek [26] Figure 3.3: The three modes of cracking fields associated with each cracking mode follows in the manner of Irwvin (1957), based on the work of Westergaard (1939) [72]. The Mode I stress solutions are presented below, referring to Figure 3.4 for notation. Mode II and Mode III solutions have a similar form. A complete derivation of these solutions in presented in [72]. a = Kr c [1-sin sin sin (2ir)1/2 cos [s 2 2 r (2irr) COS + sin 2sin J (3 Kr 09 30 =) sin = cos i cos T ar (2irr) 11/2s 2 2 and r,e = -Ty = 0. The Mode I stress relations may be expressed by the following generalized form: K1 -ij = (2rr)i/2 fii(#) (3.4) KI is the so-called stress intensity factor for Mode I which turns out to be a function of the location and geometry of the crack, the dimensions of the cracked body, and

55 -yz y front z Taken from Paris, et al. [72] Figure 3.4: Three dimensional crack tip coordinate and stress systems the manner in which the external load is applied. Note that as r -+ 0 the stress becomes infinite. The stress intensity factor may then be seen as a measure of the 'stress singularity' at the crack tip. 3.3.1 Models based on strength theory Within this type of an approach crack propagation is based on failure criteria expressed in terms of stresses (or strains). Various failure models for concrete materials were discussed in Section 2.2. Within discrete cracking representations stress levels are typically monitored at nodes along the element boundaries. When the limiting stress levels are exceeded at a given location, the corresponding nodes are separated to simulate the crack propagation. As mentioned earlier, the stresses across the newly formed crack are

56 either immediately reduced to zero or allowed to gradually diminish according to the stiffness of the link elements representing aggregate interlock. If one is using higher order elements, particularly those using isoparametric formulations, simply retrieving the stress values at the corner nodes may provide poor results. Special techniques for improving the corner node stress estimates must be employed, such as extrapolating from values measured at the Gauss points within each element surrounding a given node and then taking the weighted average with respect to elemental area. Within smeared cracking representations stress levels are monitored within the element interiors, typically at the Gauss integration points where stress retrieval is considered accurate [6]. Many smeared cracking models assume concrete to behave linear elastically right up to the point of cracking which is determined by either a maximum principal tensile stress or strain criterion. That is, a crack propagates into the next finite element when the principal tensile stress (or strain) in that element reaches a specified limit. A major objection to the use of propagation criteria based on strength theory is that the resulting analysis is inherently mesh dependent. In reference to the smeared cracking approach, the magnitude of the stress at the front of the crack 'band' is dependent on the finite element mesh size, or more precisely, it is dependent on the width of the front element of the crack band. "The smaller this width is chosen, the larger is the stress in the element just ahead of the front, and the smaller is the applied load which causes crack propagation [14]." In the limit case, one would expect crack propagation at near zero load! The above discussion applies to the strength theory based discrete cracking approach as well. In order to eliminate this crack propagation dependence on mesh size, researchers began investigating the use of fracture mechanics related techniques for describing

57 the cracking process. Some basic fracture mechanics concepts and their application towards concrete materials will be outlined in the next section. It should be noted that recent studies [35] have defended the use of stress controlled smeared cracking methods for situations where the limit strength is not soley governed by cracking (e. g. situations where the amount of reinforcing steel is sufficient to maintain post-cracking equilibrium). The failure of a stress controlled cracking analysis to achieve a stable cracked configuration would strongly indicate the need for a fracture mechanics based approach. 3.3.2 Models based on fracture mechanics concepts "The tensile strength of an idealized crystalline body is the stress which must be applied to cause it to fracture across a particular crystallographic plane [52]." A rough measure of a materials idealized fracture strength may be obtained by modeling the interaction energies of the atoms lying on opposing sides of the fracture plane. Such considerations indicate the ideal fracture strength of many materials to be on the order of E/10, where E is Young's modulus [52]. Why is it then that measured strengths of actual materials lie inbetween E/1000 and E/100? It has long been recognized that the roots of failure at these lower stress levels lie in the presence of pre-existing cracks or other defects in the material. One of the main advances attributed to fracture mechanics is the recognition that the largest defect (e. g. pore, crack) could control the specimen strength and not the overall defect concentration. In light of Inglis work on crack tip stress fields, Griffith (1921) supposed that the small defects present in a seemingly homogeneous material might cause localized stress concentrations reaching the 'ideal' fracture stress. Griffith proposed that frac

58 ture (crack propagation) would occur if the elastic strain energy released per unit increase in crack area was sufficient to provide the increase in elastic energy of the body associated with the formation of the newly cracked surfaces per unit area: dU dS da da where dU/da = elastic strain energy released per unit extension crack area (fixed grip condition), dS/da = increase in elastic energy of the body from the surface energy per unit extension of crack area. The Griffith energy balance concept is the precursor of modern day fracture mechanics. It is limited, however, in that it only considers the resistance to fracture as being provided by the elastic surface energy term (the energy required to form two surfaces) and therefore is only applicable to fracture of very brittle materials. A more general energy balance formulation is necessary to characterize the fracture process of more 'ductile' materials commonly used in engineering structures. Papers by Orowan (1945) and Irwin (1957) proposed that the energy release during fracture in a "quasi-brittle" material was to a large extent dissipated by the production of plastic flow near the crack tip. Tests done by Irwin indicated the energy expended in plastic work necessary to cause unstable crack propagation to be orders of magnitude greater than the surface energy term in Griffith's equations. "However, it appeared that the amount of plastic work in the crack tip region which preceded unstable crack propagation was independent of the initial crack length and was hence as characteristic a measure of the material's resistance to fracture as would be it's surface energy term if it were breaking in a completely elastic manner [52]." Con

59 sequentially, the methods of linear elasticity could still be used to relate the crack tip events to the applied loading. It important to note that these observations were made in the presence of small-scale yielding. That is, the amount of plastic flow at instability was very much smaller in extent than both the crack length and the width of the test sheet. Considerable research has been performed since the early 1960's to determine the applicability of linear elastic fracture mechanics methods to concrete materials. While certain aspects of these research findings are in disagreement, the following general remarks can be made: * Concrete is certainly not a brittle material in the sense that, for example, glass and many other ceramics are brittle. An elastic energy balance approach in the pure Griffith sense is therefore not applicable to concrete. * Concrete materials generally do not possess a unique fracture toughness value. That is, the portion of the energy release rate attributable to irreversable work does not remain constant with varying crack length or specimen dimensions. The generalized linear elastic fracture mechanics approach is therefore not applicable to concrete, except for very large structures which seem to be examples of small scale yielding. * It appears that energy balance concepts are applicable to fracture in concrete provided the irrecoverable work done during the crack formation process is accounted for. Consequently, recent efforts have been directed towards accounting for the work done by microcracking and aggregate interlocking in the fracture process zone during cracking. Some noteworthy contributions to this objective will be briefly outlined.

60 Dugdale model Hillerborg model Figure 3.5: Dugdale and Hillerborg cracking models [48] In [48] Hillerborg and his co-researchers proposed a concrete cracking model similar to the Dugdale-Barrenblatt model for fracture of ductile metals. The Dugdale model (constant closing pressure, a.,, within the plastic zone) and Hillerborg's model are presented in Figure 3.5. In Hillerborg's model, fracture is assumed to occur when the stress at the crack tip reaches the material tensile strength. Here, the same criterion is used for both crack initiation and crack propagation. A crack opening width versus closing stress relation (see Fig. 3.6) is introduced which allows a gradual reduction, rather than an immediate drop to zero, in stress transfer after cracking thus approximately characterizing the behavior of the microcracking zone. At crack width wl, the closing stress has fallen to zero. Because a closing stress acts on the tip portion of the crack during the opening process, energy is absorbed; The amount of energy absorbed per unit crack area in

61 a 0 I.\ ) \ Gn\ W w Crack Width 1 Figure 3.6: General variation of stress with crack width widening the crack from zero displacement up to or beyond a width w1 is equal to fo adw, which corresponds to the area under the crack opening width vs. closing stress curve in Figure 3.6. By choosing the relation in Figure 3.6 such that fo a dw = Gc, where Gc is the critical energy release rate, the cracking model is now equivalent to an energy balance approach. It is noted in [48] that since the precise form of the opening displacement vs. closing stress curve is generally unknown, one may employ mathematically simple models in order to facilitate numerical procedures. In [13] Bazant uses Rice's results [78] for energy change due to crack advance in a continuum as a starting point for the "blunt crack band" model. This model basically modifies the smeared cracking approach so that it represents an energy balance formulation. In its original form, the blunt crack band model assumes the material within the cracked volume to be capable of transferring only normal stress in the direction parallel to the cracking. More recently, the blunt crack band model has been modified to account for the strain softening of the concrete in the microcracking zone [18]. In the presence of reinforcing steel, objectivity in analysis is maintained

62 by realistically accounting for both bond slip and the work done by the release of the bond forces present in the uncracked state. Recent studies using the cracked band model modified to account for strain softening have indicated that the form of the concrete stress-strain relations may significantly influence certain aspects of the response [57]. The following trends were noted: * The post-peak behavior of the concrete stress-strain relations significantly influences the finite element analysis results. The slope of the stress-strain curve immediately following the peak stress is particularly important. * Increasing the fracture energy primarily toughens the specimen response, rather than strengthening it. That is, the ultimate load predicted by the analysis is relatively insensitive to small changes in fracture energy. Despite the theoretical evidence that indicates an energy balance approach is necessary for objective crack representation, many practicing engineers and researchers alike are reluctant to abandon the traditional strength theory based approach (most building codes are based on strength theory failure criteria). Such reluctance may be due, in part, to the fact that stress controlled finite element analyses have had reasonable success in predicting the test results for various types of specimens. Of course, the stress controlled approach is also simpler to implement in a numerical analysis. In [11] Bazant points out that it is important to consider structural size effects when interpreting test results from smaller scale specimens. Bazant has come up with a simple relation associating size effects to fracture toughness for geometrically similar structures. Considering the balance of energy at crack propagation in concrete

63 or rock materials, he hypothesized that the total release of strain energy caused by fracture, W, is the sum of two components: W = W1 + W2 where W1 is the energy released into the fracture front from the remaining part of the structure and W2 represents the energy released from the material in a constant width crack band. The constant width of the crack band, we, has been confirmed by numerous test data. Moreover, the following relation holds: wC = nda where da equals the maximum aggregate size in concrete (or grain size in rock) and n is an empirical constant (n: 3 for concrete and n; 5 for rock). Using both geometrical arguements and energy balance concepts, the following relation for predicting the stress at which the crack band propagates, AN, can be established (see [11] for details): aN = Bft = B + (3.6) (1 + Xl/o)1/2 where A = d/da is the ratio of the characteristic structure size d relative to the aggregate size da (A > n = 3 for concrete) and both Ao and B are constants for geometrically similar structures. Under the hypothesis that the total potential energy release W caused by the fracture is a function of both: 1. The length, a, of the fracture (crack band), 2. The area of the cracked zone, nd,, a dimensional analysis likewise leads to Equation 3.6. The same result is also confirmed in [98]. "The dependence of W on length a characterizes the energy that is released into the fracture front from the surrounding uncracked structure, while the dependence of W on nda characterizes the energy released from the cracking zone [11]." It is

64 z -- E -Strength or yield | fi t^rfr,> s. criterion c -^,,- Linear elastic X ' A h fracture mechanics Most existing tests Nonl U3 wNonlinear -Wi~~ sfracture mechanics 2e. 2 0 z log(size) Figure 3.7: Size effect for plain concrete structures [11] quite interesting to note that if the W depended only on a one would get the well known size effect from linear elastic fracture mechanics (i. e. aN = CVdi where C is a constant). If W depended only on nda one would obtain a simple strength criterion (i. e. aN =constant). As Equation 3.6 would predict, the true size effect in concrete structures may be represented as a gradual transition from a strength criterion for smaller structures to a linear elastic fracture mechanics estimate of the fracture strength for larger structures (Fig. 3.7). The size effect relation indicates that most of the fracture test results were obtained from specimens in the 'strength criteria' governing size range, while also indicating linear elastic fracture mechanics may become applicable for large enough structures (e. g. dams, nuclear power vessels). In light of such information, one sees that generalizing test results obtained from reduced scale concrete specimens to much larger real structures may be erroneous and in fact dangerous. The 'safety margin' may be smaller than thought for large

65 log N 2 N 2 1 1 _____ N _ N__ _ strength criterion reinforced plain yield no yield log Figure 3.8: Size effects for reinforced concrete [11] structures built using current building codes in situations where brittle fracture is possible. The same type of size effect relation is found when reinforcing steel is present at the fracture front provided the steel has not yielded (Fig. 3.8). The effect of reinforcing steel is to shift the asymptotic approach to the LEFM solution to the right. That is, the validity of a strength criterion is extended and the transition to a LEFM criterion occurs at a larger structural size. If the reinforcing steel yields, the strains become large to the extent that the concrete is no longer providing any tensile resistance and the loading ends up being resisted by the reinforcement alone. In such a case, the value of log aN again approaches a constant. A few comments are now presented on statistical size effect in concrete structures [11]. Due to concrete's heterogeneous nature, concrete strength varies randomly throughout a structure (neglecting strength dependencies on fabrication process, etc.). While this variation in strength is independent of structural size, however

66 the stress gradient usually varies inversely with structure size (e. g. beam in bending) and therefore the regions of maximum stress become larger in a larger structure. It follows that the chance of a highly stressed region containing a low strength sector increases with larger structural size. That is, apparent strength decreases with an increase in structure size. However, this trend eventually diminishes as the region of near maximum stress becomes much larger than the regions of low strength.

CHAPTER IV SOLUTION TECHNIQUES FOR NONLINEAR PROBLEMS The material in this section attempts to highlight some of the basic techniques and considerations related to solving the nonlinear problem. At the expense of some generality, the discussion will be limited to quasi-static loading conditions and material nonlinearity only type problems. Dealing with problems where inertial effects, geometric nonlinearities, or material-time dependencies are present generally involves a straitforward extension of the basic concepts presented hereafter. 4.1 Fundamental Concepts The analysis of path dependent materials generally requires an incremental application of loading (or displacement) since the material properties, and thus the structural response, depend on stress and strain history. Often, the solution procedure assumes a piecewise linearization of the nonlinear behavior over each increment. The basic approach in formulating an incremental (step-by-step) solution procedure is to assume that the solution at time t is known and that the solution at time t + At is required, where At is a time increment suitable for the problem at hand [6]. In general, the greater the nonlinear behavior, the smaller the time increment chosen. Note that even though the problem is static in nature it is convenient to refer 67

68 to the variation in external loading and the corresponding structural configuration as a function of time. The simplest of the incremental methods for solving the nonlinear problem is Euler-Cauchy forward integration. The solution procedure may be expressed in the following form: tKU R (4.1) where tK represents the tangent stiffness matrix corresponding to the structural configuration at time t, U is the vector of nodal point incremental displacements (i.e. U = t+U - tU), and R is the vector of incremental externally applied loads (i.e. R = t+tR- tR). However, the change in material properties over the time increment At is not reflected and the response will generally appear stiffer than it actually is. This solution may rapidly drift away from the true solution path since no check is made to insure structural equilibrium against the external loading after each load increment. In order to estimate the accuracy of such a procedure one would have to perform repeated analyses using successively smaller loading increments; a reliable solution becoming quite expensive. The equilibrium relations for any time t may be expressed by: tR- F = 0 (4.2) where tR is the vector of externally applied nodal loads which, in general, is the sum of several components: tR tR+ tRs + tRc (4.3) where the subscripts B, S and C denote the effects of body forces, surface forces, and concentrated loading, respectively (see [6] page 124 for further details).

69 tF is the vector of nodal point forces corresponding to the element internal stresses and may be computed according to: tF = E ( tB(m)T tr(m) tdV() (4 4) where the summation is performed over all elements, tB(m) is the strain displacement matrix for element m, tr(m) are the stresses in element m, and tV(m) is the volume of element m (for small displacement analyses tV(m) = Vm)). Note that during the solution process the equilibrium conditions expressed by Eq. 4.2 generally can not be satisfied exactly and one is left with an unbalanced residual load vector tJ computed as the difference between the externally applied loading and the loading corresponding to element internal stresses: t = tR - F (4.5) In summary, in order for nonlinear analyses of path dependent materials to provide meaningful results the solution technique should satisfy the following requirements: 1. The load is applied incrementally in order to capture the path dependent nature of the material. The load step size should reflect the degree of nonlinearity in the response. 2. Equilibrium is sought for each accepted structural configuration along the solution path (i.e. at times t, At, etc.). Generally, equilibrium as expressed by Eq. 4.2 can not be satisfied exactly and one should therefore think in terms of bounding the unbalanced residual forces ' within certain limits.

70 4.2 Solution Techniques One means of reestablishing equilibrium when using the Euler-Cauchy forward integration scheme is to simply add the current unbalanced residual forces to the next load increment. This method is attractive due to its simplicity but there is no guarantee that the true deformation path is being followed [20]. A more elegant, albeit more computationally expensive, approach is to iterate within each increment to find the equilibrium state. These equilibrium iterations are quite important when analysing reinforced concrete structures. For example, after applying a load increment the elastic energy released by cracking in a certain region must be redistributed out to the surrounding elements which may cause further cracking. Equilibrium iterations which properly converge allow cracking to stabilize before moving on to the next load increment. The stabilization of the cracking process during iteration is shown in Figure 4.1, where the data has been taken from the shearbeam example presented in Section 5.4. Two popular techniques commonly applied to nonlinear reinforced concrete analyses will be discussed: 1. Newton-Raphson iteration 2. Modified Newton-Raphson iteration Assume that the true solution to the equilibrium conditions (Eq. 4.2) at the end of the load increment may be represented as U*. The state of equilibrium may then be expressed by: i(U*) = t+AtR(U*) - t+tF(U*) = 0 (4.6)

71 12 1 10 8 0II 4o 2 0 2 4 6 8 10 12 14 I terotion Cycle Figure 4.1: Cracking during iterative process By iterating we wish to successively improve test solutions U(i) for (i = 1, 2...) until the test solution at say the nth iteration, Un), is sufficiently close to the true solution for certain convergence criteria to be met. The convergence criteria play an important role in determining both the validity and efficiency of a nonlinear analysis where equilibrium iterations are necessary. If the convergence tolerances are too relaxed the true conditions might not be accurately represented (e.g. accepting convergence before strain energy from cracking elements is allowed to redistribute as in Fig. 4.1). On the other hand, convergence criteria which are too strict may cause undue computations to be performed. In general, the solution obtained for each step should be tested to see whether it has converged within specified tolerances. It is also important to monitor, or simply limit, the iteration sequence in case the solution process begins to diverge. Further considerations regarding convergence criteria are given in [6].

72 If a first term Taylor expansion is used to evaluate successive improvements in the displacement solution, we get the following expressions which represent a NewtonRaphson solution of the equilibrium relations for time t + At in Eq. 4.2 (see [6] for derivations): t+AtK(i-l)AU(i) - t+atR - t+tF(i-1) (47) t+atU(i) - t+tU(i-1) + U(i) (4.8) where t+AtK('-l) is the tangent stiffness matrix for iteration (i 1). The initial conditions to the iteration sequence are t+AtK(o) = tK, t+tF() = tF, and t+tU(o) = tU. Note also that the external loading on the structure is assumed independent of displacement and is therefore independent of the iteration step. It is convenient to use a one-dimensional graphical representation of the structural load-deflection behavior when discussing these iteration techniques. The structure tangent stiffness at any point on the load-deflection response curve may be represented as the slope at that point. The essential character of the Newton-Raphson technique may be represented graphically in one dimension as shown in Figure 4.2b. It should be pointed out that the structural equilibrium is actually formulated in multidimensional space and that the gradient of such expressions may be thought of as a hypersurface within that space. While the Newton-Raphson method often possesses excellent convergence characteristics (quadratic convergence rate), however, the tangent stiffness matrix must be updated and triangularized each iteration step. A modified Newton-Raphson technique effectively reduces the number of stiffness updates and factorizations during the solution process. Commonly, the tangent stiffness matrix is calculated and reduced only for the initial iteration step of each load increment and then held constant during subsequent iterations (Fig. 4.2b). Therefore, each subsequent iteration requires

73 R R U U Initial Stiffness Tangent Stiffness Figure 4.2: Iteration techniques represented in one dimension only an evaluation of the unbalanced residual forces and then a back-substitution to arrive at the next test solution. Alternatively, it might be advantageous to hold the tangent stiffness matrix constant over consecutive load increments if the system shows little nonlinearity in its response. The modified Newton-Raphson iteration may be expressed by: KA(i) = t+AtR - t+AtF(i-1) (4.9) with r representing one of the accepted equilibrium conditions (e. g. at time t, t+ A t, etc.). The interval at which the stiffness matrix is updated and triangularized should depend upon the degree of nonlinearity in the system response. While this technique will require more iterations to converge relative to a pure N-R method, each iteration itself will demand much less computational effort, especially for problems containing a large number of equations. Along with the trouble associated with the expense involved in forming and solving the tangent stiffness, for certain problems there are several other shortcomings inherent in the N-R method [84]. As a result, numerous other methods for solving

74 nonlinear problems have been proposed but have not been reported here due to both time and space limitations. However, a few techniques for increasing the efficiency of the methods outlined are worth mentioning: * Employing a so-called "acceleration scheme" within the modified N-R method which scales the displacements during iteration in a manner which simulates full Newton-Raphson convergence [64]. * Automatically computing the load increment and direction within the program (load step generation) according to a "current stiffness parameter" [20]. * Adjusting the progress of the computations along the solution path according to some form of "arc length" control [34,79]. By specifying an arc length parameter (geometrically, an approximation to the length of the step in loaddisplacement space), an auxiliary equation is introduced into the analysis relating the load and displacement increments. That is, the size of the load increment becomes mutually dependent on the size of the displacement increment, or in a sense, the stiffness of the structure (a more flexible structure will demand smaller load increments). * Using a quasi-Newton method which involves updating the stiffness matrix in a simpler way after each iteration, rather than recomputing it entirely (full N-R) or holding it constant (modified N-R) [7]. The idea is to update the inverse of the coefficient matrix in order to provide a secant approximation to the matrix during the iteration step.

75 4.3 Linear Equation Solver As previously mentioned, the solution of the nonlinear problem can essentially be approximated by the step-by-step solution of incrementally linear relations. For the sake of simplicity let them be represented by: KU=R (4.10) where K is a structural stiffness matrix (in most cases symmetric positive definite), R is a load vector, and U is a vector of unknown displacements. As Eq. 4.10 is solved repeatedly in nonlinear problems, rather than just once in a standard static linear analysis, the importance of employing an efficient linear equation solver is paramount. In [61] it was estimated that approximately 30 to 50% of the computing time in the standard linear analysis is devoted to linear equation solving, while for the nonlinear case the time spent on linear equation solving will most likely be even greater. Indeed, for the building analysis presented in Chapter VI roughly 75% of the solution time was spent factorizing the coefficient matrix. Such observations fostered considerable research directed towards developing efficient equation solving routines during the 1960's and early 70's. Another concern addressed by program developers during this time was the need to develop solution techniques capable of solving very large problems with only a limited core storage capacity at their disposal. 4.3.1 Basic methods Basically, the procedures for solving Eqs. 4.10 may be classified into one of two methods, namely: 1. iterative solution techniques 2. direct solution techniques

76 Iterative solution techniques were popular during the developmental stages of the finite element method, the Guass-Seidel method being an example of an iterative solution technique. Typically, an initial estimate of the true solution vector is improved successively by the solution algorithm until convergence is obtained within a preset tolerance. The number of iterations required to arrive at a solution varies according to the choice of initial estimate, the convergence tolerance, and the condition of the coefficient matrix. Relaxation factors may be used to speed up convergence. Iterative solution techniques have been, for the most part, displaced by the more efficient direct solution techniques described next. Iterative solution techniques may still be found useful in certain instances, such as when performing step-by-step analyses when the structural conditions don't vary significantly. Here the converged solution from the previous step provides an excellent estimate of the solution to the current step. While many direct solution techniques have been developed which appear to be unique, almost all successful techniques are based, in theory, on standard Gaussian elimination. In fact, nearly all modern day large scale finite element codes use some form of Gauss elimination as an equation solver [6]. Rather than finding the inverse of the coefficient matrix, K-'1, which would be extremely inefficient, the coefficient matrix is decomposed into the following 'triangularized' form: K = LDLT (4.11) where D is a diagonal matrix, L and LT are lower and upper triangular matrices, respectively. The load vector R of Eqs. 4.10 is reduced using L; the result of which leads, through backsubstitution, to the solution U. A thorough investigation of the Gauss elimination procedure is presented in [87].

77 It is well known that judicious nodal numbering of a finite element model will result in a reduced bandwidth for the coefficient (i.e. stiffness) matrix. A substantial savings in computing time can be realized by modifying the standard Guass reduction algorithm so that the two innermost loops only count up to the limits of the maximum half-band width. Further savings can be obtained by adjusting the loop indices to reflect the variable row-by-row band width. In addition to reducing the number of computing operations, significant reductions in data storage are realized by storing only those terms within the maximum half-bandwidth (i.e. data is stored in an N by M array where N is the number of equations and M is the maximum half-band width). Still further savings in both computing operations and storage are realized if the solution procedure proceeds in a column-wise fashion. The modified Crout algorithm presented in [62] uses a column-wise solution sequence and virtually eliminates all unnecessary operations. Storage can be optimized by saving only the elements within the matrix 'skyline' into a one dimensional array (called "compacted column storage"). An additional pointer vector is necessary for indicating the location of the first nonzero element of each column. A similar algorithm, called the "column reduction method" is presented in [6]. The frontal solution method is another technique which is based on Gauss elimination, but here the coefficient matrix assembledge and the elimination process are not so clearly separated as in most other solution techniques. The assembly of the stiffness coefficients and the reduction process are performed in tandem so that at any one stage in the solution process there is a group of "active" degrees of freedom (the wavefront) awaiting static condension. The solution process procedes elementby-element through the structure adding new degrees of freedom to the wavefront

78 while others are being condensed out. The active degrees of freedom are condensed out only after their corresponding stiffness terms are fully summed. All data is kept in auxilary storage and brought into core only when needed. With effective element numbering the size of the wavefront is kept to a minimum and relatively large problems can be solved without resorting to an out-of-core solution technique. One particular advantage of the frontal method is its ability to handle the additional elements with relative ease as the size of the bandwidth is inconsequential. This would prove useful during mesh refinement of a certain region of the model. The main disadvantage of the frontal solution technique is the large amount of both program and computational overhead required to keep track of the active degrees of freedom and the reduced stiffness coefficients. If the wave front becomes large an out-of-core solution would be necessary and the bookkeeping becomes even more challenging. 4.3.2 Coping with limited storage availabilty Traditionally, the computer has had two basic areas in which it can store data. Data can be stored in "high speed" memory locations directly accessible by the computational processing unit (in core) or on auxiliary storage devices such as tapes or disk files (out-of-core). Data stored on the auxiliary, or external, devices must be read into the high speed memory before it can be used for computational purposes. The auxiliary storage devices are often called "permanent" storage while the high speed storage is often referred to as "temporary'. Both physical and practical considerations put limitations on the size of the high speed storage area. Expanding the high speed storage area in a cost effective manner is one of the ongoing areas of research in computing technology. In contrast, the

79 storage capacity of the auxiliary devices is quite sufficient for most user needs. For many finite element problems the number of simultaneous equations to be solved is quite large (even after considering the storage saving schemes previously outlined) thus prohibiting a complete in-core solution of the linear equations. Program developers worked around this limitation by going to an out-of-core solution, that is, using the external storage as a database and then moving "blocks" of data in and out of the high speed storage as necessary. Presently, even a microcomputer is able to perform large scale finite element computations, the only physical limit to problem size being imposed by the capacity of the disk or tape. One stiff penalty comes with the increased capacity of an out-of-core solver, namely the extreme inefficiency of input/output operations to and from external storage relative to the processing speed of the arithmetic unit. Significant savings in computing time can be realized by optimizing the block dimensions to reduce the number of I/O calls. In [28], a simple block-solver procedure is developed noting that the Gauss elimination procedure is valid whether the elements of the stiffness matrix [K] are individual coefficients or square sub-matrices of coefficients. A large coefficient matrix of individual elements can therefore be transformed into a equivalent matrix of smaller dimensions where each element is now a square block of individual elements. During any stage of the solution process, Gaussian elimination requires no more than three blocks of coefficients to reside in core. Knowing the capacity of the high speed memory, it follows that the block dimensions can be made as large as possible to reduce I/O costs. It should be noted that this technique, as presented in [28], requires the inversion of the block matrices along the diagonal. In [93], the band solver technique discussed in Section 4.3.1 is extended to perform

80 out-of-core operations. The coefficient matrix and load vectors are partioned so that all blocks contain the same number of equations per block, the number per block being chosen according to the available core storage and the maximum half-band width of the coefficient matrix. The SAP program uses this type of out-of-core method [93]. In [63], the modified Crout routine is extended to perform out-of-core solutions. A complete description of the partitioning methods. and a source code listing of the program is also given in [63]. Using the compacted column storage format results in more efficient use of the space within each block, thus reducing the rintotal number of blocks required for solution. It follows that both the number of I/O operations and the overall computing time to reach a solution will be less relative to the equation solvers based on the banded storage method. The ADINA program uses an out-ofcore solver similar to this one [6]. More recently, computers with "virtual memory" capability have become widely available. The limitations on the core size still exist, however the programmer does not see this restriction as the computer is able to "page" the data to and from core as necessary. This not only relieves the programmer from the burden of coding the outof-core algorithm, but also results in much faster execution times as the computer is using more efficient data transfer stratagies than the standard I/O routines.

CHAPTER V THE SNAC PROGRAM The SNAC finite element program was developed primarily for use as a tool to investigate the inelastic response of shear wall dominant buildings subjected to quasi-static loading histories. The program is based on the finite element displacement formulation which is well documented in both textbooks and literature [6,96]. Loading in the form of controlled nodal forces and/or displacements is applied incrementally with iterations within each increment to reestablish equilibrium. The program formulation allows for realistic representation of the material nonlinearities which are essential in an ultimate analysis of reinforced concrete structures. Nonlinearity is limited to material response where the models have been taken from existing literature. SNAC is coded in standard FORTRAN 77 language and was developed on Apollo computers within the Computer Aided Engineering Network (CAEN) at The University of Michigan. Program oraganization and indeed much of the coding was taken directly from DLEARN [50], an educational finite element program for dynamic linear analysis of structures. The spirit of DLEARN has been maintained in that SNAC remains in a form instructional for student use and modification. Clarity in program logic has been stressed, possibly sacrificing efficiency in certain instances. 81

82 5.1 Program Description Upon executing the SNAC program, the routine MAIN sets the total memory size allocated to data array working space and initializes various input and output files. The routine ADRIVE is then called which, in effect, drives the program global operating sequence. The essential features of routines MAIN and ADRIVE are illustrated in Figure 5.1. At this global level, the program logic is quite similar to that of DLEARN [50] and the reader is referred to this document for more specific details. For large problems nearly 100% of the computing time is spent in routine NONLIN which drives the nonlinear solution sequence. A somewhat simplified version of routine NONLIN is outlined in Figure 5.2, where: NSEQ: number of time sequences, NSTEP: number of load increments within a given time sequence, ITRMAX: maximum number of iterations allowed. All other notation and concepts presented in Chapter IV apply here, except that reference to time in the equations has been shifted by a constant, At. That is, the last accepted equilibrium configuration is now at time t - At and we are seeking a converged solution at time t. Some options available within the SNAC program have been omitted in order to improve the readability of the flow diagram. 5.1.1 Element library The current version of SNAC supports three element types, namely: 1. QUADC - Two-dimensional continuum element, 2. TRUSS - Two node truss element, 3. BEAM - Two node beam-column element.

83 Routine MAIN * dimension total memory size * initialize input/output files * call ADRIVE (program global driving routine) * close input/output files Routine ADRIVE Input global control information I Initialize time sequence variables Set memory pointers for data arrays Input nodal data Input element data i Allocate memory for global equation system I Call NONLIN (solution sequence driving routine) I Output memory pointer directory information 1 Print elapsed time summary TRUE I CONTINU RETURN to MAIN Figure 5.1: Main and global driving routines

84 DO nsq = 1, NSEQ I — Set time sequence parameters: NSTEP, At I DO n= 1, NSTEP t= t + Set external nodal force vector 'R Set imposed displacements tUD _ _ _I ( DO i1,ITRMAX > i Update internal nodal force vector 'F tF(-l) = tF(i-l) + AF(- all almn Determine residual force vector t t(i-l) = tR - tF(i-l) Check for convergence IF (f(tp(i-l)) < TOLERANCE) CONVERGE = TRUE IF CONVERGE TRU Form tangent stiffness matrix 'K tK,(-) = E AK-I) al elem Solve for incremental displacements AU tK(O-)A0i) = t.(i-l) Update total displacements 'U 'U(i) = tu(i-1) + AU(i) I CONTINUE!! Output nodal displacements Output element stress/generalized forces I CONTINUE! C" co..u RETURN~DRo Figure 5.2: Nonlinear solution sequence driving routine

85 Each of these element types is based on standard finite element formulations and possesses nonlinear material capabilities. Further details concerning these three element types are presented in Sections 5.2 through 5.4, respectively. All calls to the element routines are made through an element library routine. Operations performed by the element routines are conveniently separated into element tasks. The SNAC program divides the element calls into four tasks, namely: 1. INPUT EL Input element data, 2. ADD FINT Add element contributions into the internal nodal force vector, 3. FORM LHS Form the left-hand-side matrix to the linearized equations, 4. STR PRNT Print element stress/generalized forces to output files. In this manner, new elements may be added to the existing library with relative ease. 5.1.2 Load application The two outer loops of the flow diagram in Fig. 5.2 may be characterized as the "advancing phase" of the solution process. Two loops are provided rather than one to allow flexibility in varying the load step size during the solution. Loads may be in the form of external nodal forces and/or nodal displacements. Loading is input in the form of nodal load vectors. Nodal load vectors vary in magnitude during the solution history according to prescribed load-time functions. Load-time functions are input as discrete values over time with piecewise linearity assumed between data points. At each time t, the nodal load vectors (scaled by the appropriate loadtime functions) are superimposed to determine the external nodal force vector, tR, and the imposed displacements, tUD. Values belonging to either tR or tUD within each load vector are differentiated according to the restraint conditions specified for the corresponding degrees of freedom. Virtually any static loading history can be

86 applied to the structural model through proper selection of the nodal load vectors and corresponding load-time functions. The inner loop of the NONLIN routine flow diagram may be characterized as the "correcting phase" of the solution process. Due to material nonlinearity an imbalance between the external loading and the internal nodal forces is created, as represented by the residual force vector ti. Equilibrium iterations are then used to bring this imbalance within acceptable tolerances. The iterative procedure shown in Fig. 5.2 corresponds to a full Newton-Raphson technique where the structural tangent stiffness matrix is formed and factorized each pass through the loop. The actual implementation of this process allows the user to specify the interval at which the tangent stiffness matrix is updated and factorized. Should the correcting phase fail to converge within a specified number of iterations the solution procedure is terminated. Nodal forces incremented in the advancing phase are held constant within the correcting loop (i.e. force controlled loading). The corresponding displacements will typically increase in magnitude during iteration for a softening type response. In general, the solution procedure will reach a state where the numerical model is unable to withstand the incremented load level and the correcting phase will fail to converge. Similarly, nodal displacements incremented in the advancing phase are held constant in the correcting phase (i.e. displacement controlled loading). The corresponding forces will vary in magnitude during iterations, typically decreasing for a softening type response. This type of loading scheme allows the tracing of structural response beyond the peak load level. As briefly mentioned in Section 4.2, more sophisticated solution techniques are available for dealing with general structural response characteristics such as softening, stiffening, "snap-through" and "snap-back".

87 5.1.3 Material only nonlinearity The current version of SNAC is limited to material-only nonlinearity in the structural response. The assumption that kinematic (geometric) nonlinear effects need not be accounted for is reasonable for many types of reinforced concrete analyses, especially the case of shear wall dominant structures where displacements prove to be limited in magnitude. Disregarding the kinematic nonlinear effects implies that the current configuration of each finite element at each stage of the loading history is assumed to coincide with the respective finite element original configuration. This has important practical implications, namely: 1) the integrations required in the formulation of the element stiffness matrices and load vectors are performed over the original volumes of each respective finite element, and 2) the strain displacement matrix of each element remains constant, independent of the element displacements. These two aspects of the finite element formulation are therefore the same as in a linear analysis. Again, the nonlinearity in the stiffness formulation during loading comes from the material behavior as expressed by the constitutive relations. Causes of reinforced concrete nonlinear behavior include: 1. nonlinear behavior of concrete under multiaxial stress states, 2. tensile cracking of the concrete, 3. nonlinear behavior of the reinforcing steel, 4. bond slippage at the reinforcing steel-concrete interface, 5. doweling action of the reinforcing steel across cracks, 6. stress transfer across cracks due to aggregate interlock.

88 With the exception of Item 4, all the above items are modeled either directly or indirectly within the continuum element modeling of the shear walls. Modeling of bond slip is desirable, and in fact necessary for an objective modeling of the structural response, but has not been included here as perfect bond is assumed to exist at the concrete-steel interface. Note that bond slippage and the associated degradation in stiffness may be a critical factor contributing to the structural reponse in cyclic loading situations. However, the material models in this study have been chosen to reflect the aspects of nonlinearity considered essential for quasi-static monotonic loading conditions. The first three items listed above are indirectly modeled within the beam-column element material subroutines. No direct modeling of these items was possible as all nonlinearity is assumed concentrated in zero-length hinges located at the element end nodes. Hinge properties were determined through an integrated section analysis where certain assumptions regarding the loading pattern and damage distribution over the member length had to be assumed a priori. 5.1.4 Convergence criteria Convergence criteria play an important role in determining both the validity and efficiency of a nonlinear analysis where equilibrium iterations are necessary. If the convergence tolerances are ineffective or too relaxed the true conditions might not be accurately represented. On the other hand, convergence criteria which are too strict may cause excessive computations to be performed. The SNAC program performs two types of convergence check each pass through the correcting phase loop of the solution process, as suggested in [5,6]. First, the euclidean norm of the current residual force vector, tq(t-l), is required to be within

89 a specified tolerance, EF, of the imbalance initially created by the load increment: It|(0)) - 2<eF (5.1) where we are seeking a converged solution at time t. In certain instances, convergence of the out-of-balance forces may not provide a reliable indication of whether the associated displacements are sufficiently close to their equilibrium values. For example, consider an elastic-plastic type material with small strain hardening modulus becoming plastic during a load increment. Seemingly acceptable imbalances in force terms may correspond to quite large discrepancies in displacements. These considerations have been addressed by requiring the increment in internal energy over each iteration (i.e. the amount of work done by the residual forces acting through the associated displacement increments) to be within a specified tolerance, eE, of the initial internal energy increment. AU(i)T tgI(i-1) AU(1)T t() (5.2) If the normalized internal energy increment, as expressed by the left hand side of Eq. 5.2, becomes greater than unity the user is warned of possible solution divergence since the equilibrium iterations are creating energy in the numeric model. Upon completion of the equilibrium iterations SNAC writes the final convergence data, including the number of iterations required, the normalized norm of the residual force vector, and the normalized internal energy increment. 5.1.5 Linear equation formation/solution The linear set of equations is solved by a compacted column solver technique based on the Crout method of elimination [50]. Such an approach attempts to

90 optimize both the number of operations involved and the storage size of the stiffness coefficient matrix. 5.1.6 Constraint equations Various kinematic boundary conditions can be realized by specifying a set of constraint equations for the degree of freedom in question. The constraint equations are enforced at the element level during assembly of the stiffness matrix and load vector. Constraint equation number j takes the form: Uj = aiui i = l...n (5.3) where uj is the dependent degree of freedom, ui are the independent degrees of freedom, ac are multiplier coefficients, and n is the number of terms in the constraint equation. The ability to enforce displacement constraints in this general manner proved invaluable during the building analysis. For example, constraint equations were used to insure compatable displacements when transitioning the QUADC elements from a coarser to finer mesh, provide stiffness to the rotational degrees of freedom where the BEAM elements connect to the QUADC elements, and specify building lateral displacements according to rigid diaphram assumptions. 5.1.7 Graphics post-processing According to user request SNAC writes information to several graphics postprocessing data files. A menu driven interactive graphics program was written in Pascal language using the GPR (Graphics Primatives) routines available on the Apollo system. After loading the graphics program, the user can access any data set saved from a previous SNAC run. One can then advance through the solution sequence step-by-step to follow deformations, principal stress trends, or nonlinear

91 trends (cracking and plasticity) by pressing any one of the three mouse keys provided on the system. By pressing the rightmost mouse key the user creates a Post Script language version of the screen graphics which can be sent to a laser writer for printing. The deformation plots and plots of cracking trends seen later in this report where created in this manner. 5.2 QUADC - Two-dimensional Continuum Element The QUADC element supported by SNAC is basically a version of the twodimensional isotropic elasticity element presented in DLEARN [50] which has been modified to account for material nonlinearity. Element stiffness matrices are based on standard isoparametric finite element formulations (i.e. both the element coordinates and the element displacements are represented using the same shape functions which are defined in a natural coordinate system). The QUADC element may be used in triangular (three-node) or quadrilateral (four-node) form with the latter case used exclusively in this study. These elements are often referred to as "linear" or "lowerorder" elements as interpolations are based on linear or bilinear shape functions for the three-node and four-node cases, respectively. The reader is referred to [50] for further details regarding the properties and formulation of these elements. Triangular or quadrilateral linear elements have been used in a variety of reinforced concrete planar analyses, including [13,17,30,36,57,95]. The choice of lower order elements seems reasonable considering the types of nonlinearities involved in reinforced concrete modeling. Physical discontinuities, such as concrete cracking and localized crushing, are represented by abrupt changes in deformation in the finite element model. The use of higher order interpolations would imply artificial differentiability and smoothness of the displacement functions [20].

92 Consistent with a smeared modeling approach, material nonlinearity in the QUADC element is simulated by adjusting the element material response matrix coefficients at the appropriate integration points. Each QUADC element is identified with one of three material types as specified on its "material card" during program input. The first material type allows for linear elastic isotropic material behavior as defined by Young's modulus, E, and Poisson's ratio, v. The remaining two material types are considerably more involved and are discussed in the following sections. There is also an option (independent of material type) which allows for the smeared representation of a steel reinforcing mesh within the concrete material. The user is able to specify the reinforcing ratios of the steel mesh in the two orthogonal directions aligned with the global coordinate axes. Reinforcing mesh material nonlinearity is limited to elastic-perfectly plastic behavior in each of the two orthogonal directions. The reinforcing mesh is assumed to be incapable of resisting shearing stress. 5.2.1 Non-orthogonal cracking model w/tensile softening This model, hereafter referred to as the "non-orthogonal cracking model", was based on the work presented in several recent publications [24,25,81]. Some of the essential features of this model will be briefly outlined and the interested reader is referred to these documents for more detailed information. The smeared cracking approach treats cracked concrete as a continuum material, as if the local discontinuities due to cracking were distributed evenly over a certain tributary area of the finite element (see Section 3.2.2). Traditional approaches [77, 85] assume a mixed stress-strain relation for the cracked material (i.e. both the intact concrete and the discontinuities obey the same stress-strain relations and are

93 l Onn ann 'nn co/6s~~o nEnn ann nn Figure 5.3: Decomposition of total strain [81] therefore "indistinguishable"). The striking feature of the non-orthogonal cracking model is that it treats the concrete constitutive behavior separately from the crack interface behavior. This is accomplished by resolving the total strain, E, within the fracture zone into concrete strain, eco, and crack strain, em. Fig. 5.3 illustrates this type of resolution for strain components normal to the cracking direction. The strain, Ec, which represents- the opening of the microcracks, acts only over a limited region (i.e. the width of the fracture zone or the "crack band width" in the numeric model). Cracking strains should therefore reflect a certain dependence on the width of the fracture zone. Since the solution process proceeds incrementally, the total strain increments are resolved into concrete strain increments and crack strain increments. AC = ACo + Ac' v(5.4)

94 The usefulness of such a decomposition becomes apparent after allowing for further sub-resolutions in the concrete strain and crack strain fields. For example, the concrete strains may be composed of elastic, plastic, creep, and thermal strains [24]. Likewise, the crack strain field may be composed of the contributions due to a number of individual cracks which act concurrently at an integration point. In the above approach, nonlinear actions in the intact concrete may be combined with multiple non-orthogonal cracking in a systematic manner. Let Aect denote the vector which assembles the individual crack strain increments: Ae = (Ae AyF Aer A-I... Aecr AYcr)T (5.5) where Aecr and Ay^c are the normal and the shear crack strain increments for crack number n, respectively. The number of elements contained in Aec varies according to the number of open cracks, n, at an integration point. The SNAC implementation of this model is limited to two simultaneous cracks, although this restriction is not essential. Individual crack strain increments are related to the global crack strain increments through a transformation matrix, N, which is given below for a two-dimensional material idealization: Ac = N Ae' (5.6) cos2 81 - sin 81 cos 01... cos2 On - sin On cos On N = sin2 81 sin 0l cos 01... sin" 8n sin On cos On 2 sin 81 cos cos2 0 -sin2 1... 2 sin n cos On cos2 n - sin2 On where 0,, represents the inclination angle of the normal to crack n with respect to

95 the global x-axis. Again, the size of the transformation matrix N varies according to the number of cracks present at an integration point. The relationship between global stress increments and crack local stress increments is defined in a similar manner. First, we define the vector As' which assembles the stress increments of the individual cracks: As' = (As" At" As" t At.. as t At) (5.7) where As" and Att represent the normal and shear stress increments of crack number n, respectively. The transformation from global to local is accomplished using the same transformation matrix N listed above: As- = NTAa (5.8) At this stage, the stress-strain relations for the intact concrete and the smeared cracks need to be considered. For the intact concrete between the cracks, an incrementally linear relationship is assumed: Aa = DCAOec (5.9) where the matrix DCO expresses the constitutive properties of uncracked concrete. The SNAC implementation of this model assumes linear elastic isotropic material behavior for the intact concrete. Extensions to include nonlinear effects, such as plasticity, have proven successful [24,25], yet require additional considerations to insure compatibility with the smeared crack stress-strain relations outlined next. The vector containing the individual crack strain increments, Aecr, is related to the vector containing the individual crack stress increments, AsCr, through the following relations: Asr = DAecr (5.10)

96 D' 0 0 Dc = 0 Dc 0 (5.11) 0 0 DCr The elements of DC are 2 * 2 submatrices and the size of the matrix Dc depends on the number of cracks, n, acting concurrently at an integration point. By setting the off-diagonal submatrix terms to be zero, the relations expressed 'by Eqs. 5.10 and 5.11 assume no coupling effects occur between simultaneous cracks. This is a simplification of reality since material damage due to previous cracking will affect the crack initiation stress levels and energy absorption characteristics of subsequent cracking, especially for cracks at small angles to one another. Each diagonal submatrix element of DC expresses the crack interface stress-strain relations local to a particular crack. For crack n, the crack interface submatrix, DI', has been given the following form: DC 0 DCr Dcr 0 (5.12) 0 Gc where, referring to Fig. 5.4, we set: DC = DC for active cracks (i.e. cracks which are opening), D = DC for cracks which are closing/reopening, and Gc is the crack shear modulus. That the off-diagonal terms in Dc are zero is another simplification as in general there will be coupling between normal and shearing actions. It should be noted again

97 softening unloading/reloading nonlinear softening Figure 5.4: Crack normal softening relations [24] that the stress-strain relations normal to the crack generally need to be adjusted so that proper fracture energy is consumed during crack propagation in the numeric model. Failure to perform such modifications will generally lead to an analysis which is not objective with respect to mesh refinement. These concepts were discussed more fully in Chapter III and will be illustrated by example later in this section. The SNAC version of the above model relies on the simplicity of both a linear fracture function and a linear variation in shear stress with respect to shear strain. The solution of the linear set of equations may be subject to relatively large errors if the tangent stiffness matrix is not positive definite. To avoid numerical difficulties when using the standard Newton-Raphson technique the positive secant modulus, DC, is used to form the tangent stiffness matrix for both the active and arrested cracks. When updating the internal nodal force vector, F, the negative softening modulus, Do, is employed.

98 The stage is now set for deriving the constitutive relations for cracked concrete. Substituting Eqs. 5.4, 5.6, and 5.9 into Eq. 5.8 yields: AsC' = NTD [Ae - NAe'] (5.13) Combining with Eq. 5.10 results in a expression relating crack strains to total strain: Ae" = [DC' + NTDcoN]-lNTDcoA (5.14) Substituting this expression into: Ar = D~[Ae - NAeCr] (5.15) leads to the desired stress-strain relations for cracked concrete: Aor = [DC~ - DCON(Dc + NTDCoN)-lNTDCo]Ae (5.16) After cracking the material becomes incrementally linear anisotropic. Through Eq. 5.16, the effects of cracking may be interpreted as a reduction of the intact concrete stiffness DC. If the constitutive matrices DCO and Do are symmetric, the stress-strain relations for cracked concrete will also be symmetric. Furthermore, it has been shown that the above constitutive model obeys the principal of material frame indifference (i.e. objectivity with respect to reference frame orientation), provided the constitutive law for intact concrete obeys this same principal [24]. The SNAC implementation of the above cracking algorithm was used to numerically model the load-deflection response of a series of fracture specimens tested by Petersson at the Lund Institute of Technology, as reported in [81]. The test setup essentially consisted of a simply supported pre-notched fracture specimen under centerpoint loading. Direct tensile tests accompanied the fracture specimen study providing the material input parameters required for numeric analysis. The specimen geometry and material properties are presented in Fig. 5.5.

99 E = 30000 N/mm2 v = 0.2 ft = 3.33 N/mm2 Gf = 124 N/m = 0.001 [41 Ud I-r ---------- I x d x b = 2000 x 200 x 50mm a/d = 0.5 Figure 5.5: Fracture specimen geometry and material properties [81] To illustrate the importance of maintaining proper energy consumption during the fracture process, the analysis was performed using two different finite element models. Each model differs both in the degree of mesh refinement in the region surrounding the fracture zone, and in the slope given to the tensile softening curve (a coarser mesh will require a relatively steeper fracture function to insure proper energy consumption, as expressed by the material fracture energy Gf). Each model was "pre-notched" by providing an element-wide gap through the bottom half of the mesh. In all other aspects, the two finite element models were identical. Fig. 5.6 shows both the fine mesh and the coarse mesh models as they deform under loading. Loading was in the form of controlled displacements applied to the two central nodes on the upper face of the beam. The plots presented for each model correspond to centerline vertical displacements of 0.0, 0.4, 0.8, and 1.2 mm.,

100 respectively. The load-deflection response for the fine mesh model is compared with the experimental results in Fig. 5.7(a). As noted in [57], the use of four-node quadrilaterals with a four-point integration scheme leads to spurious cracking angles which tend to alternate in direction along the length of the fracture zone. Consequently, the structure will appear stiffer and more ductile than if the cracks were forming parallel to the plane of symmetry in a pure Mode I opening response. When employing a reduced one-point integration scheme for the elements comprising the fracture zone, cracking occurs along the plane of symmetry and the load-deflection response shows better agreement with the test results. The linear softening model does not provide the necessary degree of stress relief normal to the crack immediately after crack formation (see Fig. 5.4). Further improvements in the load-deflection response can be realized by utilizing a more general softening model (i.e. bilinear or nonlinear fracture functions) [57,81,82]. Fig. 5.7(b) compares the load-deflection response for the fine and coarse mesh models. Objective results are obtained provided the softening branch of the tensile relations has been adjusted so that the correct value of fracture energy is consumed during crack formation. On the other hand, analyses based purely on strength criteria (i.e. no softening) prove not only incapable of predicting the failure load of the specimen, but also show a marked dependence on mesh size. The capabilities of the non-orthogonal crack model are more dramatically illustrated in the following analysis of a shear critical reinforced concrete beam. This example is taken from a series of reinforced concrete beam tests recently performed at The University of Michigan [3]. The specimen under consideration can be considered simply supported and was subjected to displacement controlled loading applied

I I I I I I I I 1 56-6- U — y I I I I I I I, r,,,,,, t i, i 1 _7 1 1 t ---Lt-i1 1 I I f 1 ~ I I I jBi FINE MESH COARSE MESH Figure 5.6: Fracture specimen deformation during loading

102 1000.0 4-Point integration ~800. 0 Vi\ ~1-Point integrtion 800. 0 -\ -— I 600.0 z 0 2J 400.0 200.0 - D Experimental Scatter \ (Petersson) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Displacement (mm) (a) Numerical versus Experimental Results 1000. 0 ^,-~ Coarse mesh 800.ine mesh 600.0 - z 6 0 0 y ~ \1 j \ OBJECTIVE RESULTS 0 t0 / 9 400.0- 200.0- I \ / NON-OBJECTIVE RESULTS 0.0 — I I 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Displocement (mm) (b) Material Objectivity with Repsect to Mesh Size Figure 5.7: 'Fracture specimen load-deflection response

103 P/2 Section width = 7 in. (other dimensions in inches) shown:-.infl:: F.i.. 5.8, 7 iif whr stilymmeltry of tlthprlem. ha been exlitd. ^ —144~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 4;::) — -28 -------- ^~~~~~:::) TRUSS elements representing steel reinforcement along this level of the mesh at the one-third points along the span length (two-point loading). The structural shown in Fig. 5.8, where symmetry of the problem has been exploited. Flexural the test specimen. Material parameters used in the analysis are listed in Table 5.1. Fig. 5.9. The presence of cracking in the continuum elements is indicated by fine lines of "crack strains" in the continuum model. That is, the length of the crack lines indicates the relative magnitude of the crack strain at the integration point in ques

104 Concrete Properties fc = -5.8 ksi f = o.1lfcl EC = 4340 ksi (ACI Eq. 8.5.1) v = 0.2 D = -830 ksi* / = 0.2* Reinforcing Steel Area and Properties A, = 1.0 in2 E,= 29000 ksi f, = 75 ksi * based on estimates for fracture energy, Gf, and crack band width, h, according to Bazant and Oh [17] **,8 represents the so-called shear retention factor indicating the percentage of elastic shear capacity retained after cracking. Table 5.1: Material parameters for shear critical beam tion. The cracking trends witnessed in the test specimen have also been included in the figure, where the numbers alongside the cracks indicate crack advance in the test specimen with respect to load level (values in Kips). The numeric cracking patterns and the load levels at which they occur compare reasonably well with the actual test results. What is particularly significant about the numeric procedure is its ability to model the transition from a "flexurally dominant" cracking mode to a situation where diagonal shear cracks abruptly form and then begin to dominate the response. Traditional smeared cracking models generally restrict secondary cracking to occur orthogonal to the primary crack directions and, therefore, possibly preclude the type of "transitional behavior" witnessed in this example.

105,, i./ i - -- (load levels in Kips)............... 30,,..... 33 \I I i f \\\ EL pe. imnta Ret\\.9 ('\l~iload levelsl i P - 25.5 K P - 40.4 K P - 39.9 K P - 50.9 K Numeric Results Figure 5.9: Numeric cracking compared to experimental results

106 5.2.2 Plasticity based model w/abrupt tension cutoff A material model founded on the flow theory of plasticity was developed to simulate the nonlinear behavior of concrete in compression. Linear elastic isotropic material behavior is assumed up to an initial yield surface where upon a plasticity theory (associated flow rule) approach is used to construct the incremental stressstrain relations. The material hardens isotropically up to a limiting surface where it then begins to soften isotropically, eventually losing all load carrying capacity after violating a failure surface defined in terms of strain. Both the initial yield surface and the subsequent loading surfaces are chosen to approximate scaled versions of the biaxial test failure envelope reported in [56]. Background information regarding plasticity theory base constitutive formulations was presented in Chapter II. Hereafter, aspects of reinforced concrete modeling which require special considerations will be discussed, including the general forms of the yield function and hardening rules, cracking representation, and the crushing criterion. Mathematical descriptions of the limiting surface are generally formulated in terms of invariants of the stress (or strain) state. The following function is particularly suited to both concrete and rock-like materials: f(I7, J2, J3)= 0 (5.17) where i, is the first invariant of the Cauchy stress tensor, aij, and J1, J2 are the second and third invariants of the deviatoric stress tensor, sij, respectively (see Chapter II). A yield function in terms of only the first two invariants, I, and J2, has been utilized to model plane stress situations [27,71]. The SNAC plasticity model employs a yield function of the following form [71]: f(I,, J2) = [f(3J2) + ac]1/2 - ro = 0 (5.18)

107 where a and a3 are material parameters and (o is an effective stress measure extrapolated from uniaxial test results. Note that the von-Mises yield criterion is obtained by assuming a = 0 and 3 = 1.0. In terms of principal stresses, and setting 03 = 0 for plane stress, the relations in Eq. 5.18 become: [ea2 + a2 - 72] + a(l7 + a2)- 0a2 = 0 ((5.19) The material parameters have been chosen by matching biaxial test data at selected principal stress ratios. Using the results of Kupfer et al. [56], the material parameters become [71]: a = 0.335ao;, = 1.355 (5.20) The yield function expressed by Eqs. 5.18 through 5.20 appears as an elipse biaxial stress space, as portrayed in Fig. 5.10. The SNAC plasticity model allows the user to set the effective stress level which defines the initial yield surface. The initial yield surface indicated in Fig. 5.10 is for aolfo/f = 0.3. As the effective stress level increases during loading, the yield surface expands isotropically up until the limiting surface (aoo/f = 1.0) is met. Subsequent loading causes the effective stress level to drop and the yield surface will therefore contract. Plastic loading may continue until a failure surface (defined in terms of strain) is violated, whereupon the material is assumed to lose all stiffness and load carrying capacity. The rate of change in effective stress (and therefore the rate of expansion/contrac tion of the yield surface) during plastic loading is governed by the hardening relations. In the SNAC model, effective stress was assumed to be dependent on an effective strain measure defined by an accumulation of plastic strain increments (see Eq. 2.29). The particular form of the hardening rule is based on a parabolic idealization of the

108 -0. -1. -1.5 -1.0 -0.5 0.0 0.5 1.0 f(lo.J ([3J, +r = ~, Figure 5.10: Yield function in biaxial principal stress space concrete uniaxial response: = Eoe - E 2 (5.21) where Eo is the initial elasticity modulus, e is the total strain, and co is the total strain at peak stress fi. Substituting the elastic strain component, ce = al/Eo, into Eq. 5.21 results in the following expression relating stress to plastic strain (or effective stress to effective strain for the hardening rule): or=-Eoep+ 2Eo2ep for f< < a < f (5.22) where sp is the plastic strain component (or effective strain) and g is a parameter for marking the initial position of the yield surface. For plastic loading beyond the peak stress, fa, the hardening rule follows a linear softening branch according to a user specified slope. In order to model concrete cracking, the yield function expressed by Eqs. 5.18 through 5.20 is augmented with a simple tension cut-off surface, as illustrated in

109 Fig.'5.11 where Kupfer's test data have been included for comparison. Brittle tensile fracture is assumed to occur for any "uncracked" stress state violating the tension cut-off surface. Stresses normal to the cracking direction (i.e. stresses parallel to the principal tensile stress direction) are instantaneously released. The material response matrix then becomes orthotropic (see Eq. 3.2) with the direction of orthotropy aligned to the crack direction. The material stiffness normal to the crack direction is set to zero while the value along the crack direction either remains unaltered or is modified to reflect plasticity in the uniaxial stress state. A positive value for the shear modulus is retained to represent aggregate interlock, doweling action, and to avoid numeric difficulties associated with singularities which might arise if all integration points surrounding a given node crack in nearly the same direction. Cracks may close and reopen, however, once a crack forms its orientation remains fixed. Secondary cracking is possible in any orientation provided the primary crack has closed. In order to limit the concrete ultimate deformation capacity, a crushing criterion has been formulated in terms of strain. As suggested in [71], the yield criterion in terms of stress has been directly converted into a crushing criterion in terms of strain: f(l, J) = [3(3J2) + ac]1/2 - = (5.23) where I4 and J; are invariants of the strain state [32] and ~u is an ultimate effective strain measure. When any strain state violates the surface defined by Eq. 5.23, the material surrounding the sampling point in question is assumed to lose all its stiffness and load carrying capacity. The concrete model response in compression is shown in i Figure 5.12 for selected proportional loading paths. The initial yield surface was set at ao/fc = 0.3, the slope of the linear softening branch of the effective stress-strain curve was specified

110 Uimit Surfoe Tsion Cutoff_ Initij Yied 0 Kupfer. et.ol. 0..5 0.05 I'r/ I 1 10 I -1.5 -1.5 -1.0 -0.5 0.0 0.5 Figure 5.11: Yield function augmented with tension cut-off surface as Z = -150, and the ultimate effective strain value defining the crushing surface was set at e, = 0.0045. All other material parameters and the test data shown in the figures are from Kupfer et al. [56]. The numeric response agrees reasonably well with the test data, except for the case of a2/a1 = 0.5 where the numeric model overestimates material volume expansion at high stress levels. In the following example, the SNAC plasticity model is used to model the concrete behavior within a reinforced concrete shear wall specimen subjected to lateral loading. The particular specimen under consideration is one from a series of tests performed at the Portland Cement Association [69]. Material input parameters for the numeric model were taken from [69] and are listed in Table 5.2. The specimen geometry/undeformed finite element mesh are shown in Fig. 5.13. TRUSS elements were used to represent the flexural reinforcement within the wall boundaries. Both

111 1.5-1 0 1.0 0.0- ~,-. -," --- — ~ --- — -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 strain (tension positive) *10 -Pridcip stress ratio = 0.0 1.5 0.0 -5.0 -4.0 -3.0 -2.0 -1.0 0.0 strain (tension positive) *10-3 Principd stress ratio = 0.5 1.5 -gU 1.0 ' d E 0.5 0.0-,,, -6.0 -5.0 -4.0 -3.0 -2.0 -1.0 0.0 strain (tension Costive) *10~3 Rinb -tes rai - -. Prircipd stress ratio = tO Figure 5.12: Concrete model response in biaxial compression

Concrete Properties = -6.6 ksi Ec = 4080 ksi Reinforcing Steel Ratios and Properties pf = 0.0197 Ph = 0.0063 Pn = 0.0029 p, = 0.0135 E, = 27100 ksi fy = 65 ksi Values from Oesterle, et al. [69] pf is the ratio of main flexural reinforcement area to gross concrete area of boundary element. Ph is the ratio of horizontal shear reinforcement area to gross concrete area of a vertical section of the wall web. pn is the ratio of vertical shear reinforcement area to gross concrete area of a horizontal section of the wall web. p, is the ratio of effective volume of confining reinforcement to volume of core in accordance with Eq. A.4 of ACI 318-71. Table 5.2: Material input parameters for shear wall analysis the orthogonal steel reinforcing mesh in the web and the confining reinforcement in the boundary elements were modeled using the smeared representation option within the QUADC elements. The load-deflection response of the numeric model is compared with the experimental results in Figure 5.14. Cracking initiates along the tensile perimeter of the specimen within the first load step of the numeric solution (for a top lateral displacement, to,, of 0.1 inches). The load level initiating tensile yielding of the flexural reinforcement was 131 and 120 Kips for the numeric and experimental cases, respectively. The peak base shear level sustained by the numeric model (179 Kips) occured at top = 2.8 inches. Continued loading caused a gradual loss of specimen lateral strength up until Stop = 3.2 inches where the strength abruptly dropped off due to concrete crushing in the extreme compressive fibers of the lowest boundary element.

113 - 1 2 51-5 --- 12k — C CROSS-SECTI 1 2 CROSS-SECTION Po 8,I top 180 ELEVATION (dimensions in inches) Dimensions from Oesterle, et al. [69] Figure 5.13: Shear wall geometry and undeformed finite element mesh

114 200.0 -..,,// X 100.0 - Co. o 0.0-II M/I ~ Oesterle et al. (1979) Numeric response -100.0-. i I i I ' -1.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Top Laterd Displacement (in) Figure 5.14: Shear wall load-deflection response comparison Deformed plots of the numeric specimen, indicating both concrete tensile cracking and crushing, are presented in Figure 5.15 for the cases where Stop = 3.2, 3.3, and 4.9 inches. Forcing the numeric model to displace beyond Stop = 4.9 inches caused concrete crushing at numerous additional sampling points and the solution process was unable to regain equilibrium. A reproduced photograph of the test specimen (taken at Stop t 4.9 in.) is provided in Figure 5.16 for comparison. The numeric model does not show the strong presence of diagonal tension cracking witnessed in the test results. This is due, in part, to the numeric procedure's inability to model simultaneous non-orthogonal secondary cracking. Note that the test specimen shows some distress due to concrete crushing and spalling at locations in agreement with the numeric results.

8top = 3.2 in. ~e 1 - ''; -- 8top 3.3 in. L o^. ^.. I. <top 4.9 in. _: ~::~ ~ "', —,!y \. wall ci -- -. - ' ' \,| *wall ci C01;racking,rushing Figure 5.15: Shear wall deformation at selected load levels

116 Figure 5.16: Shear wall crack pattern during testing (Stop: 4.9 in.) [69] 5.3 TRUSS - Two Node Truss Element The TRUSS element supported by SNAC is a version of the truss-type element presented in DLEARN [50] which has been modified to allow for material nonlinearity. Material stress-strain (or element force-deformation) relations are idealized as bilinear elastic strain-hardening plastic and follow either a kinematic or an isotropic hardening rule, as shown in Figure 5.17. The TRUSS element is included as a means of representing the reinforcing bars in a discrete manner, as opposed to smearing the effects of the reinforcing steel over the element volume. Both the discrete and the smeared representations could be formulated to provide similar results in a global sense. However, the discrete modeling approach has been favored in anticipation of extending the analysis to include the effects of bond slip.

117 100.0 50 0. -50. -100.0 sttroin ' 6.0 8.0 o) Kinematic hording rule 3 100.0 stress (ks) -50 -100. b) Isotropic hardening rule Figure 5.17: TRUSS element bilinear material idealization

118 5.4 BEAM - Two Node Beam-Column Element 5.4.1 Two surface plasticity model in 1-D action space The BEAM (beam-column) element supported by SNAC is basically a reduced version of the 3-D beam-column element reported in [31,75]. More specifically, the 3-D beam-column element, which allows for four-dimensional action interaction (i.e. interaction between biaxial bending moments, axial load, and torque), has been simplified to a 2-D beam-column element possessing no action interaction capability. The basic theories and computational procedures underlying the more general 3-D model have been utilized in the formulation of the reduced model. However, certain difficulties which must be addressed when considering action interaction (e.g. stress state corrections, overlapping of the yield surfaces) do not arise in the simplified model. The essential features of the BEAM element are now briefly outlined. The interested reader is referred to [31,75] for specific details regarding both the element formulation and concepts related to multidimensional interaction surfaces. 1. The BEAM element is a discrete line element which may be used to represent structural members (e.g. beams, beam-columns, slender walls) which are arbitrarily orientated in two dimensions. The BEAM element possesses both translational and rotational degrees of freedom identified at each of its end nodes. 2. A "section" type model is used to represent element inelastic behavior. That is, inelastic behavior is defined for the cross-section as a whole through a set of predefined force-deformation relations. This is in contrast to a "fiber" type modeling approach, where a member cross-section is discretized into a number

119 of smaller areas called fibers. The material comprising each individual fiber is assumed to undergo uniform axial deformation and behave according to prescribed uniaxial stress-strain rules. Cross-section stiffness is then obtained through proper summation of the individual fiber tangent stiffnesses. Section type models are more abstract and phenomenologically dependent relative to fiber type models which are based on more deterministic quantities (e.g. relatively detailed analytical descriptions of both the cross-section geometry and the material properties). However, the section type model is attractive from a computational cost standpoint. 3. Inelastic behavior is approximated by assuming the element consists of a linear elastic beam element with a nonlinear zero-length hinge at each end, as shown in Figure 5.18. The elastic portion of the beam is defined by an axial stiffness, a flexural stiffness, and an effective shear rigidity. Hinge nonlinearity is assumed to be dependent on the bending moment acting at the corresponding end node. 4. For monotonic loading, the hinges provide a rigid-plastic-strain-hardening response. Strain-hardening stiffnesses are specified via moment-rotation relationships for the hinges. Upon load reversal, an elastic component is introduced into the hinge moment-rotation relations to simulate element stiffness degradation. Stiffness degradation is dependent upon previous plastic deformations in the hinge and is controlled through user specified coefficients [31]. Hinge response may then be referred to as elastic-strain-hardening-plastic. 5. Each hinge is comprised of two subhinges, which can be characterized as "cracking" and "yielding" subhinges. Bilinear action-deformation (moment-rotation)

121 subhingit? 600.0 400.0 200.0 0.0 -200.0 -400.0 600.0 400.0 200.0 0.0 -200. 0 -400.0 -0.05 0.00 0.05 0.10 Subhinge Rotations (radions) Subhinge Moment-Rotation Diagrams Hinge Rotation (radians) Hinge Moment-Rotation Diagram Figure 5.19: Hinge moment-rotation response during cyclic loading

122 translate in action space until the yield surface corresponding to the yielding subhinge is met. Continued plastic loading will then cause both yield surfaces to translate together in action space. 5.4.2 Hinge plastic stiffness determination All inelastic behavior is concentrated within the zero-length hinges located at the element end nodes. Each hinge consists of a two subhinges connected in series. The behavior of each subhinge follows predefined bilinear action-deformation (momentrotation) relations resulting in a trilinear action-deformation relationship for the complete hinge. While this type of model explicitly accounts for cracking, yielding, and stiffness degradation within the member, the appropriate hinge plastic stiffness values are difficult to establish, especially for general loading histories. This shortcoming is inherent in any lumped plasticity model where stiffness, strength, and hardening characteristics derived at the cross-section level are used to represent inelastic behavior at the element level. Since the distribution of inelastic activity in the member depends on the moment variation over the member length, the element stiffness characteristics should also show dependence on this same moment variation. Certain assumptions must therefore be made regarding the type of moment variation over the beam length. In some instances, basic moment profiles may adequately describe the actual loading conditions. For example, a uniform moment may be used in situations where the moment is not expected to vary appreciably over the element length. A linearly varying moment with equal and opposite end values may be used to model the moment profile acting on coupling beams between shear walls. The assumed moment

123 profile, in conjunction with the moment-curvature properties of the cross-section, can be used to determine the hinge plastic stiffness values. Cross-section moment-curvature relations were determined in this study using standard techniques. That is, the cross-section was first discretized over its depth into a number of small areas or "fibers". The strain experienced by a given fiber depended linearly on the fiber distance from the neutral axis (i.e. plane sections remain plane). Fiber response to deformation was defined by prescribed uniaxial stress-strain relations according to the type of material comprising the fiber. Figure 5.20 shows the material models used in this study to describe the concrete and reinforcing steel behaviors. For a given curvature, the sectional moment was found by summing the moment contributions of the individual fibers. Moment-curvature relations were idealized to be trilinear with breaks in slope at the point of cracking and the point of yielding. The point of ultimate behavior was determined by a concrete strain of 0.004 in the extreme compressive fiber of the section. This ultimate point was then used to define the post-yield slope of the trilinear moment-curvature relation. For the case where uniform moment acts over a member, end node equivalent rotations at cracking, yield, and ultimate (as shown in Figure 5.21) were found by integrating the corresponding curvatures over the member length. The subhinge plastic stiffness values were then retrieved by noting that the elastic beam and the two subhinges act as springs connected in series. Referring to Figure 5.21: Kp KK2 (5.24) KI - (2 Kp2= K2- K (5.25) 2 - K3

124 fC c so E a) Confined and unconfined concrete a lu y E Esh ~u E b) Reinforcing Steel Figure 5.20: Uniaxial material properties for sectional analysis

125 M My-.. 1- K3 K2 Mc c r ey Ou 0 Figure 5.21: Moment vs. rotation at element end node [31] where Kpx and Kp2 are the desired "cracking" and "yielding" subhinge plastic stiffness values, respectively. For the case of linearly varying moment with equal and opposite end values, the determination of the subhinge plastic stiffness values follows a slightly different path. Since an inflection point occurs at the member midspan, each half length of the member can be regarded as an equivalent cantilever subjected to a tip load. Using the trilinear moment-curvature relation for the section, together with the second moment area theorem from elastic beam analysis, the relation between the cantilever tip load, P, and the tip displacement, 6, was found. Referring to Figure 5.22, the subhinge plastic stiffness values were extracted in the same manner as before: KK2 L (5.26) K1 - K2 2 Kp K2K3 L (5.27) K2 - K3 2

126 PT PU Py 8cr 5y 8u 6 Figure 5.22: Tip load vs. tip deflection in equivalent cantilever [31] The effect of inelastic shear deformation was roughly accounted for by increasing the cantilever tip deflection by an appropriate factor. This factor was assumed to be the ratio of total tip deflection (bending plus shearing effects) to tip deflection considering only bending effects, where both terms in the ratio were determined for an elastic cantilever with rectangular section. Such a procedure probably underestimates the loss of stiffness due to inelastic shear deformation. However, for members possessing a small shear span to depth ratio significant reductions in subhinge plastic stiffness values are still realized. Note that this technique is not applicable to the uniform moment model since the relationship between shear and moment during the loading history is not uniquely defined as it is for the cantilever model.

CHAPTER VI CHILEAN BUILDING ANALYSIS 6.1 Background On March 3, 1985 a strong motion earthquake, magnitude M, = 7.8, occurred just off the coast of central Chile. Due to the history of seismicity in that region an extensive network of strong motion instrumentation was in place to record the event. Recorded acceleration histories indicate the strong motion occurred throughout the epicenteral region with several stations reporting durations of strong motion (> O.lg) on the order of 60 seconds. Four stations reported peak horizontal accelerations greater than or equal to 0.6g [94]. The earthquake caused considerable damage and casualties in the central region of Chile, particularly due to the failure of low rise adobe and unreinforced masonry dwellings [90]. In contrast, most of the modern engineering structures performed reasonably well considering the energy content of the input ground motion. One remarkable feature was the exceptional performance of most of the more than 400 medium to medium-high rise buildings in the coastal city of Vinia del Mar. In response to these extraordinary circumstances, a joint effort between the the University of California at Berkeley, the University of Illinois (Urbana), and The University of Michigan was formed to investigate the buildings of the Vina del Mar region. As 127

128 part of these ongoing efforts, funded by the United States National Science Foundation, several works have been published from which the following information was adopted [80,94]: * All buildings in Vinia del Mar with five or more stories were constructed of reinforced concrete. Shear walls were predominately the most common structural system for carrying gravity and lateral loads. In most cases, these walls tend to follow the architectural layout of the floor rather than being strictly limited to locations along the building perimeter or central core region. * Very few buildings suffered major damage. In most of the cases where significant structural or nonstructural damage was incurred, problems linked to either local soil-site amplication or unrepaired damage from prior earthquakes were cited as possible causes. Most buildings had only minor or superficial damage. * Ninety percent of the total 415 buildings considered in the survey were in the 5 to 10 story range. Less than three percent of these building were greater than 20 stories. * The ratio of wall to plan area for these buildings was much greater than that of customary designs in the United States. The wall to plan area index for the Chilean buildings ranged from 3 to 8 percent (roughly independent of building height) with an average value near 6 percent. By assuming the seismic base shear to be 10 percent of the building weight, it can be shown that most of the Chilean buildings were designed to produce an average shear stress of less than 5 kg/cm2 (70 psi) in the shear walls.

129 * The Chilean Code for reinforced concrete does not contain provisions for seismic design, and in particular no special detailing requirements to insure ductile behavior as is customary for buildings in seismic regions of the United States. The Chilean engineers do look to the ACI Code (Appendix A) for detailing frame type structures, but seldom use it for walls. No limits are placed on the amount of flexural reinforcement along the wall boundaries and no specifications exist regarding proper confinement of the these boundary elements. * For the most part, no clear trends were established when trying to link observed building damage to such variables as date of construction, structural type, building height, wall area to plan area index. It was concluded that more detailed studies would be required to make a better assessment of the performance of these buildings. One of the buildings included in the above survey, the Villa Real building, will be analyzed as outlined in the following sections. 6.2 Building Description The Villa Real complex consists of a ten-story residential building, a subterranean parking structure that surrounds the building on three sides, and a small commercial center adjacent to the fourth side of the building. Villa Real is located at 5 Norte/4 Poniente close to the downtown area of Vinia del Mar. The building design was finished in 1981 and the construction was completed in 1983. The main structure is a slab-girder-shear wall system, the shear walls being irregularly spaced and varying dramatically in plan from the first to second story (Figure 6.1). Stories two through eight are similar in plan, while stories nine and ten again differ. Figures 6.2 and 6.3 show building cross-section sets and indicate the

130 floor/story numbering sequences in each building direction. Additional information related to the main building dimensions, including the specified material properties, is listed in Table 6.1. The building shear walls range in thickness from 20 to 30 centimeters near the base of the structure. However, the majority of the walls in the upper stories are 20 cm thick. Each wall contains two curtains of uniformly spaced orthogonal steel reinforcing mesh with reinforcing ratios ranging between 0.0017 to 0.0040. The total building complex plan area is approximately 800 square meters at the foundation level. The main building, the parking structure, and the commercial center are connected at the foundation level by a network of foundation beams. Continuous wall footings on a 70 cm mat-type foundation support the main building. As reported in [94], structural damage due to the earthquake was minimal and was concentrated near the base of the structure. Diagonal cracking was observed in the structural walls of the first, second, and third stories with most cracks measuring less than 0.2 mm wide. 6.3 Modeling of the Structure The structural idealizations presented in this section reflect a macro-model approach to the building analysis. The building global response during loading is of primary interest and the precise details of the localized events contributing to the overall reponse are generally neglected. Loading will be applied to the structure in the E-W direction with no loads acting simultaneously in the N-S orthogonal direction (further details regarding the loading are presented in Section 6.3). It is assumed that the shear walls provide relatively little resistance to out-of-plane deformation. Therefore, only the walls in the direction

131 ~ L=:3 xll I iiI_ 1 I I I I I I - I I x I 14 4 m m I i I. a I I —.-j I - - 18.6 m Second Story Plan E N S w - ndit shr wal ovr story height -= Ind ohl am slEnd olv I L I - I I U 1 m m[ (bi - mb ( 1'm - --- - —, --- - 1A& - ^-16-.6 -T.4 m First Story Plan - Wnd1m uhw -all ow smcy neNpH ti= k*mi _ a t fim Ii Figure 6.1: Villa Real building plan views

IR 13110 l 10 9 8 8 7 8 l7 6 5 5 4 3 2 2 0 SECTiN 3-3 and 7-7 SECTION 4-4 SECTION 3-3 and 7-7 SECTION 4-4 SECTION 2-2 SECTION 5-5 SECTION 6-6 SECTION 9-9 Figure 6.2: Villa Real cross-sections (E-W direction)

10 6 4 2 11 1 SECTION D-D 10 9 8 7 6 5 4 3 2 1 0 no ni SECTION F-F SECTION F-F SECTION C-C SECTION E-E SECTION G-G Figure 6.3: Villa Real cross-sections (N-S direction)

134 Total building height = 28.3 m. (92.5') - above grade, incl. elevator housing = -3.8 m. (12.5') - below grade to foundation base Bldg. aspect ratio* = 2 to 1 East-West = 1.6 to 1 North-South Number of stories = 10 Typical story height = 2.55 m. (8'-4") Floor slab thickness = 12 to 16 cm. Concrete strength** = 300 kgf/cm2 (: 4.3 ksi) Reinf. steel strength = 4200 kgf/cm2 (x 60 ksi) - in shear walls =2800 kgf/cm2 (: 40 ksi) - in slabs Wall area/floor area = 0.072 - first floor = 0.059 - second floor * Building height to plan dimension at base. **Specified 28-day cubic resistance (R2s). Table 6.1: Villa Real dimensions and specified material properties considered will be directly modeled using finite elements. In some instances, walls orthogonal to the motion considered will act as flanges attached to the in-plane walls. These cases will be indirectly accounted for in the finite element model as outlined in Section 6.3.5. The QUAD4 element will be used to represent the concrete shear wall continuum in a state of plane stress. BEAM elements will be used to represent both the coupling beams and certain walls which have high aspect ratios and therefore lend themselves better to a discrete line element modeling. TRUSS elements are used to model the layers of flexural reinforcing steel near the shear wall boundaries. Details regarding these elements were presented in Chapter V. The undeformed finite element mesh is shown in Figure 6.4, where the QUAD4 elements comprising the shear walls have been shaded to help distinguish the BEAM

135 elements representing the coupling beams and slender walls. The bold lines in Figure 6.5 indicate the locations of the TRUSS elements representing the flexural reinforcing steel near the wall boundaries. The TRUSS element end nodes (not discernable in Figure 6.5) match the corresponding QUAD4 element nodes in Figure 6.4. In addition to the numerous assumptions inherent in the material models and input values, the following assumptions and modeling techniques were utilized. 6.3.1 Upper story idealization The upper five stories of the building structure were assumed to remain undamaged (linear elastic) during the loading history. This decision is reflected in the finite element discretization where a coarser mesh of QUAD4 elements has been used to represent the top five stories of the building. Constraint equations have been introduced to insure displacement compatability at the transition region. While inelastic deformations in the upper stories would contribute to the overall drift of the structure, generally most of the nonlinear actions contributing to the overall failure state of the structure should be confined to the first few stories of the structure. The results presented in Section 6.3 tend to bear these assumptions out. It may also be argued, in the spirit of St. Venant's principal, that the rough modeling of the upper stories will not appreciably affect the analysis in the lowermost stories. Besides the great amount of computational savings already realized by the above idealization, additional savings could also be obtained by statically condensing out all 'unnecessary' degrees of freedom within the linear elastic upper story region.

I-A 0> 4% * \ \ ~t: Figure 6.4: Building finite element model (E-W direction)

Figure 6.5: TRUSS element modeling of exural reinforcement Figure 6.5: TRUSS element modeling of flexural reinforcement

138 6.3.2 Rigid diaphram assumption The building floor slabs were assumed to act as rigid diaphrams with respect to lateral translation, but provide no resistance with respect to vertical translation and rotation. The rigid diaphram modeling was effectively carried out using constraint equations to insure each nodal point at a given floor level translated equally. That is, all lateral displacements at a given floor level were 'slaved' to a single 'master' node. This assumption also facilitated lateral load application to the structure. Whether or not the rigid diaphram assumption is an acceptable approximation to the slab behavior depends largely on the nature of the structure being analyzed. For Villa Real, similarity in plan layout from Floor 2 through Floor 9 implies little change in the general resisting behavior of each wall element (relative to the surrounding structure) over this height of the building. In addition, the walls are spaced and proportioned in such a manner that lateral load resisting stiffness seems to be well distributed. It should therefore be acceptable, as an approximation, to assume rigid floor diaphram action with respect to lateral displacements over this height. However, due to the introduction of Walls 5-5 and 7-7 in the lower stories the validity of the rigid diaphram assumption at these levels becomes somewhat more tenuous. The walls which continue over the building height will tend to carry large amounts of shear down towards the base. It is unreasonable for a wall which abruptly begins at a certain level (e.g. Wall 5-5 at Level 1 or Wall 7-7 at Level 0, as shown in Fig. 6.4) to immediately pick up its entire 'share' of the total building shear. That is, a realistic modeling of the floor slab system must somehow account for the slab's limited ability to transfer loading. A general formulation allowing for slab flexibility is desirable, yet would present considerable difficulty to implement. A relatively straitforward means to account

139 for finite slab flexibility would be to simply adjust the coefficients in the constraint equations linking 'slave' to 'master' node at the level in question. This simple flexible slab model is presented in Section 6.6.2. 6.3.3 Equivalent column modeling of shear walls Wall 3-3 and the upper portions of Walls 5-5 and 7-7 were modeled using BEAM elements. A line element modeling of these members is more appropriate due to their slender nature. A certain amount of engineering judgement was exercised when idealizing Walls 3-3 and 7-7. The coupling beams and outer portions of these walls were neglected as they are relatively flexible compared to the central core portion which has sizable flanges (orthogonal walls) framing in. Effective flange properties and widths were determined as outlined in Section 6.3.5. The two-surface (cracking and yielding) plasticity model presented in Section 5.4 was employed allowing for a trilinear moment-rotation description of the elements nonlinear behavior. Nonlinear material input parameters were determined using the section analysis outlined in Section 5.4.2 under the assumption of uniform moment acting over the length of each element and constant axial load acting for the duration of the loading history. The axial load was set equal to the gravity load attributed to the element in question. 6.3.4 Coupling beams The coupling beams within Walls 2-2 and 9-9 are represented by BEAM elements with trilinear moment-rotation capabilities. The nonlinear material properties were determined as outlined in Section 5.4.2 under the assumption of equal and opposite end moment bending and no axial load (the rigid diaphram assumption precludes

140 axial loading in the coupling beams). Rotational stiffness was effectively supplied to the coupling beam end nodes via displacement constraints. Floor slabs were assumed to contribute to the strength of the coupling beams allowing for a total effective flange width of L/4 where the slab occurs on both sides of the coupling beam and L is the length of the coupling beam clear span. If the coupling beam were subjected to only single curvature bending over its length, the effective flange width should be greater than L/4 due to continuity of the slab and beam reinforcement well beyond the clear span length. However, double curvature bending presents conditions nearly analogous to those considered by the ACI Building Code when specifying an effective flange width of L/4 for simply supported beam bending in single curvature. Likewise, where the slab is found only on one side of the coupling beam, the effective flange width was set at L/12. The entire effective flange width was assumed to act when bending caused compression in the top fibers of the coupling beam, while only the coupling beam width contributes to strength when bending caused the top fibers to be in tension. For the latter case, slab reinforcing steel within the effective flange would normally contribute the strength of the section, however, for the effective flange lengths arising in this analysis these contributions were found to be negligible. 6.3.5 Orthogonal walls The effects of orthogonal walls framing into the shear walls in the lower stories of the structure were accounted for in the finite element modeling. The finite element discretization was chosen to reflect the presence of orthogonal walls, as can be seen in Figure 6.4. The elements which appear somewhat narrower over the story height indicate the locations where orthogonal walls frame in.

141 In general, these orthogonal walls act as flanges which may greatly increase the strength and stiffness of the walls in the direction being analyzed. Since we are performing an analysis were tensile cracking and other nonlinearities are allowed to occur, the width and material properties of the flange will nessecarily differ depending on the direction of loading (compression or tension in the flange). Actually, a realistic representation of the progress of inelastic actions in these flanges during loading requires a three-dimensional analysis. However, the analysis is limited here by the framework of a two-dimensional model. Without relying on some sort of unconventional technique, one must therefore anticipate the near ultimate state of the three-dimensional phenomena and use this information to determine the input parameters for the two-dimensional case. However unavoidable, this modeling is still objectionable since it is unable to account for the path dependent nature of reinforced concrete during loading. For flanges acting in compression, the element thickness (or effective flange width) is determined by assuming the compression field created in the flange due to actions in the wall spreads out at a 1:1 slope away from the wall. The size of the effective flange increases in this manner up until the limits set by the ACI Code are reached [1]: b, + 6hf be = min b (bw + S)/2 for slab on one side of the member (6.1) b, + 1/2 b, + 16hf b= min S for slab on two sides of the member (6.2) 1/4

142 where: be: effective flange width, b,o: web width, hf: flange thickness, S: spacing between webs, I: member span length. The reinforcing steel contained in the effective flange tip regions is neglected in compression. For flanges acting in tension, the element thickness (effective flange width) is not increased since we anticipate tensile cracking of the concrete in this region. On the other hand, the reinforcing steel in the tension flange may significantly contribute to the ultimate strength of the wall section. By assuming the tension field created in the flange due to actions of the wall spreads out at a 2:1 slope away from the wall near ultimate [74], the amount of flange reinforcing steel which is mobilized near ultimate conditions can be approximated. This area of steel is then lumped together with the flexural reinforcing steel in the wall as an input parameter. In this study, the size of the tensile region defined by the 2:1 slope is limited only by the obvious restrictions caused by adjacent tension zones due to actions in other walls and the flange length itself. Both these limitations and the actual height of the walls kept the tension region size within reason and it was therefore unnecessary to define any other practical limitations. It should be noted that any attempt to define effective widths of the compression and tension flanges is somewhat clouded by the possibility of large components of earthquake motion orthogonal to the direction being considered. Due to actions

143 caused by out-of-plane motion, a portion of flange area thought to be in compression may actually be in tension, or visa versa. However, this is not so critical for the analysis presented here since the orthogonal walls which frame into the lower few stories are not tall and would therefore not be expected to undergo large overturning moments. 6.3.6 Fixed base assumption The structure is assumed to be fixed at its base to an infinitely rigid environment. The flexibility of the surrounding soil and any effects related to soil structure interaction will therefore not be accounted for in this analysis. Neglecting these concerns is not so objectionable for several reasons: 1. Due to the high degree of strength and redundancy of the structural walls in the subterranean level of the structure, its not unreasonable to assume the base would act against the surrounding soil in nearly a rigid body fashion. The structural walls above the base would therefore not see any high degree of differential movement with respect to one another which might occur given a flexible basement level. 2. The aspect ratio (height/plan length) of the building in the direction of motion considered is approximately 2800/1400 = 2.0. Considering a nearly rigid subterranean level, the stresses on the surrounding soil due to the overturning moment of the building during loading should be well distributed and relatively small in magnitude. 3. For the type of loading considered (quasistatic monotonic), again assuming a nearly rigid subterranean level, the main consequence of providing a flexible soil

144 structure would be an increase in building drift and therefore some contribution to the secondary P-6 effect. Problems related to soil-structure resonance could only be detected in some sort of dynamic analysis. 6.3.7 Reinforcing mesh representation The two curtains of steel reinforcing mesh contained in each shear wall were accounted for in the following manner. For walls modeled using continuum elements, the effects of the reinforcing mesh in each orthogonal direction were assumed to be smeared' over the element volume as described in Section 5.2. Nonlinear behavior was assumed to be elastic perfectly plastic. For walls modeled using the line element approach, the vertical component of the reinforcing mesh was included when determining member input properties using the section analysis outlined in Section 5.4. In the same manner, reinforcing steel details called for in the building plans were roughly accounted for in the continuum model by smearing the steel quantities out over the region in question. Locations where such special detailing was considered include: the regions just above the openings in the first story of Wall 2-2 and around the reentrant corners on the right side of Walls 4-4, 6-6, and 9-9. 6.3.8 Nonstructural elements While the presence of nonstructural elements was considered in the building gravity load calculations, no attempt was made to model their contribution to the structure's overall stiffness and strength reserves. Studies have indicated that nonstructural components may greatly influence structural load resistance behavior. However, the degree in which these nonstructural components participate in the overall building response depends largely on the structural type and the character of the

145 Kitchen Total Level Slab Beam Terrace & Bath Walls Stairs Others* (T) 0 108.6 24.0 ---- 1.6 165.3 4.7 4.7 308.9 1 112.6 28.9 2.8 2.4 95.6 4.7 10.8 257.9 2 106.2 24.2 2.8 2.4 91.5 4.7 12.2 244.0 3 104.8 22.2 2.8 2.4 88.7 4.7 12.2 237.8 4 104.8 21.7 2.8 2.4 85.9 4.7 12.2 234.5 5 104.8 21.7 2.8 2.4 85.9 4.7 12.2 234.5 6 104.8 21.7 2.8 2.4 85.0 4.7 12.2 233.5 7 104.3 21.7 2.8 2.4 83.7 4.7 12.2 231.8 8 92.9 26.9 2.8 2.4 75.7 4.7 12.2 217.5 9 77.7 23.2 3.9 1.7 60.7 4.7 7.5 179.4 10 58.4 34.3 19.0** -- 26.7 2.4 1.4 142.2 2522 * -- Partitions and windows. * -- Corresponds to weight of roofing and elevator Table 6.2: Villa Real building gravity load summary nonstructural components. The structural system of the Villa Real building is stiff to begin with and the effects of the various nonstructural elements (e.g. limited brick and wood partitions, an elevator casing, stairways) should therefore be relatively less influencial on the building response. 6.4 Modeling of the Loading 6.4.1 Gravity loading The gravity loads assigned to each floor level were based on a set of calculations performed by students at Santa Maria University in Valparaiso, Chile and are presented in Table 6.2. The wall, stair, and partition loadings lumped to a particular level were determined by transfering the corresponding weights lying within one-half the story heights above and below the same level.

146 Gravity loads were assumed to act uniformly over each respective floor level. These gravity loads varied in magnitude from approximately 1300 kg/m2 at the first floor level to roughly 900 kg/m2 for Levels 2 through 8. This range of values agrees reasonably well with the 1000 kg/m2 uniform gravity load deemed typical for the other Chilean buildings investigated in the joint effort. Uniform gravity loads were distributed amongst the supporting walls assuming that a wall or beam attracts the load lying within a certain "tributary area". In most cases, tributary areas were defined by boundaries bisecting the 90~ angles at the slab span corners. The manner in which gravity loads acted in the finite element analysis depended on the type of idealization representing the wall. For walls modeled using the continuum element approach, the gravity load magnitude was further divided into equivalent nodal loads. Nodal load values were determined by dividing the total load lumped to the wall by the wall length, and then multiplying this value by the tributary length assigned to each node along the wall length. This conversion will generally not represent a true uniform surface load, but it is close enough for the purposes of this analysis. As discussed in Section 6.3.5, orthogonal walls frame into the in-plane walls at various locations in the lower stories. It was therefore necessary to increase the gravity load on the 'thickened' elements at these locations to avoid any inconsistency. The uniform load per length of orthogonal wall was multiplied by the length of the 'flange' portion and the result was lumped to the corresponding nodes. For slender walls modeled using the discrete line element approach, the effect of gravity load was incorporated at the section analysis stage (see Section 5.4). That is, the gravity load on these elements was not directly applied to the finite element model. but rather it was used in the determination of the BEAM element material

147 properties. Note that the finite element idealization of the structure specifies only horizontal link-type coupling between the several walls. Therefore, these BEAM elements will not experience any axial load variance during the loading history. The gravity load for a BEAM element in a given story was determined by summing the gravity load contributions at every floor level above the story in question. 6.4.2 Lateral loading In addition to gravity loading, the building was subjected to in-plane lateral loadings which were monotonically increased (proportionally) until structural failure criteria were violated. Failure was defined to occur when either the building structure begins to lose its load carrying capacity or the lateral displacements become large so that secondary effects become significant. The lateral load variance will be assumed to take place in a "quasi-static" manner (i.e. as if the loads are applied very slowly relative to the natural periods of the structure). While inertial effects of the structural masses are not explicitly accounted for in the finite element formulation, they may be implied by assuming a lateral loading distribution which approximates the fundamental mode response of a structure during truly dynamic motion. Studies have shown the fundamental mode response to be a reasonable approximation to the building response, except for slender buildings where higher mode responses become relatively more significant. The lateral load distribution over the building height was therefore based on a triangular acceleration profile as an attempt to approximate the building fundamental mode response. Additional discussions related to the lateral load distribution profile are presented in Section 6.4.3 and Section 6.6.5.

148 Another consequence of the quasi-static load variance assumption is that strain rate effects on the material behavior not explicitly accounted for in the finite elememt formulation. It has been reported that strain rate does affect the material properties of both plain concrete [22] and reinforced concrete [59]; both strength and initial stiffness increasing with increasing strain rate. However, these strain rate effects are greatest only for loading within the elastic range or the initial excursion into the inelastic range. Subsequent excursions into the inelastic region show that strain rate has little influence on stiffness or strength of the material. Severe earthquake motion may be typified by relatively few strong load cycles producing stress histories well into the inelastic range in reinforced concrete. Previous earthquakes or even foreshocks to the strong input motions will also produce inelastic activity in the reinforced concrete. In light of such considerations it has been suggested that strain rate effects on material behavior may be neglected for most nonlinear analyses of reinforced concrete structures [88]. It should be noted, however, that for blast or impact type loading situations strain rates may be at least two orders of magnitude higher than those due to earthquake excitations and the material may have seen relatively little previous inelastic activity (e. g. nuclear power plants, dams). For such situations strain rate effects might considerably influence the resulting structural behavior. 6.4.3 Uniform Building Code lateral load distribution The Uniform Building Code implies a heightwise distribution of lateral forces based on a triangular horizontal acceleration profile (horizontal accelerations linearly proportional to the elevation above ground). An additional lateral force is added at

149 the roof level to account for higher mode contributions to the building response in longer period structures (T > 0.7 sec.). The Applied Technology Council-3 Recommendations use a similar approach. For short period structures (T < 0.5 sec.) a triangular acceleration profile is again used. For long period structures (T > 2.5 sec.) the heightwise lateral force distribution is based on an acceleration profile proportional to the square of the elevation above the base (a parabolic acceleration profile). Structures with intermediate period length are subjected to lateral load profiles based on acceleration profiles varying according to an intermediate power (1 to 2) of the height above the base. The exponential variance of the acceleration profile is an attempt to account for the higher modes contributions to the overall building response. Preliminary studies have estimated the Villa Real building fundamental period to be slightly greater than T = 0.4 seconds for motion directed along either building principal axis [23]. Other studies have indicated the Chilean buildings are generally quite stiff and N/20 may be a reasonable approximation to the fundamental period, where N is the number of building stories. Either approach indicates the Villa Real fundamental period to be within the short period range implying the 'code' use of a triangular, or nearly triangular, acceleration profile to establish the lateral load distribution for use in the static analysis. One goal of this study was to compare the total base shear on the building model at failure with that specified by the Uniform Building Code. The UBC specified base shear may be determined using UBC Formula (12-1), as given below. Assuming the building fundamental period to be 0.5 seconds, the calculated base shear for Villa Real is 380 metric tons (which is approximately 15% of the building dead load plus partition weight). This base shear is distributed over the building height in

150 accordance with UBC provisions and the resulting lateral load profile is shown in Figure 6.6. V = ZIKCSW (6.3) where: V: total lateral force (or base shear) on the building, Z: Zone Coefficient (= 1.0 for Zone 4), I: Occupancy Importance Factor (= 1.0 for Villa Real), K: Horizontal Force Factor (= 1.33 for Villa Real), C: numeric coefficient related to building period, S: Soil Profile Coefficient (= 1.2 for Profile Type S2), W: total dead load plus partition weight for building. 6.5 Analysis Results The building model was analyzed numerous times using various load increment sizes and convergence tolerances. Results from the analyses proved to be relatively insensitive to reasonable adjustments in these parameters about the selected values. The load vector magnitude and load-time functions were scaled in a manner which allowed the lateral load profile to vary linearly over time and equal the UBC distribution at time t = 1.0. Both "solution time" and "normalized base shear" (i.e. V/VUBc) may therefore be referred to interchangably. By successively decreasing the load step size as the failure state was approached, the failure load of the numerical model was established. Hereafter, reference to structural failure will be understood as failure of the numerical model.

151 30.0 - 25.0- \ /.20.0- o / f 15.0- 5.0 0.0.0 0.0 10.0 20.0 30.0 40.0 50.0 60.0 Laterd Force (metric tons) Figure 6.6: Uniform Building Code lateral load distribution ILoof level lateral displacement versus normalized base shear is plotted for both directions of loading in Figure 6.7. For loading acting towards the East (positive displacement), structural failure occurred at time t = 2.79 (i.e. at a base shear of V = 2.79VUBC = 1060 metric tons which equals roughly 42% of the building gravity load in magnitude). The building drift at failure was approximately 1.5%. For loading in the opposite direction, failure occurred at time t = 3.10 or approximately 47% of the building gravity load. The building drift at failure was approximately 1.4% for this case. For both loading cases failure occurred suddenly and numerous attempts to extend the solution life proved unsuccessful. Note that certain input parameters differed slightly for each direction of loading to account for the varying properties of the tension and compression "flanges" created by out-of-plane walls connecting within the lower stories.

152 4.0 -3.0 i 2.0 -2 - 0 -3.0 -4.0- I -4.0 \, ---- -60.0 -40.0 -20.0 0.0 20.0 40.0 60.0 Roof Loteral Displacement (cm) Figure 6.7: Roof lateral displacement vs. normalized base shear Due to the amount of information generated by the building analysis, only the case of loading for one of the above directions will be examined in greater detail. Except where otherwise noted, the information presented in the following sections corresponds to lateral loading causing motion towards the East direction. Loading in the opposite direction produced qualitatively similar results. 6.5.1 Global ductility measures Measures of structural ductility are usually expressed in the form of a ductility ratio, I. In general terms, ductility ratio may be defined as the ratio of the maximum (useful) deformation in the inelastic range, 6, to the deformation corresponding to a predefined yield limit, 6y, as shown by Eq. 6.4. Ductility ratio takes on various meanings depending on whether the loading is monotonic or cyclic and whether the

153 deformation considered is in the form of displacements, rotations, or curvatures. = (6.4) by In order to obtain a quantitative indication of the building overall performance for the type of loading considered, a ductility ratio defined by roof level lateral displacements is considered here. The maximum roof level displacement is set by the analysis results. However, the yield displacement may vary according to any number of yield definitions. Ductility ratios based on two different definitions yield are presented here, as shown in Fig. 6.8. In Figure 6.8(a), the yield definition has been based on energy absorption concepts. By choosing the yield definition through equating the two shaded areas, both the original curve and the idealized elastic-plastic curve absorb the same amount of energy over the loading history. Using this approach, the lateral roof displacement at yield equals 18.2 cm. (0.64% drift index) and the ductility ratio becomes i = 42.2/18.2 = 2.3. The equal energy approach is often used to mark the yield point in cases where the response clearly flattens before failure, as is the case for many frame type structures. However, wall dominant structures show no clear plateau, but rather maintain a steady hardening reponse up until failure. In such cases, the equal energy approach may push the yield point out to unreasonable displacements. Figure 6.8(b) presents an alternate approach for yield point definition which has been applied to shear wall type structures [73]. A yield load, Fy, is determined as the load level corresponding to a 1% drift over the building height. A straight line directed from the origin through the response curve at 0.75Fy will then define the yield displacement at the load level Fy. Using this method, the yield displacement becomes 14.6 cm. (0.51% drift index) which results in a ductility ratio of y = 42.2/14.6 = 2.9.

154 3.0-.0 u ' Equol Are's V) E 00 z I 0. 0.0 10.0 20.0 30.0 40.0 50.0 Roof Loteral Displacement (cm) (a) Yield Definition using Equal Energy Approach 3.0 0 ' i; E /.0, I TI / 0.0 - — '. -, - '. -. ' - -.-.- - 0.0 10.0 20.0 30.0 40.0 50.0 Roof Lateral Displocement (cm) (b) Yield Definition using Park's Method (./30 t002. 0. 005 IotLtrl Drlcmrt (m {/). ~ (b / il ei~inusn aksMt Figure 6.8: Building Displacement Ductility Ratios

155 It should be noted that the above building ductility ratios give only a general measure of building performance and provide little indication of the ductility demands required at the member level. However, considering the high values found for the building strength these ductility ratios indicate good performance. 6.5.2 Nonlinear trends The SNAC program produces a condensed listing of all state changes in the element material subroutines over the loading history. This listing is updated on a per iteration basis and therefore allows tracing of the model response within each load increment. While information retrieved within the load increment does not represent an accepted equilibrium configuration, it is helpful when following the trends leading towards/away from equilibrium. In Fig. 6.9, roof level lateral displacement versus normal base shear is again plotted, this time indicating the load levels at which selected nonlinear actions commence. For example, between times t = 0.4 and t = 0.85 all shear walls in the finite element model begin to crack at their extreme tensile fibers (for coupled walls, cracking for both wall portions is indicated). "Wall reinforcement yielding" indicates yielding of both the flexural and mesh reinforcement in the extreme tensile regions of the walls. The initiation of the various nonlinear actions is not completely uncoupled as appears in Figure 6.9. Initial cracking of the coupling beams actually ranges from t = 0.15 all the way up to t = 0.75 if the beams in the first two levels are included. Likewise, initial yielding of the coupling beams actually occurs between t = 0.9 and t = 2.0 if coupling beams in the lower two levels are considered. The behavior of the beams in the first two stories was not included because the time of their cracking/yielding appears to be isolated from that of the other beams and because

156 3. 0 -2.5 2.0 /, Wall Reinforcement Yielding 0 ' 1.5 / g 1 * / Coupling Beam Yielding.0 -1. Shear Wall Cracking 0.5 -I-y_ Coupling Beam Cracking 0.0 I i ' i ' l l l 0.0 10.0 20.0 30.0 40.0 50.0 Roof Level Laterd Displacement (cm) Figure 6.9: Nonlinear actions during loading history the lower level beams should have relatively less effect on the models global response. Additional nonlinear actions which are not indicated in Figure 6.9 include: * The second story of Wall 2-2 experiences concrete compressive nonlinearity (i.e. an effective compressive stress greater than 0.5f ) as early as t = 1.00. Initial nonlinearity of the concrete in the extreme compression fibers of the remaining walls ranges from t = 1.50 to t = 1.90. Beginning at t = 2.55 and leading up to failure a limited number of the wall elements showed concrete behavior in the post-peak compressive range. * From t = 2.35 to t = 2.75 compressive yielding of the flexural and mesh reinforcement commences in the extreme compression zones of the walls. * Cracking initiates in the slender walls modeled using the BEAM elements (Wall 3-3 and the upper portions of Walls 5-5 and 7-7) between t = 0.30 and t = 0.90.

157 Yield moments were exceeded in these walls beginning at t = 1.35 through t = 2.35. All cracking and yielding of the beam elements was limited to the lower half of the building structure even though nonlinear behavior was permitted over the entire building height for these elements. The trends leading up to structural failure are illustrated more dramatically by Figures 6.10 through 6.13, where selected nonlinear actions in the major walls are shown for times t = 0.5, 1.0, 2.0, and 2.79 (failure), respectively. Coupling beam ends in the cracked and yielded states are indicated by a circular outline and a solid circle, respectively. The presence of cracking in the continuum elements is indicated by fine lines centered at the element integration points. These crack indicators do not represent actual cracks (i.e. material discontinuities) but rather they indicate the magnitude of "crack strains" in the continuum model. That is, the length of the crack lines indicates the relative magnitude of the "crack strain" at the integration point in question. Scaling the crack indicators to ultimately extend beyond the element boundaries, while disturbing at first glance, allows one to follow crack pattern development at lower load levels. Concrete crushing in the continuum elements is indicated by a small solid circle centered at the corresponding integration point. The crack angles portrayed in the preceding four figures have not been adjusted to reflect the rotation of the structure in the deformed plot. The difference is negligible for all but the plots in Figure 6.13 where the deformation is largest. Figure 6.14 shows the nonlinear actions at t = 2.79 in an undeformed plot for comparison. The difference in crack representation between Figures 6.13 and 6.14 is not significant, but the latter figure tends to show the presence of a diagonal compression field somewhat better.

I --- —-_ 0 -- -_ () ~ )* ( ( ( ( ( C C -4) -4) I F Iq I I -A Cn 00 \ wall cracking * wall cruhing bamun cracking 0 brm yildig Figure 6.10: Wall nonlinear trends at t = 0.50

I.. i < I I6 I ---- < >~ i i O II ( 4 < 4 4 O:0 - IF -- 0- --- 0 4 ) I -- 4), 7.. I - I I C' C0 cn q~ I \ wall cracking * wall crushing 0 beam cracking * bum yilding Figure 6.11: Wall nonlinear trends at t = 1.00

\ wall cracking * wall crushing 0 beam cracking * beam yilding Figure 6.12: Wall nonlinear trends at t = 2.00

\ wall cracking * wall cruahing 0 beam cracking 0 beam ywding Figure 6.12: Wall nonlinear trends at t = 2.00

r' 4 4 4 -— i,,., Ei -\ ' ~- ~. -— Z —. N \ T --- \ wall cracking * wall crushing 0 beam cracking 0 beam yieding fligure 6.14: Wall nonlinear tlrends at t = 2.79 (unldeforined plot)

163 The final converged state for loading in the opposite direction (towards the West) occured at t = 3.10. Selected nonlinear actions in the major walls occuring at this load stage are shown in Figure 6.15. The nonlinear trends for t < 3.10 appear quite similar to those already shown in Figures 6.10 through 6.12 for loading towards the East. 6.5.3 Interstory drift indices Lateral displacement profiles over the building height are shown in Figure 6.16 for several stages of the loading history. At lower lateral load levels, the drift profile exhibits a slight reversing of curvature in the upper stories due to the influence of the coupling beams in Walls 2-2 and 9-9. This behavior is more evident in Figure 6.17 where both tangential and total interstory drift indices are plotted for several stages of the loading history. Total interstory drift indices were computed by dividing interstory displacements by the corresponding story height and then converting to percentage units by multiplying by a factor of 100. Tangential interstory drift indices are equal to the total drift indices minus the percentage of drift due to rotation of the corresponding story floor level. The story drifts indices at V = 0.05VUBC have been included as a representative of the structural model in an entirely elastic state. At lower load levels the effects of the coupling beams on the drift profile are most pronounced. At V = VUBC the maximum total drift index is slightly greater that 0.10%. This value clearly indicates the structure's stiff nature as it is well within the 0.50% and 1.50% limits set by the Uniform Building Code and the Applied Technology Council Recommendations-3 [4], respectively. The building drift index (i.e. roof displacement/building height*100) at this load level measures 0.08%. The degree of

li64: t ~ s.: / // 3/ trr -.; MS/ l 0 l /^ ' '1,6' i', I I I If i/i, I t ~ ~ ~ I I i I I \/1/ - -- * o. "ifs *o.d

165 30 25 20 Laterd Dispcement (cm) Figure 6.16: Lateral drift profiles at several load levels reversed curvature in the drift profile has somewhat diminished due to crack partial yielding of the coupling beams. At V = 2.79VUBC (failure) the maximum interstory drift and building drift indices measure approximately 2% and 1.5%, respectively. At this load level there is no discernable reverse curvature. From the Story 5 upward the tangential interstory drift appears negligible as if this portion of the structure were undergoing rigid body motion. This behavior is partly due to the elastic idealization of the upper half of the structure. 6.5.4 Neutral axis migration Figure 6.18 shows the shifting of the neutral axis during the loading history for the second story portion of Wall 4-4. The location of the neutral axis was found

166 10 8 6 4 2 -1 10 8 I1 I ~~~~~~~~ _- I I I L ] C T=rigentid Drift - Totd Drift.0 (a) 0.0 1.0 20 3.0 4.0 *105 Interstory Drift Ratios at V=0.05V, I 6 4 I I I _I I -I - ~ I I I I 2 0 -0 1] = TngeTntid rift - TO Drift ).02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 b) hterstory Drift Ratios at V=tOOVuc 10 8 I I 6 4..m ~ rT I = TC - Tc o o0.5 1 5 U U U U U U I * *. 2 0 rgentid Drift:td Drift I 2 25 (c) terstory Drift Ratios at V=279Vuc Figure 6.17: Interstory drift indices at several load levels

167 by interpolating between vertical strain values registering at the element sampling points. Due to the presence of gravity loading the neutral axis is located outside the wall section for very low lateral load levels. With increasing loading, the neutral axis tends to migrate towards the compressive end of the wall due to both the relative decreasing effect of the gravity load and the nonlinear actions in the wall (and may reverse direction after compressive softening). One would expect wall cracking and wall yielding to noticeably affect the location of the neutral axis. Within this analysis, the rigid diaphram coupling of all walls in the structural model, plus the fact that each wall cracks/yields at different stages in the load sequence, tends to smooth out the variance in the neutral axis position as plotted in Figure 6.18. Near failure, the depth to the neutral axis is only about 14% of the total wall depth. The other walls modeled using the continuum element approach showed similar behavior. The above results were entirely predictable and have been included here only to show how naturally the continuum element approach accounts for variation in the neutral axis location. One concern of researchers employing the equivalent column approach to modeling shear walls is the method's inability to account for shifting of the neutral axis during the loading history. This' deficiency may greatly influence the analysis results when walls are coupled together in a vertical sense (e.g. via coupling beams). 6.5.5 Failure state characterlzation Lateral loading on the structure was applied in increments of At = 0.05 ( or AV = 0.05VuBc) over much of the loading history. Approaching the failure load this load increment size was reduced first by a factor of 10 and then 100. The final converged state during the solution process actually occurred at t = 2.7935. For

168 1.2 1. 0 -c 0.8 -< \: -- 0. 6 -Z -o 0.4 - N 0.0,., 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Normalized Base Shear (V/Vuc) Figure 6.18: Migration of wall 4-4 neutral axis during loading -__sequent loading, the structural abruptly loses its load resisting capacity. Since loading was applied in a "force control" sense, the solution procedure was unable to regain equilibrium at t = 2.7940. Following the damage occurring within this critical load increment provides information useful in characterizing the ultimate failure mode of the structure. Figures 6.19 through 6.21 graphically portray the progressive failure of the structure during iteration cycles 2, 4, and 6, respectively. Deformations have been scaled by the same factor as in previous plots of this type. It must be emphasized that the plots shown in these figures do not represent accepted equilibrium states, rather they show the structure falling away from equilibrium. Both the information presented in Figures 6.19 through 6.21 and additional program data indicate the progressive crushing of the concrete in the compression regions of the first and second story shear walls to be the ultimate failure mechanism of the

169 structure. Crushing is prevented at various locations in the lower stories due to the confinement provided by orthogonal walls which act as boundary elements. With the exception of crushing in Wall 4-4, all crushing occurs in "unconfined" regions. Nearly identical trends were observed for lateral loading in the opposite direction (i.e. for t > 3.10 during westward lateral loading). Building drift at the final converged state, t = 2.7935, was found to be 1.55%. During the next load increment, the building drift grows rapidly from 1.55% at iteration cycle 2 to 1.78% and 2.03% for iteration cycles 4 and 6, respectively. Beyond the sixth iteration cycle, displacements become tremendously large and most of the program information becomes meaningless. 6.6 Additional Considerations 6.6.1 Concrete tensile strength Adjustments to the tensile softening relations for preserving correct fracture energy during crack propagation are valid provided the element size is smaller than a certain limiting size [17]. For large structures, this restriction can be problematic since the analysis would require a prohibitively large number of elements. Alternatively, proper fracture energy can be preserved by assuming a vertical stress drop (i.e. no tensile softening) and replacing the actual tensile strength, f, with an equivalent tensile strength, feq. This approach, as outlined in [17], was used to determine an equivalent concrete tensile strength for use in the building analyses. The resulting equivalent strength value was equal to approximately 62% of the actual tensile strength (i.e. feq = 0.62ft). The procedure entailed making various assumptions, including value of the actual concrete tensile strength, the maximum aggregate size used in the walls, and the direction of crack propagation through the mesh.

\ wall cracking * wall crushing 0 beam cracking * beam yielding Figure 6.19: Post-converged state for Iteration Cycle 2

\ wall cracking *0 wall crushing 0 beam cracking 0 beam yielding Figure 6.20: Post-converged state for Iteration Cycle 4

Figure 6.21: Post-converged state for Iteration Cycle 6

173 3.0 2.5- / - -I, 2.0 0 1.5 *.. 0.0 0.0 10.0 20.0 30.0 40.0 50.0 Roof Lateral Displocement (cm) Figure 6.22: Sensitivity of building response to concrete tensile strength Building load-deflection response using the equivalent tensile strength is compared to the response using the actual concrete tensile strength in Figure 6.22. As expected, the equivalent tensile strength model shows a slightly more flexible response. Both models predict ultimate conditions at nearly the same base shear level. 6.6.2 Floor slab flexibility The building floor slabs were assumed to act as rigid diaphrams with respect to lateral translation, but provide no resistance with respect to vertical translation and rotation. The rigid diaphram modeling was effectively carried out using constraint equations to insure each nodal point at a given floor level translated equally. However, as pointed out in Section 6.3.2, certain situations require a more realistic modeling

174 of the floor system which accounts for the slab's limited ability to transfer loading. A relatively straitforward means to account for finite slab flexibility would be to simply adjust the coefficients in the constraint equations linking 'slave' to 'master' node at the level in question. As a rough approximation, these coefficients can be determined in the following manner as indicated in Figure 6.23. Walls continuous above Level i (Wall A) will translate equally with the Level i 'master' node. Walls beginning at Level i (Wall B) will carry their 'share' of the lateral load applied at Level i, but will not immediately pick up their 'share' of the lateral load corresponding to Levels j, k...n. For lack of more specific data, it is assumed that lateral load transfer to Wall B varies linearly from zero up to a point where complete lateral load transfer occurs. The point of complete transfer is estimated by assuming the floor slab to transmit compression at a 1:1 slope away from the wall. The displacement of Wall B due to Fi, added to the displacement caused by the transfer of load due to Fj, Fk,...F, will define the constraint equation coefficients for Wall B as shown in Figure 6.23(b). This type of procedure was used to simulate finite floor slab flexibility in the regions surrounding Walls 5-5 and 7-7 which begin abruptly in the lower stories of the structure, as shown in Fig. 6.4. These modifications to the building model did not appreciably affect the overall building response, as evident in Figure 6.24. 6.6.3 Lateral load distribution profile For the Villa Real building, both UBC and ATC-3 imply the use of a triangular, or nearly triangular, acceleration profile to establish the lateral load distribution for use in the static analysis. Note that these buildings 'codes' predict the same lateral load distribution regardless of structural type so long as the building fundamental pe

175 Wall B k-x a) Walls connected by quasi-flexible floor system aOi [ U level i Umaster node level i Wall A 1.0 Fj Wall B n F i __ --- — _._ - O — Wall Length X L b) Constraint equation coefficients at level i Figure 6.23: Flexible diaphram idealization

176 15 10 Rigid Diophrom 5- * Quasi-flexible Diaphram 0- -** \ 0.0 10.0 20.0 30.0 40.0 50.0 Lateral Displocement (cm) Figure 6.24: Building response with quasi-flexible diaphram riod T, the story heights hi, and the lumped story weights wi remain the same. This seems unreasonable since the vibration mode shapes of a moment-resisting frame system may differ considerably from those of wall dominant systems. Both tests and analytical procedures investigating frame-wall and braced frame systems have shown the distribution of seismic forces over the building height can appear quite different when comparing the profile at time of maximum axial-flexural (overturning) strength demand and maximum base shear strength demands [21]. At the time of maximum base shear demand the lateral load distribution may appear closer to uniform (rectangular) over the building height. It follows that a triangular distribution profile should be a more strict check on the effects of overturning moment (flexure in the walls) while the rectangular profile should be a more severe test against shear mode failure of the structure.

177 Because of such considerations, the results due to a uniform (rectangular) load distribution over the building height have been compared to the former results using the triangular acceleration profile. As shown in Figure 6.25, the rectangular load distribution pattern appears to be a less severe test of the building lateral strength. This result was anticipated as the individual walls are somewhat slender and the failure mode from the earlier runs was clearly flexurally dominated. The lateral displacements at failure for the rectangular profile are only slightly greater in the lower levels and virtually identical at roof level compared to the UBC type load profile. Plots of the nonlinear trends at failure show similar similar cracking and coupling beam behavior and will not be duplicated here. Again, the global failure is preceded by concrete crushing at the extreme compression end of Wall 2-2 in the second story. It is remarkable that the failure mode and displacement are nearly identical for each of the assumed loading profiles. The material models inability to account for nonorthogonal secondary cracking certainly enhances the chances of a flexurally dominated failure mode. Flexural cracks form early in the response and then effectively 'lock out' the formation of diagonal shear cracking during increasing diagonal compression later in the response.

178 4.0 3.5.. 3.0 - 2.5 () 2.0 - D ' / I 1.5 0 i o - z Triaonular acceleration distribution RectonguIa lo load distribution 0. 5 0.0 I ' I.!. I. 0.0 10.0 20.0 30.0 40.0 50.0 Lateral Displacement (cm) Figure 6.25: Load distribution profile response comparison

CHAPTER VII SUMMARY 7.1 SNAC Program Development A finite element method based computer program has been developed for predicting the progressive inelastic actions leading up to the failure state of shear wall dominant buildings subjected to inplane quasistatic lateral loading. The program was named SNAC (Simple Nonlinear Analysis of reinforced Concrete structures) and is coded in standard FORTRAN 77 language. The SNAC program utilized many of the subroutines and modern programming concepts of DLEARN [50], an educationally orientated program for the dynamic linear analysis of structures. Program logic stressed clarity and ease of modification. SNAC program development was based on a macro-model approach to nonlinear reinforced concrete analysis. Within this type of approach the global structural response is of primary importance and the details of highly localized phenomena are generally neglected. Due account was given to the inelastic material behavior of reinforced concrete through application of basic material models taken from the wealth of literature on the subject. One of the most notable features of the SNAC program is the two-dimensional continuum element modeling of the shear walls in a state of plane stress. This is seen 179

1U as a necessary improvement over the traditional beam-column approach to modeling of shear wall type structures. The discrete line element approach to modeling shear walls is unable to directly account for such multi-dimensional phenomena as nonlinear strain distribution over the wall length and neutral axis migration during loading. These concerns become more critical as the wall aspect ratio becomes smaller. Beamcolumn elements have been included in the program library for modeling coupling beams and walls with high aspect ratios which lend themselves better to a discrete line element approach. The continuum element material subroutines were validated through comparison with test results. The SNAC program had success in matching certain aspects of various test results including those of unreinforced concrete fracture specimens, concrete specimens under biaxial loading, deep beams under shear type loading, and isolated shears walls subjected to lateral loading. The beam element material subroutines were developed to provide results corresponding to a section-type analysis, where the deformation state and concentration of nonlinearity over the beam length were specified a priori. 7.2 Chilean Building Analysis Observations Most of the more than 400 medium to medium-high rise buildings in the Chilean coastal city of Vifia del Mar performed exceptionally well during the March 3, 1985 earthquake experienced in that region. In response to these extraordinary circumstances, the United States National Science Foundation has sponsored an interuniversity research program to investigate the buildings of the Vifia del Mar region. Universities belonging to the research program include the University of California at Berkeley, the University of Illinois (Urbana), and The University of Michigan.

181 One of the Vina del Mar buildings, the Villa Real building, has been analyzed using the SNAC program. Inplane lateral loading applied to the analysis model was assumed to vary quasi-statically and increase proportionally over the loading history. The lateral load distribution over the building height was based on a triangular horizontal acceleration profile roughly approximating the building fundamental mode response during dynamic excitations. Due to the short period nature of the structure, an identical lateral load distribution would be specified using either the 1985 Uniform Building Code provisions or the Applied Technology Council-3 Recommendations. Additional details regarding the Villa Real building analysis are presented in Chapter VI of this report. The following paragraphs present a brief summary of the analysis findings. 1. The Villa Real building possesses considerable overstrength against lateral loading relative to the Uniform Building Code design levels. Failure of the numerical model occured at a base shear roughly three times the Uniform Building Code specified seismic base shear. Damage occurring in the building model at the UBC base shear level was limited to cracking in the coupling beams and minor cracking in the shear walls. 2. The Villa Real building appears to be quite stiff with respect to lateral load resistance. The maximum interstory drift index resulting at a base shear level V = VUBC is only 0.10%. This value is roughly an order of magnitude less than building code allowable drift limits and helps substantiate the short building fundamental period values, T = 0.5 seconds, obtained in preliminary analyses. 3. The ultimate failure mechanism of the structure proved to be the progressive crushing of the concrete in the compression regions of the first and second

182 story shear walls. Failure occurred abruptly indicating a lack of structural redundancy with respect to this type of failure. Crushing is prevented at various locations in the lower stories by the presence of orthogonal walls which act as boundary elements. With the exception of crushing in Wall 4-4, all crushing occurs in "unconfined" regions. 4. The failure mode and load level proved to be insensitive to variations in the concrete tensile strength and the degree of floor slab rigidity in the lower levels. Varying the lateral load distribution profile did not affect the ultimate failure mode of the model. However, significantly higher base shear levels were realized at failure for the rectangular load profile relative the UBC specified load profile. The preceding analysis has been for only one of the Chilean buildings. Results from the analyses of other buildings within the Vifia del Mar building inventory are anticipated. At that time generalizations can be made. 7.3 Future Considerations An equivalent lateral force method based approach to building analysis provides useful yet limited information regarding building performance during earthquake type motions. Within the program developed in this study, there is the need for more general material models which are able to account for both non-orthognal secondary cracking and compressive softening in the shear walls. Such capabilities are necessary for detecting possible diagonal shear type failure later in the loading history. Such models have been developed [24,25], but have not been implemented here. Beyond the deficiencies in the material modeling, the present analysis framework suffers from two major shortcomings. Firstly, equivalent lateral force methods are

183 unable to account for the amount of ground motion energy transfered to the structure during the loading history. This consideration becomes especially critical when investigating cases where long duration motion or structural resonance might occur. These problems can be addressed by incorporating both the inertial effects of the reactive masses and cyclic loading capability into the analysis scheme. Conceptually, this goal can be accomplished through straightforward extensions of the analysis program developed in this study. The major difficulty centers around the development of material models able to accurately account for the various actions causing stiffness degradation during cyclic loading. Secondly, two-dimensional analyses are unable to directly account for the effects of out-of-plane motion on the building response. In general, earthquake excitations contain all six compononents of ground motion. Linear analyses may somewhat account for out-of-plane motion using superposition techniques. However, such methods are not applicable to nonlinear analyses. Inevitably, building analysis programs must be formulated to perform full three-dimensional simulation of building nonlinear response to general input motion. Such capabilities would also reduce the ambiguities related to modeling the effects of structural components framing in from out-of-plane directions. Both the degree of difficultly in 3-D concrete material modeling and the increase in computational expense associated with extending the analysis another dimension are major obstacles to realizing these goals. The latter of these two concerns may be seen as only a temporary disadvantage, however, as the benifits from the current trends in supercomputing (e.g. reduced cycle times, vector and parallel processing) should soon become more widely available to researching engineers.

BIBLIOGRAPHY 184

185 BIBLIOGRAPHY [1] American Concrete Institute Committee 318 (1983): Building Code Requirements for Reinforced Concrete, ACI 318-83, Detroit. [2] Agrawal, A. B., Jaeger, L. G., and Mufti, A. A., "Crack Propagation and Plasticity of Reinforced Concrete Shear-Wall Under Monotonic and Cyclic Loading," Finite Element Methods in Engineering, Proceedings of the 1976 International Conference on Finite Element Methods in Engineering, Y. K. Cheung and S. G. Hutton, eds., University of Adelaide, Australia, 1976. [3] Al-Nahlawi, K., "Shear Strength of Reinforced Concrete Beams with Low Reinforcing Ratios," Ph.D. Dissertation (in preparation), Department of Civil Engineering, The University of Michigan, Ann Arbor, MI, 1989. [4] Applied Technology Council, "Tentative Provisions for the Development of Seismic Regulations for Buildings," Report ATC 3-06, NBS Special Publication 510, NSF Publication 78-08, June 1978. [5] Bathe, K. J., "ADINA - A Finite Element Program for Automatic Dynamic Incremental Nonlinear Analysis," Report 82448-1, Acoustics and Vibration Laboratory, Department of Mechanical Engineering, MIT, Cambridge, MA, 1975, rev. 1978. [6] Bathe, K. J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1982, 735 pp. [7] Bathe, K. J., and Cimento, A. P., "Some Practical Procedures for the Solution of Nonlinear Finite Element Equations," Computer Methods in Applied Mechanics and Engineering 22, North-Holland Publishing Co., 1980, pp. 59-85. [8] Bathe, K. J., and Ramaswamy S., "On Three-Dimensional Nonlinear Analysis of Concrete Structures," Nuclear Engineering and Design 52, North-Holland Publishing Co., i979, pp. 385-409. [9] Bazant, Z. P., "Instability, Ductility and Size Effects in Strain Softening Concrete," Journal of the Engineering Mechanics Division, ASCE, Vol. 102, No. EM2, April 1976, pp. 331-344. [10] Bazant, Z. P., "Comment on Orthotropic Models for Concrete and Geomaterials," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 109, No. 3, June 1983, pp. 577-593.

186 [11 Bazant, Z. P., "Size Effect in Blunt Fracture: Concrete, Rock, Metal," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 110, April 1984, pp. 518-535. [12] Bazant, Z. P., and Bhat, P. D., "Endochronic Theory of Inelasticity and Failure of Concrete," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 102, No. EM4, August 1976, pp. 701-722. [13] Bazant, Z. P., and Cedolin, L., "Blunt Crack Band Propagation in Finte Element Analysis," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 105, No. EM2, April 1979, pp. 297-315. [14] Bazant, Z. P., and Cedolin, L., "Fracture Mechanics of Reinforced Concrete," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 106, No. EM6, October 1979, pp. 1287-1306. [15] Bazant, Z. P., and Cedolin, L., "Finite Element Modeling of Crack Band Propagation," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 109, No. 1, 1983, pp. 69-92. [16] Bazant, Z. P., and Kim, S. S., "Plastic Fracturing Theory for Concrete," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 105, No. EM3, June 1979, pp. 407-428. [17] Bazant, Z. P., and Oh, B. H., "Crack Band Theory for Fracture of Concrete," Materials and Structures, (RILEM, Paris), Vol. 16, 1983, pp. 155-177. [18] Bazant, Z. P., and Oh, B. H., "Rock Fracture Via Strain Softening Finite Elements," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 110, No. 7, July 1984, pp. 1015-1035. [19] Bazant, Z. P., and Tsubaki, T., "Slip-Dilitancy Model for Cracked Reinforced Concrete," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 106, No. ST9, September 1980, pp. 1947-1966. [20] Bergan, P. G., and Holland, I., "Nonlinear Finite Element Analysis of Concrete Structures," Computer Methods in Applied Mechanics and Engineering 17/18, North-Holland Publishing Co., 1979, pp. 443-467. [21] Bertero, V. V., "Lessons Learned from Recent Earthquakes and Research and Implications for Earthquake-Resistant Design of Building Structures in the United States," Earthquake Spectra, Vol. 2, No. 4, 1986, pp. 825-858. [22] Bicanic, N., and Zienkiewicz, O. C., "Constitutive Model for Concrete Under Dynamic Loading," Earthquake Engineering and Structural Dynamics, Vol. 11, 1983, pp. 689-710. [23] Bolander, J., Feng, S. P., and Fustok, M., "Analysis of the Villa Real Building," (Report summarizing a preliminary study of the Villa Real Building),

187 Department of Civil Engineering, The University of Michigan, Ann Arbor, MI, February, 1987. [24] de Borst, R., "Smeared Cracking, Plasticity, Creep, and Thermal Loading - A Unified Approach," Computer Methods in Applied Mechanics and Engineering 62, North-Holland Publishing Co. 1987, pp. 89-110. [25] de Borst, R., and Nauta, P., "Non-orthogonal Cracks in a Smeared Finite Element Model," Engineering Computing, Pineridge Press Ltd., Vol. 2, March 1985, pp. 35-46. [26] Broek, D., Elementary Engineering Fracture Mechanics, Second Edition, Sijthoff and Noordhoff, 1978, 437 pp. [27] Buyukozturk, O., "Nonlinear Analysis of Reinforced Concrete Structures," Computers and Structures, Vol. 7, 1977, pp. 149-156. [28] Cantin, G., "An Equation Solver of Very Large Capacity," Internation Journal for Numeric Methods in Engineering, Vol. 3, 1971, pp. 379-388. [29] Cedolin, L., and Nilson, A. H., "A Convergence Study of Iterative Methods Applied to Finite Element Analysis of Reinforced Concrete," Internation Journal for Numeric Methods in Engineering, Vol. 12, 1978, pp. 437-451. [30] Cedolin, L., and Bazant, Z. P., "Effect of Finite Element Choice in Blunt Crack Band Analysis," Computer Methods in Applied Mechanics and Engineering 24, North-Holland Publishing Co. 1980, pp. 305-316. [31] Chen, P. F. S., and Powell, G. H., "Generalized Plastic Hinge Concepts for 3D Beam-Column Elements," Earthquake Engineering Research Center Report No. UBC/EERC 82-20, University of California, Berkeley, CA, November, 1982. [32] Chen, W. F., Plasticity in Reinforced Concrete, McGraw Hill Co., New York, NY, 1982, 474 pp. [33] Chen, A. T. C., and Chen, W. F., "Constitutive Relations for Concrete," Journal of the Engineering Mechanics Division, ASCE, Vol. 101, EM4, August 1975, pp. 465-481. [34] Crisfield, M. A., "A Fast Incremental/Iterative Solution Procedure that Handles 'Snap-Through'," Computers and Structures, Vol. 13, Pergamon Press, Great Britain, 1981, pp. 55-62. [35] Darwin, D., "Concrete Crack Propagation - Study of Model Parameters," U. S. Japan Seminar on Finite Element Analysis of Reinforced Concrete Structures, ASCE, New York, NY, 1985, pp. 184-203. [36] Darwin, D., and Pecknold, D. A., "Analysis of RC Shear Panels Under Cyclic Loading," Journal of the Structural Division, ASCE, Vol. 102, ST2, February 1976, pp. 355-369.

188 [37] Darwin, D., and Pecknold, D. A., "Nonlinear Biaxial Stress-Strain Law for Concrete," Journal of the Engineering Mechanics Division, ASCE, Vol. 103, EM2, April 1977, pp. 229-421. [38] DiTommaso, A., "Evaluation of Concrete Fracture," Fracture Mechanics of Concrete: Material Characterization and Testing, A. Carpinteri and A. R. Ingraffea (eds.), Martinus Nijhoff Publishers, The Hague, 1984. [39] Gerstle, K. H., Aschl, H., Bellotti, R., Bertacchi, P., Kotsovos, M. D., Ko, H. Y., Linse, D., Newman, J. B., Rossi, P., Schickert, G., Taylor, M. A., Traina, L. A., Winkler, H., and Zimmerman, R. M., "Behavior of Concrete Under Multiaxial Stress States," Journal of the Engineering Mechanics Division, ASCE, Vol. 106, EM6,December 1980, pp. 1383-1403. [40] Gerstle, K. H., Aschl, H., Bellotti, R., Bertacchi, P., Kotsovos, M. D., Ko, H. Y., Linse, D., Newman, J. B., Rossi, P., Schickert, G., Taylor, M. A., Traina, L. A., Winkler, H., and Zimmerman, R. M., "Strength of Concrete Under MultiAxial Stress States," Proceedings, Douglas McHenry International Symposium on Concrete and Concrete Structures, ACI Publication SP-55, American Concrete Institute, Detroit, 1978, pp. 103-131. [41] Glemberg, R., and Samuelsson, A., "A General Constitutive Model for Concrete Structures," Publication 83:4, Department of Structural Mechanics, Chalmers University of Technology, G6teborg, Sweden, 1983. [42] Glucklich, J., "Fracture of Plain Concrete,"- Journal of the Engineering Mechanics Division, ASCE, Vol. 89, EM6,December 1963, pp. 127-137. [43] Grootenboer, H. J., Leijten, S. F. C. H., and Blaauwendraad, J., "Numerical Models for Reinforced Concrete Structures in Plane Stress," Heron, Vol. 26, No. Ic, 1981. [44] Hall, W. J., "Earthquake Engineering Guidelines and Special Considerations," Proceedings of Earthquake and Earthquake Engineering: The Eastern United States, Vol. 1, pp. 53-84, Knoxville, TN, September 14-16, 1981. [45] Han, D. J., and Chen, W. F., "A Nonuniform Hardening Plasticity Model for Concrete Materials," Mechanics of Materials 4, North-Holland, 1985, pp. 283 -302. [46] Han, D. J., and Chen, W. F., "Constitutive modeling in Analysis of Concrete Structures," Journal of the Engineering Mechanics Division, ASCE, Vol. 113, No. 4, April 1987, pp. 577-593. [47] Hegemier, G. A., and Read, H. E., "On Deformation and Failure of Brittle Solids: Some Outstanding Issues," Mechanics of Materials 4, North-Holland, 1985, pp. 215-259.

189 [48] Hillerborg, A., Modeer, M., and Petersson, P. E., "Analysis of Crack Formation and Crack Growth in Concrete by Means of Fracture Mechanics and Finite Elements," Cement and Concrete Research, Vol. 6, 1976, pp. 773-782. [49] Hillemeier, B., and Hilsdorf, H. K., "Fracture Mechanics Studies on Concrete Compounds," Journal of Cement and Concrete Research,.Vol. 7, 1977, pp. 523 -536. [50] Hughes, T.J.R., The Finite Element Method - Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1987, 803 pp. [51] Kachanov, L. M., Fundamentals of the Theory of Plasticity, MIR Publishers, Moscow, 1974. [52] Knott, K. F., Fundamentals of Fracture Mechanics, Butterworths, London, England, 1973. [53] Krajcinovic, D., and Fonseka, G. U., "The Continuous Damage Theory of Brittle Materials: Part 1 General Theory," Journal of Applied Mechanics, Vol. 48, December 1981, pp. 809-815. [54] Krajcinovic, D., discussion of "A Nonuniform Hardening Plasticity Model for Concrete Materials," by Han and Chen, Mechanics of Materials 4, NorthHolland, 1985, pp. 309-310. [55] Kupfer, H. R., and Gerstle, K. H., "Behavior of Concrete Under Biaxial Stresses," Journal of the Engineering Mechanics Division, ASCE, Vol. 99, EM4, August 1973, pp. 852-866. [56] Kupfer, H. R., Hilsdorf, H. K., and Rusch, H., "Behavior of Concrete Under Biaxial Stresses," American Concrete Institute Journal, Vol. 66, No. 8, August 1969, pp. 656-666. [57] Leibengood, L. D., Darwin, D., and Dodds, R. H., "Parameters Affecting FE Analysis of Concrete Structures," Journal of the Structural Division, ASCE, Vol. 112, No. 2, February 1986, pp. 326-341. [58] Liu, T. C. Y., Nilson, A. H., and Slate, F. O., "Stress-Strain Response and Fracture of Concrete in Uniaxial and Biaxial Compression," American Concrete Institute Journal, Vol. 69, No. 5, May 1972, pp. 291-295. [59] Mahin, S. A., and Bertero, V. V., "Rate of Loading Effects on Uncracked and Repaired Reinforced Concrete Members," Earthquake Engineering Research Center Report No. EERC 72-9, University of California, Berkeley, CA, 1972. [60] Malvern, L. E., Introduction to the Mechanics of a Continuous Medium, Prentice Hall, Inc., Englewood Cliffs, NJ, 1969. [61] Melosh, R. J., and Bamford, R. M. "Efficient Solution of Load Deflection Equations," Journal of the Structural Division, ASCE, Vol. 95, No. ST4, April 1969, pp. 661-676.

190 [62] Mondkar, D. P., and Powell, G. H., "Towards Optimal In-Core Equation Solving," Computers and Structures, Vol. 4, Pergamon Press, Great Britain, 1974, pp. 531-548. [63] Mondkar, D. P., and Powell, G. H., "Large Capacity Equation Solver for Structural Analysis," Computers and Structures, Vol. 4; Pergamon Press, Great Britain, 1974, pp. 699-728. [64] Nayak, G. C., and Zienkiewicz, 0. C., "Note on Alpha-Constant Stiffness Method for the Analysis of Nonlinear Problems," Internation Journal for Numeric Methods in Engineering, Vol. 4, 1972, pp. 579-582. [65] Nayak, G. C., and Zienkiewicz, 0. C., "Elasto-plastic Stress Analysis. A Generalization for Various Constitutive Relations Including Strain Softening," Internation Journal for Numeric Methods in Engineering, Vol. 5, 1972, pp. 113-135. [66] Nelissen, L. J. M., "Biaxial Testing of Normal Concrete," Heron, The Netherlands, Vol. 18, No. 1, 1972. [67] Ngo, D., and Scordelis, A. C., "Finite Element Analysis of Reinforced Concrete Beams," American Concrete Institute Journal, Vol. 64, No. 3, March 1967, pp. 152-163. [68] Ortiz, M., "A Constitutive Theory for the Inelastic Behavior of Concrete," Mechanics of Materials 4, North-Holland, 1985, pp. 67-93. [69] Oesterle, R. G., Aristizabal-Ochoa, J. D., Fiorato, A. E., Russell, H. G., and Corley, W. G., "Earthquake Resistance Structural Walls - Tests of Isolated Walls - Phase II," PCA R/D Series 1629, Portland Cement Association, Skokie, IL, October 1979. [70] Owen, D. R. J., Prakash, A., and Zienkiewicz, 0. C., "Finite Element Analysis of Nonlinear Composite Materials by Use of Overlay Systems," Computers and Structures, Vol. 4, 1974, pp. 1251-1267. [71] Owen, D. R. J., Figueiras, J. A., and Damjanic, F., "Finite Element Analysis of Reinforced and Prestressed Concrete Structures Including Thermal Loading," Computer Methods in Applied Mechanics and Engineering 41, North-Holland Publishing Co. 1983, pp. 323-366. [72] Paris, P. C., and Sih, G. C., "Stress Analysis of Cracks," ASTM Special Technical Publication No. 381, 1965. [73] Park, R., "Ductility Evaluation from Laboratory and Analytical Testing," Stateof-the-Art Report presented at the Ninth World Conference on Earthquake Engineering, Tokyo-Kyoto, Japan, August 2-9, 1988. [74] Paulay, T., "The Design of Ductile Reinforced Concrete Structural Walls for Earthquake Resistance," Earthquake Spectra, Vol. 2, No. 4, 1986.

191 [75] Powell, G. H., and Chen, P. F. S., "3D Beam-Column Element with Generalized Plastic Hinges," Journal of the Engineering Mechanics Division, ASCE Proceedings, Vol. 112, No. 7, July 1986, pp. 627-641. [76] Przemieniecki, J. S., "Matrix Structural Analysis of Substructures," AIAA Journal, Vol. 1, 1963, pp. 138-147. [77] Rashid, Y. R., "Analysis of Prestressed Concrete Pressure Vessels," Nuclear Engineering and Design, Vol. 7, No. 4, April 1968, pp. 334-344. [78] Rice, J. R., "Mathematical Analysis in the Mechanics of Fracture," Fracture, An Advanced Treatise, H. Liebowitz, ed. Vol. 2, Academic Press, New York, 1968, pp. 191-250. [79] Riks, E., "An Incremental Approach to the Solution of Snapping and Buckling Problems," International Journal of Numeric Methods in Engineering, Vol. 17, 1979, pp. 529-551. [80] Riddell, R. R., Wood, S. L., and De la Llera, J. C., " The 1985 Chile Earthquake: Structural Characteristics and Damage Statistics for the Building Inventory in Vina del Mar," Civil Engineering Studies, Structural Research Series No. 534, University of Illinois, Urbana, April 1987. [81] Rots, J. G., Nauta, P., Kusters, G. M. A., and Blaauwendraad, J., "Smeared Crack Approach and Fracture Localization in Concrete," Heron, Vol. 30, No. 1, 1985. [82] Rots, J. G., "Fracture Toughness and Fracture Energy of Concrete," Proceedings of the International Conference of Fracture Mechanics of Concrete, Lausanne, Switzerland, October, 1985, F. H. Wittman, ed., Elsevier Science Publishers B. V., Amsterdam, 1986. [83] Saenz, L. P., discussion of "Equation for Stress-Strain Curve of Concrete," by Desayi and Krishnan, American Concrete Institute Journal, Vol. 61, No. 9, September 1964, pp. 1129-1235. [84] Simons, J. W., and Powell, G. H., "Solution Strategies for Statically Loaded Nonlinear Structures," Earthquake Engineering Research Center Report No. UBC/EERC 82-22, University of California, Berkeley, CA, November, 1982. [85] Suidan, M., and Schnobrich, W. C., "Finite Element Analysis of Reinforced Concrete," Journal of the Structural Division, ASCE, Vol. 99, ST10, October 1973, pp. 2109-2122. [86] Schnobrich, W. C., "Behavior of Reinforced Concrete Structures Predicted by the Finite Element Method," Computers and Structures, Vol. 7, Pergamon Press, Great Britain, 1977, pp. 365-376. [87] Strang, G., Introduction to Applied Mathematics, Wellesley Cambridge Press, Wellesley, MA, 1986.

UNIVERSITY OF MICHIGAN 192 3 9015 02229 1051 [88] Task Committee on Finite Element Analysis of Reinforced Concrete Structures, State of the Art Report on Finite Element Analysis of Reinforced Concrete, American Society of Civil Engineers, New York, NY, 1982, 545 pp. [89] Tasuji, M. E., Slate, F. 0., and Nilson, A. H., "Stress-Strain Response and Fracture of Concrete in Biaxial Loading," American Concrete Institute Journal, Vol. 75, No. 7, July 1975, pp. 306-312. [90] "The Chile Earthquake of March 3, 1985," Earthquake Spectra, Vol. 2, No. 2, 1986. [91] -Uniform Building Code, 1985 Edition, International Conference of Building Officials, Whittier, Ca. [92] U. S. Japan Seminar on Finite Element Analysis of Reinforced Concrete Structures, American Society of Civil Engineers, New York, May 1985. [93] Wilson, E. L., Bathe, K. J., and Doherty, W. P., "Direct Solution of Large Systems of Linear Equations," Computers and Structures, Vol. 4, Pergamon Press, Great Britain, 1974, pp. 363-372. [94] Wood, S. L., Wight, J. K., and Moehle, J. P., "The 1985 Chile Earthquake: Observations on Earthquake-Resistant Construction in Vina del Mar," Civil Engineering Studies, Structural Research Series No. 532,_ University of Illinois, Urbana, February 1987. [95] Yuzugullu, 0., and Schnobrich, W. C., "A Numerical Proceedure for the Determination of the Behavior of a Shear Wall Frame System," American Concrete Institute Journal, Vol. 70, July 1973, pp. 474-479. [96] Zienkiewitz, 0. C., The Finite Element Method, Third Edition, McGraw-Hill Book Co. New York, 1977. [97] Zienkiewitz, 0. C., Villiappan, S., and King, I. P., "Elasto-Plastic Solutions of Engineering Problems 'Initial Stress' Finite Element Approach," International Journal for Numeric Methods in Engineering, Vol. 1, 1969, pp.75-100. [98] Zaitsev, Y. V., and Kovler, K. L., "Notch Sensitivity of Concrete and Size Effect, Part 1: Effect of Specimen Size and Crack Length by 3-Point Bending," Journal of Cement and Concrete Research, Vol. 15, 1985, pp. 979-987.