Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UF00090202/00001
## Material Information- Title:
- Optimization of reinforced concrete frames using integrated analysis and reliability
- Series Title:
- Optimization of reinforced concrete frames using integrated analysis and reliability
- Creator:
- Soeiro, Alfredo V.,
- Place of Publication:
- Gainesville FL
- Publisher:
- University of Florida
- Publication Date:
- 1989
## Subjects- Subjects / Keywords:
- Data lines ( jstor )
Design optimization ( jstor ) Geometric lines ( jstor ) Lagrangian function ( jstor ) Linear programming ( jstor ) Mathematical variables ( jstor ) Mechanical systems ( jstor ) Steels ( jstor ) Stiffness ( jstor ) Subroutines ( jstor )
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Alfredo V. Soeiro. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 001515144 ( alephbibnum )
21969366 ( oclc )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

OPTIMIZATION OF REINFORCED CONCRETE FRAMES USING INTEGRATED ANALYSIS AND RELIABILITY By ALFREDO V. SOEIRO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1989 ACKNOWLEDGEMENTS I want to express my sincerest gratitude to Dr. Marc Hoit for monitoring my research, for the inventive ideas and for his constant support. I am specially thankful to Dr. Fernando Fagundo for the productive discussions, his friendship and his vigorous encouragement. I owe to these two my best recollections from the University of Florida. My sincerest appreciation is extended to Dr. Prabhat Hajela for the teachings and the careful reading of my dissertation. My indebtment goes also to Dr. Clifford Hays and Dr. John Lybas for the useful conversations and their activity in my committee. I am also grateful to Dr. David Bloomquist for his support to initiate my geotechnical reliability research and his continuous disposition to help. I would also like to acknowledge my appreciation to the Fulbright Comission and to the Department of Civil Engineering of the University of Florida for their financial aid and the opportunity to study the fascinating area of structural optimization. My thanks go also to the Department of Civil Engineering of the University of Porto, specially to Dr. Adao da Fonseca, for giving me the possibility to research and study in the USA. My sincere appreciation and best remembrances go to my friends in the Gainesville Portuguese community and to my colleagues Jose, Joon, Lin and Prasit that helped smooth the life contours created by the research work. Finally, my gratitude goes to my wife, Paula, for her work, her patience and her support throughout the whole period during which this dissertation was completed. iii TABLE OF CONTENTS Page ACKNOWLEDGEMENTS. ................................... ii LIST OF TABLES. .......... ... ..... . ...... ....... vi LIST OF FIGURES.. .. .............. ............... vii ABSTRACT ............. .. ........ ..... .. .. . .. Viii CHAPTERS 1 STRUCTURAL OPTIMIZATION....................... 1 Introduction.......................... ....... 1 Historical Background....................... 3 Methods ...... ............................... 6 Typical Applications......................... 8 Study Objectives.............................. 15 Summary... ..................... ............. 17 2 INTEGRATED OPTIMIZATION OF LINEAR FRAMES...... 21 Original Research........................... 21 Augmented Lagrangian Function................ 22 Unconstrained Minimization Techniques......... 25 Final Results........................ ..... .... 28 Further Improvements.......................... 32 3 NONLINEAR REINFORCED CONCRETE ELEMENT......... 37 Introduction ........................ ....... 37 Element Modeling Survey...................... 38 Beam Element with Inelastic Hinges............ 40 Beam Element Stiffness........................ 49 4 STRUCTURAL ELEMENT RELIABILITY ............... 54 Introduction............... . .. .. .. .. 54 Two Dimensional Space Example................ 60 Reinforced Concrete Element Reliability....... 69 5 SYSTEM RELIABILITY.......................... 74 Introduction.... .............................. 74 System Reliability and Optimization........... 75 Methods................... ................ '77 Generation of Failure Modes................. 82 Beta Unzipping Method ....................... 90 Page 6 PROCEDURE IMPLEMENTATION..................... 97 Introduction .......................... ... 97 Augmented Lagrangian Formulation.............. 98 Generalized Reduced Gradient.................. 108 Reliability................................ 114 7 EXAMPLES....................................... 119 Introduction......................... ..... 119 Result Verification.......................... 120 Debug Frame............................ ....... 121 Compared Frame............. .................... 131 Building Frame................................ 136 8 CONCLUSIONS AND RECOMMENDATIONS............... 139 Linear Material Behavior...................... 139 Nonlinear Material Behavior................... 141 Future Work ................................ 142 APPENDICES A AUGMENTED LAGRANGIAN SUBROUTINES.............. 145 B GENERALIZED REDUCED GRADIENT EXAMPLE.......... 189 C GENERALIZED REDUCED GRADIENT SUBROUTINES...... 195 REFERENCES......................... ............ 230 BIOGRAPHICAL SKETCH................................... 236 LIST OF TABLES Table Page 7.1. Debug frame (GRG): linear version results......... 124 7.2. Debug frame: Augmented Lagrangian version.......... 126 7.3. Debug frame (GRG): yielding stiffness results..... 127 7.4. Debug frame (GRG): secant stiffness results....... 129 7.5. Debug frame: element moments...................... 130 7.6. Compared frame: initial steel area reinforcement...................... 133 7.7. Compared frame results............................ 135 7.8. Building frame results............................ 138 LIST OF FIGURES Figure Pae 1.1. Implicit optimization............................ 5 1.2. Element optimization........................... 10 1.3. Truss optimization ............................. 11 1.4. System optimization ............................. 13 1.5. Geometry optimization ...... .................... 14 2.1. Pattern Search......... ....... .................. 27 2.2. Cantilever beam .................................. 29 2.3. One bay frame ................................... 31 2.4. Gradient method.................................. 34 3.1. Element model................. .............. ... 41 3.2. Material behavior.............. .................. 43 3.3. Reinforced concrete section...................... 45 3.4. Element deformation diagrams..................... 48 3.5. Curvature integration............................ 50 3.6. Secant spring stiffness.......................... 52 4.1. Design safety region............................. 61 4.2. Probabilistic functions.......................... 64 4.3. Safety checks ...................................... 66 4.4. Reliability index................................. 68 5.1. System models .................... ............... 78 5.2. Failure graph ................ .............. ... 83 5.3. Element displacements definition................. 85 5.4. System failure modes ............................ 91 5.5. Combinatorial tree............................... 96 6.1. Augmented lagrangian function plot............... 104 6.2. Augmented Lagrangian version flowchart........... 106 6.3. Generalized Reduced Gradient version flowchart... 113 6.4. Bilinear elastic-plastic diagram................. 117 7.1. Displacement verification ....................... 122 7.2. Debug frame ................................ .......... 123 7.3. Compared frame................................... 132 7.4. Building frame................................... 137 B.1. Integrated optimization example.................. 191 vii Abstract of the Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy OPTIMIZATION OF REINFORCED CONCRETE FRAMES USING INTEGRATED ANALYSIS AND RELIABILITY By ALFREDO V. SOEIRO August 1989 Chairman: Dr. Marc I. Hoit Cochairman: Dr. Fernando E. Fagundo Major Department: Civil Engineering Simultaneous analysis and design were considered in the optimization of reinforced concrete frames. Frame elements had rectangular cross sections with double steel reinforcement. Design variables were the section dimensions, the area of steel reinforcement and the structure global displacements. Equality constraints were the equilibrium equations and inequality constraints were generated by element reliability requirements, code reinforcement ratios and section dimension bounds. Optimization strategies were based on the Augmented Lagrangian formulation and on the Generalized Reduced Gradient method. Reliability of the frames was considered at the element and system level. An element failure function was defined using moment forces and flexural strength. The random viii variables considered were flexural strength of concrete and external loads. System reliability was evaluated at the mechanism level using combinations of the elementary failure mechanisms. Optimization of the frames considering material nonlinear behavior was also investigated. Inclusion of this property was performed using a one-component model for the reinforced concrete element. Inelastic rotational springs were added to the ends of the linear elastic element. The element matrix was obtained by condensation of element elastic stiffness and secant spring stiffness. Three frames were researched. Respective results using linear material behavior were discussed. In these three cases the optimal solutions were found. Element reliability constraints were active and system reliability was satisfied. The integrated formulation was validated in the linear behavior range. The nonlinear material behavior results were presented for the smaller frame. CHAPTER 1 STRUCTURAL OPTIMIZATION Introduction Optimization is a state of mind that is always implicitly present in the structural engineering process. From experience engineers learn to recognize good initial dimension ratios so that their preliminary designs demand small changes through the iterative process and that elements are not overdesigned. The motivation behind this attitude is to create a structure that for given purposes is simultaneously useful and economic. Structural optimization theory tries to rationalize this methodology for several reasons. The main one is to reduce the design time, specially for repetitive projects. It provides a systematized logical design procedure and yields some design improvement over conventional methods. It tries to avoid bias due to engineering intuition and experience. It also increases the possibility of obtaining improved designs and requires a minimal amount of human- machine interaction. There are, however, some limitations and disadvantages when using design optimization techniques. The first one is the increase in computational time when the number of design variables becomes large. Another disadvantage is that the applicability of the specific analysis program that results from the optimization formulation is generally limited to the particular purpose to which it was developed. A common inconvenience is that conceptual errors and incomplete formulations are frequent. Another drawback is that most optimization algorithms have difficulty in dealing with nonlinear and discontinuous functions and, hence, caution must be exercised when formulating the design problem. Another factor of concern is that the optimization algorithm does not guarantee convergence to the global optimum design; yielding on most occasions local optimum points. These facts lead to the conclusion that optimization results may often be misleading and, therefore, should always be examined. Therefore, some authors suggest that the word "optimization" in structural design should be replaced by "design improvement" as a better expression to materialize the root and outcome of this structural design activity (1). Nevertheless, there is an increasing recognition that it is a convenient and valuable tool to improve structural designs has been increasing among the designers community. Historical Background Throughout time there have been various attempts to address structural optimization. The earliest ideas of optimum design can be found in Galileo's works concerning the bending strength of beams. Other eminent scientists like Bernouilli, Lagrange, Young, worked on structural optimum design based on applied mechanics concepts (2). These pioneering attempts were based on a close relation to the thoughts and accomplishments of structural mechanics. They started with hypotheses of stress distribution in flexural elements and ended with material fatigue laws. The accepted first work in structural optimization discusses layout theory, or structural topology. The paper focused on the grouping of truss bars that creates the minimum weight structure for a given set of loads and materials. The author of this primary achievement was Maxwell, in 1854, and Michell developed and publicized these concepts in 1904 (3). The practical application of these theorems was never accomplished since significant constraints were not included in the original works. Some procedures widely used by structural designers are nothing more than techniques of structural optimization. A well known example is the so-called Magnel's diagram (4). It is used to find the optimal eccentricity of the cable that leads to the smallest prestressing force without exceeding the limits imposed on the stresses in prestressed concrete beams with excess capacity. This is a typical maximization problem in a linear design space, where the design variables are the eccentricity and the inverse of the cable prestressing force. The objective function is the value of the inverse of the cable prestressing force, and is to be maximized. The constraints represent the allowable stresses in tension and compression at the top and bottom of the cross-section. The problem is solved using a graphic representation of the problem, as shown in Figure 1.1, but could be solved numerically using the Simplex method. Numerical optimization methods and techniques have been widely researched and used in the operations research area, commonly known as Mathematical Programming. The practical application of these theories has been carried out in several areas for some decades like management, economic analysis, warfare, and industrial production. Lucien Schmit was the first to use nonlinear programming techniques in structural engineering design (5). The main purpose of structural optimization methods was to supply an automated tool to help the designer distribute scanty resources. Presently, anyone who wants to consider optimum structural design must become familiar with recent synthesis approaches as well as with accepted analysis procedures. 5 Magnel' s Diagram Optimum Pair P-e optimum Feasible region P Initial prestressing force: e Eccentricity of cable; e*- maximum cable eccentricity; a),b) minimum 1/P; c).d) maximum 1/P. Figure 1.1 Implicit optimization. 1/P Methods In the last twenty years researchers have made considerable advances in developing techniques of optimum design. Research and exploration of these methods were mainly developed in the aeronautical and mechanical industries, where the need for more economical and efficient final products was extremely important. More recently, with the availability of increasing computer capabilities, civil engineering researchers and designers have increased their participation in structural optimization following the lines defined by the other engineering disciplines. Optimization methods are, nevertheless, common to these different engineering design areas and are mainly divided in two groups. These are commonly known by the names Optimality Criteria and Mathematical Programming (6). Another area in structural optimization researched by a few scientists is based on duality theory concepts, and is an attempt to unify the two basic methodologies (7). Optimality Criteria methods are based on an iterative approach where the conditions for an optimum solution are previously defined. The concept can be used as the basis for the selection of a structure with minimum volume. This methodology derives from the extreme principles of structural mechanics and has been limited to simple structural forms and loading conditions. The formulation can be mathematically expressed as follows: xk+l = (p (xk,uk+l) where x is the vector of design variables, uk+ is an estimative of lagrangian multipliers and p is an adequate recurrence relation. Estimation of the lagrangian multipliers is made using the active constraints, those inequality or equality constraints with value close to zero. Recurrence relation p and lagrangian multipliers represent the necessary conditions for optimality known as Kuhn-Tucker conditions. On the other hand, the Mathematical Programming approach establishes an iterative method that updates the search direction. It seeks the maximum or minimum of multivariable function subject to limitations expressed by constraint functions. The iterative procedure may be defined as follows: Xk+l = xk + ak dk where ak is the step size and dk is the search direction. The search direction is obtained through an analysis of the optimization problem and the step size depends on the one- dimensional search along that direction. Methods of the second class may be divided in two areas. These areas are transformation methods, like penalty functions, barrier functions and method of multipliers, and primal methods, such as sequential linear and quadratic programming, gradient projection method, generalized reduced gradient and method of feasible directions. Typical Applications In structural optimal design applications there are several types of problems. They address different targets in structural design such as the best configuration for a truss or the cross sections of a prestressed concrete beam. There are four main properties of any structure that may be focused by structural optimization. These are mechanical or physical properties of the material, topology of the structure, geometric layout of the structure and cross- sectional dimensions. Main types of applications are optimization of elements, truss bars, flexural systems, continuum systems, geometry and topology (8). In the case dealing with element optimization, the search is done with a reduced number of variables and the use of code provisions transformed adequately to the optimization formulation. Element forces are found, element cross section is optimized, updated element forces are computed and the process is repeated until there is convergence. For instance, the optimal design of steel wide flange sections may have as design variables the width and thickness of the web and flanges. Constraints may be obtained in an explicit form, as the evaluation of the objective and constraint functions does not require matrix structural analysis. The minimization technique may be chosen as any one of the available direct search methods (9). Examples of design variables in element optimization are presented in Figure 1.2. Optimization of truss bar sections has been thoroughly studied due to the simplicity of truss structural optimization problems. There is a decline of interest since they are now rarely used in present structural engineering. Each bar is represented by one variable and the global stiffness matrix terms are linear functions of these variables. Of the various improvement techniques one is based on variable linking, consequently reducing the size of the problem. Another technique to decrease the size is based on constraint deletion, where inactive constraints are temporarily kept out of the optimization process. There are various formulations for the analysis model based on plastic analysis, force or displacement method (10). An example of the formulation used for truss optimization is presented in Figure 1.3. The problem of system optimization is commonly addressed using design sensitivity analysis and explicit approximations of constraint functions. The intent is to improve the performance of the chosen algorithm. Design sensitivity analysis is the calculation of the analytical derivatives of the objective and constraint functions with respect to the design variables. This information about the change in the value of a constraint related to the changes STEEL SECTION DESIGN VARIABLES Flange width Flange thickness Web height Web thickness CONCRETE SECTION DESIGN VARIABLES Width Height Top reinforcement Bottom reinforcement Figure 1.2. Element optimization. M ///~///// 15 ft 5 ft 12 ft 10 ft 10 ft 12 ft 5 t Minimize Z LiAi subject to Fi < Fc Fi < Ft where Li length of truss bar i Ai area of truss bar i Fi stress in truss bar i Fc allowable compressive stress Ft allowable tensile stress Figure 1.3. Truss optimization. in the design variables, contributes to the reduction of the exact analyses required during the optimization process. Explicit approximations of the constraint functions using a first order Taylor series expansion are widely used in Optimality Criteria and Mathematical Programming methods. In large and continuum systems some other techniques are used. For example, the sequential optimization of substructures or decomposition using model coordination techniques are used to improve the performance (11). An example of a type of system optimization is illustrated in Figure 1.4. Geometric and topologic optimization creates geometric design variables that are, for instance, the coordinates of nodes in a finite element mesh or the pier location for a continuous bridge. In certain cases where the areas of the elements have zero as lower bounds, the unnecessary elements can be eliminated by the optimization algorithm. Sometimes the concept of separate design spaces, one for joint coordinates and the other for cross sectional element sizes, is used when trying to reduce the size of the design space considered at any stage (12). An example of optimal configuration is presented in Figure 1.5. In large optimization problems it is usual to use multilevel optimization techniques where the structural designer has to coordinate and optimize at several levels of the design process. This technique is also useful when the main goal is to find the optimum geometry besides optimizing Optimization of a Two span prestressed beam XI X2 XI to X6 Section geometry X7 to X9 Eccentricities of draped cable X1O Prestressing force Figure 1.4. System optimization. OPTIMIZATION OF TRUSS GEOMETRY Initial Configuration S. S. S. S. S. S. S. S. S. Optimal Configuration Figure 1.5. Geometry optimization. Load Load the structural elements. Design variables that control the geometry are often handled better when considered separately from the set of sizing variables (13). Study Objectives The main objectives of the present work are to combine adequately optimization and reliability concepts and to test the performance of the integrated approach to reinforced concrete frames. Reliability requirements are imposed at the element and the system level. At element level a maximum probability of failure is imposed for each element and at the system level a minimum reliability index is imposed for the failure mechanisms. The material behavior of the reinforced concrete elements is separated in two phases. The first considers linear material behavior and the second includes the concrete and steel nonlinear behavior. Structural frame optimization problems have usually been formulated based on the cycling between two distinct phases, analysis and optimal design. This methodology is the classical approach in structural optimization. The first phase consists in an initial sizing or structure definition. In the second phase, a structural analysis is performed and in the third phase, the structure is resized or redefined using Mathematical Programming or Optimality Criteria methods. The cycling between phases two and three is interrupted when the termination criteria are met. The research option summarized here combines phases two and three into one only stage. This is accomplished by the addition of the global displacements to the set of design variables. This addition implies that the equilibrium equations, solved explicitly in the cycling approach, are added to the set of constraints as equality constraints. These new equality constraints will be solved iteratively while in the cycling approach the solution is obtained using a Gauss type decomposition. The main objective behind the adoption of this strategy was to experiment this formulation where the variables related with element stiffness definition and the displacement variables are in the same design space. For that reason the simultaneous optimization and iterative solution of equilibrium equations could be more efficient than the classical nested approach. The application of this formulation was initially performed with elastic linear frames subjected to static loading. The constraints consisted of limiting the global displacements and the element stresses, besides the additional set of qualities representing the equilibrium constraints. The optimization method used consisted of unconstrained minimization of an augmented lagrangian function of the initial objective function and the equality and inequality constraints (14). Summary Results obtained with the integrated approach were encouraging and proved that the method was acceptable for elastic design purposes with displacement and stress constraints. Despite the fact that optimum values were obtained there was however an increase in the size of the problem. This modification of the problem size was due to the fact that the number of variables and the set of constraints augmented. The final type of optimization problem considered in this work was the minimization of the total cost of a reinforced concrete plane frame submitted to static loading considering the actual stress-strain diagram for concrete and the elastic-plastic behavior of the reinforcing steel. A typical element had a constant rectangular section and doubly reinforced with equal amount of flexural steel on both sides. Width and height of the cross sections had prescribed lower bounds, representing practical requirements and an adequate ratio between the height and the width. The amount of steel was limited by lower and upper bounds dictated by the minimum and maximum reinforcing steel ratios requested by the Building Code Requirements for Reinforced Concrete, commonly known as ACI 318-83. Inequality constraints considered included maximum values for the global displacements and a minimum reliability index for the element flexural failure function. Displacements allowed were based on serviceability requirements like cracking and relative story drift. The reliability indices were based on usual values of probability of failure used in design codes. Only the flexural behavior of the frames was analyzed since it is the most important for usual structures and the members were modeled as beam elements. Inelastic behavior of the structure due to the material nonlinearities imposes a change of the global stiffness terms independently of those dictated by the alterations of the dimensions during the optimization search. For that reason, the reinforced concrete element was modeled as a linear elastic beam with nonlinear rotational springs at each end. Rotational spring stiffness was considered infinite when the moment was below the yielding moment. Above that value the element stiffness was inverted to its flexibility and the inverse of the secant spring stiffness value was added to the corresponding diagonal terms. Spring stiffness was calculated using the secant value of the bilinear moment-rotation diagram corresponding to the current global rotation. Values of the yielding and ultimate moments were obtained by integrating the actual stress-strain diagram for the compressive force in the concrete. The corresponding rotation at a hinge was calculated by integrating the curvature diagram along the element. Element reliability was evaluated using a Level 2 method, i.e., an approximation to the evaluation of the exact probability of failure. The statistical variables considered were those assumed to have greater influence on the final result. These were the compressive strength of concrete and the external loads, assumed as normal distributed variables. The corresponding reliability index was calculated for constraint evaluation using the ultimate moment obtained from the integration of the respective strain diagram. Optimization techniques tested were based on the Augmented Lagrangian and the Generalized Reduced Gradient methods. The optimization problem was run, and after termination, the structure probability of failure was compared with the assigned value. If the result was not satisfactory, the process was restarted with updated values of the element reliability indices for the members involved in the most probable collapse mechanism. Evaluation of the system reliability was divided in two phases. First phase consisted of the identification of the elementary collapse mechanisms. In the second phase these elementary mechanisms were linearly combined to generate all significant mechanisms. System reliability was calculated considering the frame as a series system where each element is one of these mechanisms with higher probability of failure. Generation of the fundamental collapse mechanisms was made using Watwood's method (15). The automatic procedure consisted of using the geometric configuration of the frame and external loading to find all the.one degree of freedom failure mechanisms. The reliabilities of these mechanisms was calculated using the corresponding failure functions System reliability was evaluated using the Beta unzipping method (16). The elementary mechanisms were linearly combined to obtain other failure mechanisms. The corresponding failure functions were created and the associated reliability indices calculated. In each set of combinations only those in the closeness of the minimum were considered for the next combination (17). CHAPTER 2 INTEGRATED OPTIMIZATION OF LINEAR FRAMES Original Research Integrated formulations for structural optimization problems has received little attention in the published literature. The works of L. Schmit and R. L. Fox are considered the pioneering work as applied to integrated structural optimization (18). The concept of this structural synthesis problem is to combine the design variables with the behavioral variables. The immediate consequences of this concept are that the problem has a larger number of design variables and the traditional nested analysis-optimization process is avoided. This approach has not been popular since past performance was not comparable to the iterative techniques based on Optimality Criteria and Mathematical Programming concepts. In the integrated formulation the equilibrium constraints generate a large additional number of equality constraints. Several researchers have recently adopted the integrated approach with encouraging results. These recent attempts have been motivated by new solution procedures 21 considered more adequate for this type of formulation and by computer hardware development. An example is the optimization of elements with stiffness and strength properties proportional to the transverse size of the elements with linearization of the displacement constraints (19). Another algorithm uses the incremental load approach and conjugate gradient methods to optimize a structure subjected to nonlinear collapse constraints (20). In this case the stiffness matrix is approximated using the element- by-element technique (21). A more recent work uses a new solution technique based on Geometric Programming theory (22). In this formulation the equilibrium constraints are the sum of geometric terms that are function of the design variables. This chapter describes research that was conducted to study the integrated analysis approach for portal frames with linear behavior and static loading (23-26). The initial phase addressed only constraints on the displacements. Stress constraints were added on a second phase. Throughout this part of the study the frame elements had continuous prismatic rectangular cross section. Augmented Lagrangian Function The optimization technique of cycling unconstrained minimizations of a penalty function, based on an pseudo- objective augmented lagrangian function, was chosen as the solution scheme (27). The design variables were the areas and inertias of each element and the global displacements. Since it is a planar frame there are three degrees of freedom for each joint in the structure. The merit function used was the volume of the structure. In frames made with one material, volume is generally considered to be proportional to the structure cost. This value was calculated as the sum, for all elements, of the product of the element area times the respective length. The set of inequality constraints was generated by the structure physical behavior and material properties. Limits were imposed on the global displacements and, in the final stage, the element stresses were also bounded. The compatibility and equilibrium requirements were guaranteed by the additional group of equality constraints. This set was given by the product of the stiffness matrix and the vector of global displacements from which the vector of external global loads was subtracted. A brief description of the problem variables and respective formulation for a typical planar frame is the following: Structural parameters n structural elements; m number of global degrees of freedom; R vector of static external loads; D vector of bounds of m; Design variables Xk, k=l,3,...,2n-1 --- area of element (k+l)/2; xj, j=2,4,...,2n --- inertia of element j/2; xi, i=2n+l,2n+2,...,2n+m --- global displacements; Objective Function f(K) = Z lpxkk p=l,n where lp length of element p; Equality Constraints H(x) = K x* R where K global stiffness matrix; x* displacement vector; Inequality Constraints G(x) = x* D < 0 Augmented Lagrangian Function L(x,u,v) = f(x) + u H + P H H + v G' + P G'G' where u, V lagrangian multipliers; P penalty factor; G' maximum of (G, -y/2P}. The optimization procedure consists of several cycles of unconstrained minimization of the pseudo-objective function. The values of the lagrangian multipliers are kept constant during each cycle of the unconstrained minimization. At the end of an unconstrained minimization cycle, the multipliers are updated using an appropriate rule (12). The procedure is repeated for successive cycles until there is no significant change of the objective function. At this point the primal and dual optima have been found and the algorithm stops. Unconstrained Minimization Techniques Initially the technique used for the unconstrained minimization of the augmented lagrangian function was a zero order method referred to as the Hooke and Jeeves method or Pattern Search. The classification as a zero order method means that it does not utilize any information about the form or shape of the function. After the phase when stress constraints were added, a first order method, Steepest Descent, was tested as an improvement in the algorithm's performance (27). The technique is based on the gradient of the function that indicates the direction with the highest slope at a given point. Second order methods were determined inappropriate because the pseudo inequality constraints, g', have discontinuous second derivatives. Hooke and Jeeves method is an iterative procedure where each step may involve two kinds of moves. The first type of moves explores the local configuration of the pseudo- objective function along the directions of the design variables. The investigation is done within a prescribed step size from the current temporary design point. Each variable is investigated one at a time. The value of the step size is increased or decreased with success or failure in the exploration. This search along the coordinate directions will eventually lead to a smaller value of the pseudo-objective function. Otherwise the optimum has been reached and the exploration stops. Once all variables have been searched, a pattern move is attempted. The pattern direction is defined by the starting and final points of the variable search and a move is made along that direction. The process of exploration and pattern moves is repeated until there is no significant improvement of the pseudo-objective function. A graphic example is presented in Figure 2.1. The initial point of the variable search, 1, and the final point of that cycle, 4, define a pattern direction that yields a better design point, 5. HOOKE and JEEVES I Initial Point 6 Final Point X24 4/5 Pattern Move Function Contours Figure 2.1. Pattern Search. A computer program was written in accordance with the previous statements and discussions. The structure of the program was conceived by taking into account future inclusions of other types of constraints, changes in the minimization techniques, element replacements and extension to nonlinear and dynamic problems. Hence the program was divided into separate subprograms for the independent tasks (26). Final Results The performance and accuracy of the formulation described above was evaluated. Test examples for that purpose were structures with an explicit optimal configuration or simple frames. In the isostatic examples the optimal explicit solutions could be obtained and compared to the computer results. For the other structures, several runs were made with different initial design points and the optimal configuration was determined. Minimum values were imposed for the dimensions of the cross sections, represented by lower bounds of the areas and moments of inertia. The optimization results show the final values of the displacement variables as the exact solutions for the equilibrium equations. The final area and moment of inertia are also the expected optimal values. Results of a cantilever beam are presented in Figure 2.2. XI.X2 X :1 area of beam :2 inertia of beam 3 horizontal tip displacement 4 vertical tip displacement :5 tip rotation VARIABLE INITIAL FINAL Xl (in2) 1.0 6.55 X2 (in4) 1.0 78625 X3 (in) 0.4 0.500 X4 (in) 0.4 0.353 X5 (rad) 0.4 0.006 Figure 2.2. Cantilever beam. / Penalty factors used in these runs were of an order of magnitude greater than that of the objective function and constraints. They were kept constant during each optimization cycle. Scaling was also mandatory since the various terms of augmented lagrangian function have different orders of magnitude. The adopted scaling method consisted of using the inverse of the initial value of the expressions concerned. Initial guesses for the design variables were also important for the algorithm performance. The closer these initial designs were to the optimum, the faster the convergence rate. An updated version of this algorithm was created with the addition of stress constraints. The results of the structures used to test this addition illustrated the adequacy of the method for this type of problems. Again, for the cantilever beam with the explicit solution, the optimum results were obtained. For the frame, the final answer corresponded to what was expected and convergence was obtained. Final mass distribution resembles that previously attained just with displacement constraints. The geometry and related values are presented in Figure 2.3. A tapered cantilever loaded at the tip was compared with the results obtained using a recursive relation between the dimensions and displacements (12). The two sets of results, those from the reference and those from the program run, are very close. The maximum absolute difference 10 Kin 1OO00K 10 Kin ^1 ELEMENT INITIAL FINAL Area (in2) 1.0 25.4 Inertia (in4) 1.0 120224 Area (in2) 1.0 179 2 Inerti. (in4) 1.0 5912 Area (in2) 1.0 35.1 Inertia (in4) 1.0 17058 Figure 2.3. One bay frame. between the correspondent section dimensions is less than five percent. Further Improvements In subsequent developments, some other improvements were added to the algorithm that used the augmented lagrangian formulation. The first consisted of eliminating from the search those constraints that were inactive. Those constraints whose value did not show a change when the line search was along one of the design variable, were skipped from recalculation. This savings in computational effort allowed a reduction of forty per cent of the total run-time. This feature was discarded when the gradient search method was implemented. With this technique the changes in the design variables were done simultaneously, all constraints were altered and selective recalculation was no longer possible. Another significant improvement was achieved by starting the solution with feasible displacements. The displacement variables were calculated at the beginning of the program corresponding to the initial loading and frame dimensions. This led to the situation where the equality constraints were exactly satisfied at the start of the iteration procedure. This addition was kept in the version using the gradient search. Work was also done on selecting the initial cross section dimensions. Rules of thumb were found to expedite calculations to obtain acceptable initial values. The method of steepest descent makes use of the gradient of the pseudo-objective function. The gradient vector represents the line along which there is the highest variation of the pseudo-objective function at the actual design point. Moving in the direction defined by the negative of the gradient vector is expected to decrease the value of the pseudo-objective function. This direction is called the steepest descent. A graphical representation of the method is displayed in Figure 2.4. Since the explicit formulation of the gradient of the pseudo-objective function was not practical to obtain, the gradient vector was obtained using a finite difference technique. To obtain the minimum point along the gradient direction another design point along that line is found such that it has a higher pseudo-objective function value. Then, the optimum value should lie in this interval and a line search is performed using the golden section method. The gradient vector was normalized to avoid numerical ill-conditioning. For the same reason, constraints and the design variables were also scaled. Numerical difficulties are predictable if just one of the constraint function, or the objective function, is of different magnitude than the rest of the terms or its rate of change is considerably different from the others. Scaling factors for each constraint were evaluated as the ratio between the gradient STEEPEST DESCENT 1-Initial Point 4-Final Point X2 Function Contours Figure 2.4. Gradient method. of the objective function and the gradient of that constraint. Scaling of the design variables was also tried. The normalization of the design variables consisted of applying scaling factors that reduced them to a single order of magnitude. The results obtained with this unconstrained minimization technique were inferior to those using the Hooke and Jeeves method. The apparent reason was the shape of the surface generated by the augmented lagrangian function. Around a relative local optimal point, where the equality constraints are satisfied, the variation of the augmented lagrangian function was very abrupt. Consequently, any line search performed starting at a relative optimal point would invariably return to the same initial point. When using a set of design variables that was not a relative local optimum, the gradient search would still not converge. The reason for this lack of convergence was the numerical error created by the steep slope of the function. This fact could not be avoided despite the several combinations of the constraint and variable scalings aimed at smoothing the shape of the augmented lagrangian function. Another phase of research consisted in using a mixed method for the search. In a first phase, Hooke and Jeeves was used to obtain a better second point than the starting design point. This second point was then used to apply the gradient search. The procedure was repeated with the 36 consequent updates of the lagrangian multipliers. This mixed method did not present any improvement over the Hooke and Jeeves method. The important conclusion from the results of this mixed strategy was that convergence could only be obtained when enough iterations of the Hooke and Jeeves phase were completed. Consequently, the adopted unconstrained minimization method for the optimization of the augmented lagrangian function in the linear static formulation was the Hooke and Jeeves method. CHAPTER 3 NONLINEAR REINFORCED CONCRETE ELEMENT Introduction Reinforced concrete elements are made of two different materials, concrete and steel. Concrete is the massive component, has a high compressive strength and fails easily when submitted to tension. Steel is embedded whenever tensile strength is required. For that reason the additional steel bars are commonly designated as reinforcing steel. Adequate combination of these two materials originates a symbiotic composite material that has been widely used (28). These elements are designed with bending, compression and torsion requirements for code and safety compliance. In some cases tension is also allowed. Concrete and steel have nonlinear stress-strain diagrams. Consequently, when material nonlinearities are included, modeling of the behavior of any composite element is very difficult (29-30). When loads produce a tensile stress greater than the maximum allowable value for the concrete cracking results. When reinforcing steel stress 37 reaches the yielding value there is a large strain and section curvature increase. Geometric nonlinearities are then created by extra rotations of flexural elements from the cracking and steel yielding. A basic assumption in nonlinear analysis of reinforced concrete frames is that the element rotations with relation to the line defined by the nodes, chord rotations, are small and the theory for straight elements may be applied with some adaptations. The most popular analysis techniques are based on incremental loadings of the structure and are known by the initial stiffness and tangent stiffness methods. A technique based on the application of the entire load at a single step is known by the secant stiffness method. This last technique was chosen for the analysis of the structure since it is more adequate to the optimization formulation. Element Modeling Survey In the last three decades there have been many attempts to create a simplified beam model of the inelastic reinforced concrete element (31-33). The main objective for this research has been to advance a solution providing precise results within reasonable computational and memory storage limits. The study has a significant importance for the analysis of reinforced concrete structures submitted to dynamic loads (34-35). In these examples the moments at the ends are close to the ultimate allowable values. This closeness implies that the concrete and steel stresses are in the nonlinear intervals of the stress-strain diagrams. The frame behaves as if inelastic plastic hinges have formed due to concrete cracking and steel yielding. Initial studies in this area addressed simple structures with moment-rotation relationship conditioned by the moments at the beam extremities. This produced the one- component model with nonlinear rotational springs at the ends. Later, another theory assumed a bilinear moment resistance with two parallel elements, one to simulate yielding and the other to represent strain hardening. Several variations of these two theories have been developed and experimentally tested (36). Recent improvements in computer software led to sophisticated modeling of reinforced concrete elements using nontraditional finite element techniques. A simple approach to this type of problem is based on the theory of damage mechanics (37). The beam element is modeled as a macroelement divided in models with explicit and accurate behavior. The behavior of the whole structure is then extrapolated from the small elements. These types of models have been tested thoroughly to ascertain its reliability and accuracy (38). These evaluations, made mostly by comparison of computer program results with experimental test data, provided a great deal of information for further enhancements and refinements. The option for this study had to fall on a element model that is a compromise between the accuracy required and the cyclic nature of the optimization process (39). Repeated evaluation of the element stiffness due to the changes of the physical properties of the elements is required. For this reason it is highly desirable to choose a model with low computational requirements. Beam Element with Inelastic Hinges Given the available solutions for the model of the reinforced concrete element, the one-component model was chosen as shown in Figure 3.1. It is a simple idealization that doesn't increase the total number of elements of the structure. This model has shown to accurately model the nonlinear behavior of reinforced concrete, even for dynamic loadings (40). Some basic assumptions and simplifications were made for the definition of the model. For example, the fact that concrete cracks under tensile loading, causing local nonlinear behavior, was not accounted for. Time dependent properties of the concrete were not considered. Shear effects were not included in this formulation. The loads were considered applied at the nodes and elements with loads in the span can be approximated by a discrete number of elements with nodal loads. The unique element internal action considered was flexure. Yielding of the reinforcing steel may only take place in the hinges at the element ends. Strain hardening One-Component Model Reinforced Concrete Element Linear Elastic Element i/G)~v I :_Q1y with Secant Stiffness Figure 3.1. Element model. Sprin and related altered element stiffness are simulated by the linear element with nonlinear rotational springs at the extremities. Inelastic rotations of reinforced concrete hinges at the element ends are determined as a function of the respective moment-curvature relationship for each element. These curves are redefined every time any element sectional properties changes during the optimization process since the ultimate and yielding moments also change. A typical moment-curvature diagram for reinforced concrete elements is bilinear. It is obtained assuming material stress-strain curves that are parabolic-linear for the concrete and bilinear for the reinforcing steel as shown in Figure 3.2 (28). The stress in the concrete is designated by fc and the stress in the steel reinforcement is represented by fs. The algorithm used to compute the moment corresponding to a certain strain diagram is an iterative Newton based iteration that determines the depth of the neutral axis guaranteeing equilibrium of the internal forces. Then, after determining the internal coupled forces the related moment is computed. All reinforced concrete elements are doubly reinforced with equal areas of steel on both sides. This assumption is valid for columns and acceptable for beams since in continuous frames there are moments of different sign along the beams. Evaluation of the moments for each reinforced concrete section was based on the exact internal equilibrium equations as follows: B f C c C A _ 0.002 0.004 Concrete Stress-Strain Diagram fs 'S Steel Stress-Strain Diagram Figure 3.2. Material behavior. Cc + Cs = Ts where Cc compressive force in the concrete and is equal to the area under stress-strain curve corresponding to concrete strain Ec; Cs compressive force in the steel area As corresponding to steel strain Ecs; Ts tensile force in the steel area As corresponding to steel strain Es (ES y Ey). Typical element moments necessary to define the bilinear moment-curvature diagram were the yielding and ultimate values. These characteristic values were determined considering the corresponding section strain distribution, the stress-strain diagrams for concrete and steel, the location of neutral axis and the moment of the internal forces as shown in Figure 3.3. The compressive force of the concrete is given at any time by Cc = a fcm b kd where Ieca a = fc/(fcmEca)dec; 0 fcm maximum flexural concrete stress; Eca concrete strain at the top compression fiber; b element cross section base; kd distance of neutral axis from top compressed fiber. SECTION CHARACTERISTICS AS -v_1^ h' As EC yCC yp-8 Geometry Strain Diagram Forces Figure 3.3. Reinforced concrete section. The force in the compressed steel is given by Cs = As fcs where As steel area; fcs stress in compressed steel. The force in the steel under tension is determined by Ts = As fy where fy yielding steel stress. For instance, the internal ultimate moment is given by the moments of these three internal forces about the top compressed fiber. For that reason a parameter 0, that defines the centroid of the concrete compressive stress diagram, is introduced as *Eca Ieca S=1 ec fc dec /(Eca fc dEc) 0 0 These parameters, a and n, when the ultimate concrete strain is defined as ec = 0.004, become a = 2/3 (region AB) n = 3/8 (region AB) a = 0.9 (region BC) n = 0.51851 (region BC) where the regions AB and BC are defined in Figure 3.2. The section flexural strength, Mi, may be defined as Mi = Cs d'+ Cc 0 kd Ts d where d'- distance of Cs to top compressed fiber; d distance of Ts to top compressed fiber. Element curvatures corresponding to these yielding and ultimate moments are obtained assuming that plane sections remain plane after deformation and there is no strain hardening of the reinforcing steel. These formulas are as follows: #y = (Ey + Eca)/d u = (Esa + Ecu)/d where -y yielding curvature; Ey Es / fy; Eg 29x106 psi; fy yielding stress of reinforcing steel; Eca maximum concrete compressive strain; u ultimate curvature; Esa actual tensile strain of steel; Ecu ultimate compressive strain of concrete. These section characteristics define section diagrams as shown on Figure 3.4. The value of the ultimate rotation was given by the integration of the curvature along the element. Two types of curvature diagrams were considered Moment Curvature Diagram Linearized Diagram Moy u 0D Moment Rotation Diagram Figure 3.4. Element deformation diagrams. 1 u , for integration. The first one was when moments at element ends had the same rotational direction and the second when the rotational directions were opposite. In both cases a simplified method was used to integrate the curvature along the element to find the corresponding rotation since the moments at the other end were kept constant. Yielding rotation for any node of the element was calculated assuming the yielding moment at that node and keeping the other moment unchanged. The same method was applied for the calculation of the ultimate rotation where a modified curvature diagram was used as schematically exemplified in Figure 3.5. Beam Element Stiffness The elastic element chosen has a stiffness derived in classical terms. End rotational springs had variable stiffness depending on element moments at the nodes. A large value was assigned to the secant spring stiffness when moments were below the yielding value assured a linear behavior. The secant stiffness value obtained from the moment rotation diagram was used for moment values above yielding. The strain hardening ratio of the linearized moment rotation diagram was computed as the difference between ultimate and yielding moments divided by the difference between the ultimate and yielding rotations. A MOMENT DIAGRAM Mi M Mi Mj j Mi Moment at node i Mj Moment at node j My Yielding moment CURVATURE DIAGRAM ] (u Ultimate curvature y Yielding curvature (i- Curvature at node j Figure 3.5. Curvature integration. graphical description of these definitions is presented in Figure 3.6. The element modified stiffness was derived from the condensation of elastic stiffness matrices of the linear elastic element and the rotational spring elements. To condense the two matrices the first step consisted of inverting the sum of the corresponding flexibility matrices concerning the independent element rotational degrees of freedom. The next step was the expansion of this element stiffness to include the axial displacements, uncoupled from the spring rotations, and the other dependent element degrees of freedom. The main steps of this step are the following: 1/Ksi 0 -1 -1 0 1/3 -1/6 + 3EI/L 1/Ksj -1/6 1/3 -1 0 0 1 0 0 [ a ] = 0 1/L 1 0 -1/L 0 0 1/L 0 0 -1/L 1 [ Kmod ] = [ a ]t [ Ks* ] [ a ] Ks secant stiffness matrix with element rotations; Ks* expanded secant stiffness matrix with uncoupled axial stiffness; Ksi stiffness of spring at node i; Ksj stiffness of spring at node j; [ K ] where Moment Mu My Rotation Spring Moment-Rotation Diagram Mu Ultimate moment Yielding moment 10e30 (Mu My)/(u - Oy) Ksec Spring stiffness for M > My Figure 3.6. Secant spring stiffness. My - Kl - K2 - E element modulus of elasticity; I element moment of inertia; L element length; a expansion matrix; Kmod modified element matrix. After evaluating the modified element stiffness matrix it was transformed from the local coordinates to the global coordinates by the use of the corresponding rotation matrix. The values of the terms of this element stiffness matrix were then used to compute the corresponding updated equality constraint values. The process was similar to assembling a structure global matrix using a location matrix relating the element degrees of freedom with the structure global degrees of freedom. CHAPTER 4 STRUCTURAL ELEMENT RELIABILITY Introduction Design and checking of structures in the field of Civil Engineering has been traditionally based on deterministic analysis. Adequate dimensions, material properties and loads are assumed and an analysis is carried out to obtain the required evaluation. Nevertheless, variations of all these parameters and questions related to the structural model may impose a different behavior than expected (41). It must be emphasized that if there were no uncertainties related to the prediction of loads, materials and structure modeling, then the respective safety would be more easily guaranteed. For these reasons the use of probabilistic principles and methodologies in structural design has been increasing. Design for safety and performance should consider the conflict between safety and risk. The objective of probability concepts and methods is to develop a framework where the effects of these uncertainties are considered. Structural reliability has received the attention of several 54 researchers and, consequently, it is introduced into almost all recent structural codes worldwide. It is a relatively young structural science that evolved in the same way as other new areas where theoretical studies dictate the general principles for systematic treatment of problems. There are however practical difficulties in obtaining enough statistical data and handling the sophistication of the probabilistic methods. For these reasons the analytical processes involved in the determination of structural reliability were grouped in different working levels (42). These working levels depend on the problem considered and the desired accuracy for the reliability evaluation. There are three basic levels and the classification increases with the sophistication of the method used and the amount of statistical data that is manipulated. Level 1 uses a methodology that provides a structural member with an adequate structural reliability by the specification of partial safety factors and characteristic values of design variables. This is the method currently used in structural design codes (43). Level 2 includes all methods that control the probability of failure at certain points on the failure boundary defined by a limit state equation (44). Level 3 groups all techniques that perform a complete and exact analysis of the structure taking into account the joint probability function of all the variables involved (45). In this chapter, the technique used to analyze the structural reliability of each reinforced concrete beam element is described. Due to the nature of the problem, where optimization and reliability evaluation are performed simultaneously at the element level, a Level 2 method was chosen. Since the concepts of limit state design and probability of failure are intimately connected with structural reliability, a brief description is also included. Concept of limit state may be described as that state beyond which a structure, or part of it, can no longer fulfill the functions or satisfy the conditions for which it was designed. Namely, the structure is said to reach a limit state when a specific response parameter attains a threshold value. Examples of ultimate limit states are the loss of equilibrium of a part or the whole of the structure considered as a rigid body, failure or excessive plasticity of critical sections due to static actions, transformation of the structure into a mechanism, buckling due to elastic or plastic instability, fatigue, excessive deflections and abundant cracking. Modern codes divide limit states into two main groups. Ultimate limit states, corresponding to the maximum load- carrying capacity, and serviceability limit states, related to the criteria governing normal use and durability (46). For each of these groups the importance of damage is different and is represented by the adopted respective probability of failure. For instance, in reinforced and prestressed concrete, code checks for the ultimate limit states are based on element forces, except in the plastic analysis where the design variables are the loads. In cases where fatigue is involved, stresses are also the control variables. The service limit states are the cracking limit state and the deformation limit state. In this work only the ultimate flexural limit state and the global deformation limit state are addressed since they are the more relevant for the optimization study. Acceptable risks of failure for any structure are affected by the nature of the structure itself and its expected application. These are dependent on social and local variations. It is common for structural engineers to balance the contradiction between the economy and safety of the structure. This particular aspect is the main reason why it is so appealing to combine reliability and optimization in structural design. Probabilities of failure used in limit state designs vary with the risk of loss of human lives, the number of lives affected and economic consequences. In ultimate limit states the range of probability of failure adopted is between 10-4 and 10-7 over a 50 year expected design life. In serviceability limit states the probability of failure varies between 10-1 and 10-3. A criterion proposed is as follows (41): pf = 10-5 U T / L where U 0.005 ........ Places of public assembly, dams; 0.05 ......... Domestic, office, industry, travel; 0.5 .......... Bridges; 5 ............ Towers, masts, offshore structures; T life period of the structure(years); L number of people involved. These values must be interpreted carefully. For example, the value of 10-3 means theoretically that, on the average, out of 1000 nominally identical buildings, one will crack or deform excessively. It is evident that in civil engineering 1000 identical buildings rarely occur, even neglecting the fact that a statistically significant number require samples at least 10 to 20 times larger. Moreover, the determination of these low probabilities requires extrapolations of statistical properties that are experimentally known only around the mean values of the random quantities. For these reasons, the probabilities of failure in civil engineering have no real statistical significance and they must be considered not as deterministic quantities but just as conventional comparative values. In consequence of the above considerations, the differences between the methods used in each of the three levels are rather operational than conceptual. There are no rigid boundaries between them. They are used in accordance with the required accuracy and the nature of the problem to be studied. Level 3 methods require a complete analysis of the problem and also the integration of the joint distribution density of the random variables extended over the safety domain. They remain in the field of research and are used to check the validity of approximations, idealizations and simplifications performed in the other two levels. Level 2 methods use random variables characterized by their known or assumed distribution functions, defined in terms of important parameters as means and variances. This avoids the multidimensional integration of the previous method. These methods may be used by engineers to solve problems of special technical and economical importance. Code committees engaged in drafting and revising standard codes of practice use them to evaluate the partial safety factors. It is possible that computational developments in the near future will allow for such methods to be more commonly used by the practicing engineer. The probabilistic aspect of the problem in the Level 1 methods is represented by characteristic values of the random variables involved. With these characteristic values partial safety factors are derived using Level 2 methods. They are used by most engineers where reliability theory and probabilistic methods are the basis of their code provisions. These Level 1 methods could be replaced by the Level 2 methods if an agreement was obtained in the following issues: selection of basic random variables for each specific problem, their distribution types and relative statistical parameters; form of the various limit state equations and choice of models; operational reliability levels to be adopted in different design situations. It must be emphasized that the advantage of Level 1 schemes over Level 2 are their great operational simplicity due to the use of fixed and constant partial safety factors for a given class of design situations. The main disadvantage of Level 1 is the selection of partial safety factors for a given structural class in such a way that the efficiency of the method proposed is satisfactory. It must assure that the deviation of the reliability of a design made on the basis of the adopted coefficients from the desired reliability level laid down in the code is acceptable. Two Dimensional Space Example Let R and S be two random variables, where R defines strength and S the load. Then the limit state function z shown in Figure 4.1 is defined as r z r- 0O ( z>0 ) ijil SAFE D' (z UNSAFE Safe and Unsafe Design Regions Figure 4.1. Design safety region. z = r s where r resistance function; s load function. The domain D (z>0) is the safe domain and D'(zO0) is the failure domain. The probability of failure, pf, is the probability that a point (R,S) belongs to D'. Once the statistical distributions of the random variables R and S are known, the numerical solution of the corresponding equation will determine pf. Assuming that both variables R and S have a Gaussian distribution, and further defining rm and sm as the mean values, and aR and aS as standard deviations of R and S, respectively, the random variable Z will also be normal and its statistical parameters are defined as zm = rm sm az = (a2R + a2S) where zm mean value function; az standard deviation function. Defining Fz as the cumulative normal distribution function, the probability of failure may be calculated as pf = P(Z A graphic representation of these functions is presented in Figure 4.2. Introducing the standardized variable u and the reliability index as u = (z zm) / aZ S= zm / Cz = (rm sm) / (a2R + a2S) then the probability of failure may be expressed as pf = Fu(-zm / az) = Fu(-B) An important concept widely used in structural safety when considering random variables is the Central Safety Factor. It relates the mean values and coefficients of variation of R and S to determine a probabilistic safety factor (44). It is a simplistic way of establishing some influence on the design variables of the respective random characteristics. To consider a more detailed study a Level 2 method is applied in the element reliability evaluation. In this method safety checks are made at a finite number of points of the failure boundary. A graphic representation in a two dimensional space is presented in Figure 4.3. In the case where this check is made at only one point, the parameter to be determined is the minimum distance between the origin of the system of the standardized variables to the boundary of the safety domain. Z mean Probability Density Function Cumulative Density Function Figure 4.2. Probabilistic functions. It is possible to associate this distance with a precise meaning in terms of reliability. A technique derived from this concept is the Lind-Hasofer Minimum Distance method illustrated in Figure 4.3 (47). Let X (X1, X2,..., Xn) be the vector of the basic random variables of a given structural problem that may be assumed to be statistically uncorrelated, involved in a given structural problem. Let z = g(xl, X2,...., Xn)= 0 be the boundary of the safety domain. The values of X belonging to the failure domain will satisfy the inequality z = g(x) < 0 The method consists in projecting the function z in the space of standardized variables defined as ui = (Xi xmi) / axi Measuring, in this space, the minimum distance B of the transformed surface g (ul, u2,....., un) from the origin of the axes. A design is regarded reliable if B > B*, where B* is prescribed by an appropriate code provision. In geometrical terms, the hypersphere having radius B* and with center at the origin of the axes ui is required to lie within the transformed safety domain. The justification for such a method is that most of the joint probability density of the variables involved will be concentrated in Ul Z-= 0 Z (UI,U2) > 0 Figure 4.3. Safety checks. the hypersphere having radius 8*, and that consequently it will be associated with values of vector X belonging to the safety domain. Mathematically, the problem to be solved is to find B = min (E u2i) In a great number of cases the safety boundary domain is linear, and one can write an expression for z as follows: z = g (xl, x2,...., Xn) = b + Z aixi Then, B can be immediately determined as follows g (ul, u2,..., un) = b + Z aixmi + Z aiaxiui = 0 and the distance of this hyperplane to the origin is B = Z (ai.xmi + b) / (Z a2iax2i) Expressing in terms of the standardized variables is equivalent to replacing the hypersurface by the hyperplane passing through P*, the point of minimum distance between the two geometric elements. A graphical illustration of this approximation in a two dimensional space is presented in Figure 4.4. Finally, the probability of failure, pf, and Z = 0 Z (U1,U2) < 0 >0 Figure 4.4. Reliability index. the reliability index, B, are within certain approximations related by Pf = 1 p(B) where p is the function of standardized cumulative normal distribution. Reinforced Concrete Element Reliability The element actions considered in analysis are only the moments at the member ends. These are the points of maximum value since only concentrated nodal loading is considered. The failure function z is then defined as z = r s = Mi Me where Mi ultimate internal resisting moment; Me maximum external element moment. The external moment at the section is obtained from the element displacements using the condensed element stiffness matrix defined in the previous chapter. The expressions to obtain the value of Mi were defined in the previous chapter. The random values chosen in this study were the characteristic strength of the concrete, f'c, and the maximum external moment in the element, Me. All other variables of the expression defining Mi could be taken as random but concrete strength was chosen due to the high coefficient of variation. Thus, the flexural failure function is linear and the respective reliability of failure can be easily calculated. Compressive strength of concrete is influenced by a large number of factors grouped in three main categories, namely materials, production and testing. Material variability depend on the cement quality, moisture content, mineral composition, physical properties and particle shape of aggregates. The production factors involve the type of watching, transportation procedure and workmanship. Testing includes sampling techniques, test methodology, specimen preparation and curing (48). It is difficult to evaluate correctly the importance of these three groups of factors. Their importance is certain to vary for different regions and construction projects. It has been found that the distribution of concrete compressive strengths can be approximated by the normal (Gaussian) distribution (49-50). Characteristic concrete compressive strengths obtained from a sampling of test data leads to a conclusion that for strength levels between 3,000 and 4,000 psi, the coefficient of variation is constant. For strengths beyond that range the standard deviation is constant (51). Since the values in reinforced concrete frames used are generally within the first interval the statistical value considered was the variance of f'c. The average standard variation for 68 a good quality control testing at the construction site is 550 psi. Using a 3500 psi specified compressive strength of concrete, f'c, the required average compressive strength of concrete is the larger of the following (51): f'cr = f'c + 1.34*0 = 3500 + 1.34 550 = 4237 psi or f'cr = f'c + 2.33*0 500 = 3500 + 2.33*550 500 = 4282 psi The coefficient of variation of f'c for this range of characteristic compressive strength is then given by V = a/f'c = 550/4282 = 0.128 and consequently the coefficient of variation of the concrete compressive flexural strength was adopted to be 0.15. External loads have different coefficients of variation for the different types of loads (52-53). For most design and construction in the United States a good estimate for the coefficient of variation of dead loads is 0.10. For the live loads the coefficient of variations are very high and range from 0.39 to 1.04. For that reason and since the building codes prescribe large values for live loads that exceed the mean value a single coefficient of variation of 0.15 was adopted for the combination of dead loads plus live loads. Basic variables considered, fc and Me, are assumed to have a probability function with normal distribution. This assumption is correct for the characteristic compressive strength of concrete but it does not hold for all external loads that create Me. In the case where a statistical refinement of the basic variable Me is required, there are techniques available to address the problem (47). Since flexural failure function, z, is linear the reliability index B of each element can be calculated for any given external moment, section and material properties. Denoting the basic variables fc as xl, and Me as x2, and eliminating the other parameters involved in the equation, the flexural failure function takes the form z = al xI + a2 x2 + b where al = a 0 b (kd)2; a2 = -1; b = As fcs d'- As fy d. Standardizing the variables xl and x2 leads to a replacement of the basic variables ul = (xI l)/a1 U2 = (x2 P2)/a2 where M1 mean value of fc; p2 mean value of Me; a1 standard deviation of fc; a2 standard deviation of Me. Replacing the standard normal variables in the flexural failure function the expression assumes the following form: z = alalul + a2a2u2 + alj1 + a292 + b Then the reliability index for each element is given by the distance from the standardized failure function to the origin of the standardized basic variables as follows: B = (al,1 + a242 + b) /(alal + a2a2) CHAPTER 5 SYSTEM RELIABILITY Introduction Optimum structural design techniques are mainly based on deterministic assumptions. There is no doubt that some of the design variables should be considered including their random nature (54-55). Of course system reliability problems are more complicated than element reliability problems. This is evident since it must consider all multiple element failure functions, the several failure modes and, in some cases, the correspondent statistical correlation. Another reason for including reliability considerations in structural optimization procedures is that, in some instances, the optimal solutions found have less redundancy and smaller ultimate load reserve than those solutions obtained with traditional design techniques (56-57). There is no doubt that the combination of optimum design techniques and reliability-based design procedures creates a powerful tool to obtain a practical optimized solution. The objective is to find a balanced design 74 between all those that satisfy the optimization constraints and at the same time will have the lowest allowable probability of failure (58). The strategy employed to evaluate the system reliability is described in the rest of the chapter. The elementary failure mechanisms of the structure are determined using Watwood's method. Then the system reliability is approximated using the Beta unzipping method, which consists of determining the relevant collapse mechanisms through linear combinations with fundamental mechanisms. The theory related with these techniques is tentatively described. System Reliability and Optimization A possible inclusion of the system probability of failure is to attribute a cost to system failure. This option originated a formulation based on the minimization of the total cost with the traditional optimization constraints (59). The objective function is as follows: Minimize Ct = Co + Cf Pf where Ct cost of the structure; Co initial cost of the structure; Cf cost of failure; Pf probability of structure failure. This option is not commonly used for inhabited structures since it is difficult to evaluate the economic value of a structural failure where human life losses are expected. A more popular alternative is to include an additional constraint representing the maximum probability of failure allowed for the structure (60). The constraint for the system reliability will be of the type Pf(X) < Pm where Pf probability of system failure; X vector of design variables; Pm allowable probability of system failure. When performing structural optimization one may consider serviceability and ultimate limit states. This possibility leads to another type of formulation where the objective function and constraints for these limit states are considered simultaneously (61). This type of problems are called reliability-based optimization and can be summarized as follows: Minimize Co subject to Gi(X) < 0, i=l,m Pu Puo Ps p Pso where Gi optimization constraints; m number of behavior constraints; Pu probability of ultimate system failure; Puo maximum probability of ultimate system failure; Ps probability of serviceability failure; Pso maximum probability of serviceability failure. The option adopted consisted of adding a constraint on the system failure. The value of the system failure at the end of the optimization cycle is compared with the target value. If it is not satisfactory the element requirements are modified and the optimization is restarted. Methods In determinate structures the collapse of any member will lead to system failure. The probability of system failure can be obtained as the probability of the union of member probability failures (16). These types of structural systems are called series systems or weakest-link systems. Redundant structures will fail only if all redundant members collapse. If this condition does not arise, whenever a member fails there will be a redistribution of loads in the system. These types of structures are called parallel systems. Graphic examples are presented in Figure 5.1. PARALLEL SYSTEM Load Truss Bars Load Figure 5.1. System models. Truss Bars SERIES SYSTEM ~1 Series systems with n elements have n failure modes. Parallel systems with n elements have more than n failure modes. These failure modes in parallel systems are dependent on whether the failure type of the elements is brittle or ductile (62-63). For redundant brittle systems the failure of an element and consequent redistribution of the loads will provoke the system failure. In these cases the system behavior may be considered to be generally identical to that of as a series system. Probability of failure of a series system can be considered as the union of the elements probability of failure Pfs = P(Ui(Zi: 0)|i=l,n) where U union of events; Pfs probability of system failure; Zi failure function of element i. If the element failure functions are not correlated then the evaluation of Pfs is relatively easy and may be assumed as Pfs = 1 Vi=1n(l P(ei=0)) where r product; ei = 0 if element is in a failure state, ei = 1 if element i is in a non-failure state; P(ei=0) probability of failure of element i. When there is correlation between element failure functions then the calculations become more complicated and time consuming. To avoid the exact evaluation, approximation and bound techniques have been developed (64- 65). The best known is the simple bounds. In this approach the upper bound for the probability of system failure assumes that all element failure functions are uncorrelated and the lower bound is obtained assuming full dependence between the element failure functions. If a more sophisticated bounding technique is necessary the Ditlevesen bounds may be used (17). The drawback is that this sophistication implies the calculation of event joint probabilities. A similar simplified approach to that used in series systems may be adopted to find the simple bounds for the failure of a parallel system. In the case of parallel systems the lower bound corresponds to the case where there is no dependence between the elements failure and the upper bound corresponds to full dependence between all elements failure (66). Exact evaluation of the probability of system failure is very difficult to obtain if the system has more than three elements. To solve a general problem, approximation or bounding techniques are used. For instance, for redundant ductile systems there is a large list of procedures, most of them with limited application (67). Some methods for redundant systems involve the determination of all collapse modes and their respective probability of failure. To obtain all collapse modes the fundamental mechanisms are determined and a Monte Carlo simulation is performed to generate all others. Afterwards the respective probabilities of failure are determined. This approach, although accurate, is very demanding in computational effort if the system is complex, and consequently is used mostly to validate the performance of other methods. In redundant ductile systems a variation of the Monte Carlo approach PNET or Point Estimate of System Collapse Probability is used. This consists in linearly combining the fundamental failure modes with the coefficients as variables. An objective function representing the reliability index of that combination is minimized and the most probable failure mechanisms are defined. Concerning redundant structures with brittle or ductile elements, other approximation and bounding techniques have been developed and studied based on graph theory. Two of those approaches for obtaining the probability of failure are the failure mode approach and the stable configuration approach (68). Both methods require the determination of all possible failure modes and the use of algorithms based on graph theory. To exemplify the determination of all possible failure modes the initial step is to build a directed network, or directed graph, with all possible events involving element failures that will lead to a collapse. Each node represents a stable configuration and each branch corresponds to a element failure. Each path is a set of branches connecting the initial and final nodes. A cut of the graph is a set of branches containing only one branch from every path. A simple example is presented in Figure 5.2. Methods based on the determination of fundamental failure mechanisms using practical simplifications from graph theory have been implemented (69). The Beta unzipping method and the branch and bound method are two examples. The principal advantages are that they are precise and easy to use. The Beta unzipping method finds the significant collapse mechanisms using combinations of fundamental mechanisms and rejecting those combinations that are outside a prescribed interval. The branch and bound method selects all failure paths that have high probabilities of occurrence. Although less exact, the Beta unzipping method was chosen due to its simplicity and performance. Generation of Failure Modes To define all failure mechanisms, the first step consists of determining the set through manipulation of elementary failure mechanisms. To obtain these, the method 2 3 1Tru Load ss Bars FAILURE GRAPH 2F F Bar Failure Figure 5.2. Failure graph. adopted was conceived by Watwood (15). It is an automatic tool to generate all failure mechanisms with one degree of freedom, or elementary failure mechanisms, of a given frame. The set of these mechanisms and all their linear combinations constitute all possible collapse configurations (70). The technique is relatively simple to use since the input data for this method is the same for traditional elastic analysis like joint and element information. Elementary failure mechanisms are dictated by the geometry of the structure and potential hinge locations created by the external load configuration. Hinge locations are considered at the end of each member. In the case where there are loads in the middle of the element, they are also considered at the points of concentrated or discretized loads. The element axial collapse is not considered in this formulation although it was included in the original version. Element global displacements of a planar frame form a vector with six variables, {S}. Using a cartesian referential set of axes x and y the displacements, S1 to Sg, may be represented as in Figure 5.3. Element deformation parameters may be defined by three independent quantities S'1 displacement about node i; S'2 rotation of node i; S'3 rotation of node j. S2 1 SI i c^\s p 4) SsI Element Displacements Independent Element Displacements Rigid Body Displacements Figure 5.3. Element displacements definition. a@ - 0 i * ^ *' - s5t SS4' When a mechanism is formed each element moves as a rigid body. The rigid body motion of an element of a planar frame can be defined by three parameters. They can be expressed in terms of the global coordinates x,y as S'4 translation in the x direction; S'5 translation in the y direction; S'6 rotation about node i. Two sets of three independent displacements, rigid body parameters and element deformations, create the transformed coordinate vector, (S'}. A relation can be established between local global coordinates and transformed coordinate vector represented by a linear transformation [T]. (S) = [T] {S'} where 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 1 [T] = 1 0 0 1 0 0 0 0 0 0 1 -L 0 0 1 0 0 1 L element length. For any elementary failure mechanism the element deformations, S'1, S'2, S'3, must be zero. This is only for elements that do not have plastic hinges. To materialize this condition, a matrix Ck is introduced for each element 87 k. This matrix is created with the first three rows of matrix T-1 for the kth element. The global condition matrix, C, is a block diagonal matrix consisting of the Ck matrices as follows: -1 Ck = 0 0 C1 C -- 0 O 0 0 -1/L -1/L 1 0 0 1/L 0 1/L Using the previous matrices and vectors relation now holds the following C (S) = (S'd) where (s) = - first element - second element - nth element and (S'd) = r S'1 S'2 S'3 S'1 S'2 S'3 - first element - nth element Compatibility between the element displacements, {S} and the structure global degrees of freedom {r} can be established {S) = [Q] [A] {r} where [Q] rotation matrix; [A] compatibility matrix. From previous equations the following expression holds [C] [Q] [A] {r} = {S'd} or [B] {r) = {S'd) An elementary mechanism of the structure is a solution of the homogeneous system [B] {r} = 0 If the structure configuration is not a mechanism there is no solution for the system except the trivial solution. To obtain a mechanism, releases of the global degrees of freedom must be introduced. Two releases per element are added corresponding to the hinges at the ends or points of application of concentrated or discretized loads. Each release corresponds to an addition of an external global degree of freedom. Addition of external degrees of freedom is done by replacing a row in matrix [Q] [A] with zeros. The changed rows correspond to the element degrees of freedom S3 and S6, the node rotations. For each row that is replaced, a unit column vector is added to the matrix [Q] [A] with a 1 in the row that has been replaced. The dimensionality of (r} is increased by the number of rows replaced in [Q] [A]. The total is a set of extra columns with a dimension that is twice the number of elements. The homogeneous system becomes [C] ([Q] [A])* {ra) = [B'] {ra) = 0 where ([Q] [A])* matrix with extra 2n columns; ({a) vector of increased global degrees of freedom. Matrix [B'] is not square and has a greater number of columns than the number of rows. The solution of the system of homogeneous equations above may be obtained using a technique similar to that when solving for redundant unknowns in the force method. Difference between number of rows and number of columns is the number of independent solutions, that coincides with the number of elementary mechanisms. Suppose the rank of [B'] is the number of columns, m, and that the number of columns is p. In this case one can find a matrix [D], nonsingular with dimensions p by p such that [B'] [D] = [I I 0] where [I] identity matrix, m by m; [0] null matrix, m by (p-m). Last columns of [D] are independent solutions of the homogeneous system of equations since they are orthogonal to the rows of [B']. To obtain [D], a reduction is performed on the columns of [B'] that is conceptually identical to a Gauss-Jordan reduction (15). The solution of such a system of equations is illustrated in Figure 5.4, where all elementary failure mechanisms for a two story frame are presented. Beta Unzipping Method Advantages of the Beta unzipping method, as stated before, are important. It can be used for reliability estimation of planar and spatial trusses and frames made with ductile or brittle elements. The probability of failure can be evaluated with different levels of accuracy. It is also a method that can be easily implemented for automated calculations. Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 5 Mode 7 Mode 8 Figure 5.4. System failure modes. |

Full Text |

172
C*************************************************************** ndof=iqh do 20 j=l,ndof P(j)=r(j) 20 continue C** ****************** ******c***** ********************* ************ c ordering theta and r matrices c*************************************************************** do 710 k=l,numec-l jflag=0 do 720 i=l,ndof if(abs(rb(i, k)).gt.0.0)j flag=i 720 continue if (jflag.eq.0)then do 730 l=k+l,numec do 740 li=l,ndof if(abs(rb(li,l)).gt.0.0)then do 750 lj=l,ndof temp(lj,1)=rb(lj,1) rb(lj ,l)=rb(lj, jflag) rb (1 j j flag) =temp (1 j 1) 750 continue do 760 lj=l,2*nel temp(1j,1)=theta(1j,1) theta(lj,l)=theta(lj,jflag) theta(1j,j flag)=temp(1j,1) 760 continue go to 733 endif 740 continue 733 continue 730 continue endif 710 continue c*************************************************************** c normalizing theta and r vectors c************************************************************** do 810 i=l,numec do 820 j=l,2*nel if(abs(theta(j,i)).ne.1.O.and. * theta(j,i).ne.0.0)then fact=abs(1./theta(j,i)) do 830 jj=l,2*nel theta(jj,i)=theta(jj,i)*fact 830 continue do 840 jj=l,ndof rb(jj ,i)=rb(j j ,i) *fact 840 continue go to 734 endif 820 continue 734 continue 810 continue c************************************************************** c transpose theta and r matrices c************************************************************** 144 substantially the size of the problem. However, the benefits of this change could be significant. Changing the nonlinear analysis method from secant stiffness approach to a tangent stiffness approach could be another solution to the lack of convergence. In this case a two stage process would be adopted. The first would consist of a linear optimization up to the formation of a hinge followed by a phase with a sequence of incremental loading and optimization procedures until convergence was obtained. 81 ductile systems there is a large list of procedures, most of them with limited application (67). Some methods for redundant systems involve the determination of all collapse modes and their respective probability of failure. To obtain all collapse modes the fundamental mechanisms are determined and a Monte Carlo simulation is performed to generate all others. Afterwards the respective probabilities of failure are determined. This approach, although accurate, is very demanding in computational effort if the system is complex, and consequently is used mostly to validate the performance of other methods. In redundant ductile systems a variation of the Monte Carlo approach PNET or Point Estimate of System Collapse Probability is used. This consists in linearly combining the fundamental failure modes with the coefficients as variables. An objective function representing the reliability index of that combination is minimized and the most probable failure mechanisms are defined. Concerning redundant structures with brittle or ductile elements, other approximation and bounding techniques have been developed and studied based on graph theory. Two of those approaches for obtaining the probability of failure are the failure mode approach and the stable configuration approach (68). Both methods require the determination of all possible failure modes and the use of algorithms based on graph theory. CHAPTER 4 STRUCTURAL ELEMENT RELIABILITY Introduction Design and checking of structures in the field of Civil Engineering has been traditionally based on deterministic analysis. Adequate dimensions, material properties and loads are assumed and an analysis is carried out to obtain the required evaluation. Nevertheless, variations of all these parameters and questions related to the structural model may impose a different behavior than expected (41). It must be emphasized that if there were no uncertainties related to the prediction of loads, materials and structure modeling, then the respective safety would be more easily guaranteed. For these reasons the use of probabilistic principles and methodologies in structural design has been increasing. Design for safety and performance should consider the conflict between safety and risk. The objective of probability concepts and methods is to develop a framework where the effects of these uncertainties are considered. Structural reliability has received the attention of several 54 99 both unitary costs by the cost of concrete and the result is as follows: where f(x) = (x Xj + xk 10) Lp Lp length of element p, p = l,m. Equality constraints, one for each global degree of freedom, are defined as hq 2r (kqr x3n+r) Rq/ <3"rl,...,m where kqr global stiffness coefficient; x3n+r global displacement r; Rq external force q. Inequality constraints that control the maximum global displacements and impose a minimum element reliability are as follows where 9i x3n+r dr i1,...,m gj = relj betaj, < 0 j=l,...n dr maximum absolute value of global displacement r; relj reliability index of element j; betaj minimum reliability index prescribed for element j. 164 c************************************************************* x=47./60.*b*fc y=0.004*es*aste-aste*fy z=-0.004*es*co*aste if((y*y).It.(4.*x*z))then y=sqrt(4.*x*z) end if vkd=(-y+sqrt(y*y-4.*x*z))/(2.*x) epcs=0.004*(vkd-co)/vkd if(epcs.ge.epsy)then epcs=epsy endif c*************************************************************** C CONCRETE FORCE IN REGION AB c*************************************************************** alphal=2./3. ccab=alphal*b*0.5*vkd*fc c* ************************************************************ * c CONCRETE FORCE IN REGION BC c****************************************************** ******* alpha2=0.9 ccbc=alpha2*b*0.5*vkd*fc (^* ************************************************************* * C DISTANCE OF CENTROID TO TOP IN AB C************************************************************** gama1=0.875*vkd c******** *************************************** **************** C DISTANCE OF CENTROID TO TOP IN BC C***************************************************** ********* gama2=0.259255*vkd c********************************************************** C COEFFICIENTS FOR FAILURE FUNCTION c************************************************************ al=(ccab*(dd-gamal)+ccbc*(dd-gama2))/fc a2=-l. c*********************************************************** c COSINE DIRECTORS c******************************************** ************** tetal=al*sigmal*fc teta2=a2*sigma2*vm c************************************************************ C INDEPENDENT TERM c*********************************************************** fps=0.004*es*(vkd-co)/vkd bi=aste*fps*(dd-co) c*********************************************************** C RELIABILITY INDEX Q*********************************************************** beta(kel)=(al*fc+a2*vm+bi)/sqrt(tetal*tetal+teta2*teta2) c************************************************************ C ULTIMATE MOMENT AND ROTATION c************************************************************ vmu(kel)=al*fc+bi phiu=0.004/vkd if((4.*phiy).It.phiu)then vmu(kel)=(vmu(kel)-vmy)/(phiu-phiy)*3.*phiy+vmy 6 Methods In the last twenty years researchers have made considerable advances in developing techniques of optimum design. Research and exploration of these methods were mainly developed in the aeronautical and mechanical industries, where the need for more economical and efficient final products was extremely important. More recently, with the availability of increasing computer capabilities, civil engineering researchers and designers have increased their participation in structural optimization following the lines defined by the other engineering disciplines. Optimization methods are, nevertheless, common to these different engineering design areas and are mainly divided in two groups. These are commonly known by the names Optimality Criteria and Mathematical Programming (6). Another area in structural optimization researched by a few scientists is based on duality theory concepts, and is an attempt to unify the two basic methodologies (7). Optimality Criteria methods are based on an iterative approach where the conditions for an optimum solution are previously defined. The concept can be used as the basis for the selection of a structure with minimum volume. This methodology derives from the extreme principles of structural mechanics and has been limited to simple structural forms and loading conditions. The formulation can be mathematically expressed as follows: 73 where ^1 - mean value of fc? ^2 - mean value of Me; <*1 - standard deviation of fc; cr2 standard deviation of Me. Replacing the standard normal variables in the flexural failure function the expression assumes the following form: z = alalul + a2a2u2 + alM-i + a2p.2 + b Then the reliability index for each element is given by the distance from the standardized failure function to the origin of the standardized basic variables as follows: 6 = (axm + a2u2 + b) /(a 1o1 + a2o2)% 60 engineers where reliability theory and probabilistic methods are the basis of their code provisions. These Level 1 methods could be replaced by the Level 2 methods if an agreement was obtained in the following issues: selection of basic random variables for each specific problem, their distribution types and relative statistical parameters; form of the various limit state equations and choice of models; operational reliability levels to be adopted in different design situations. It must be emphasized that the advantage of Level 1 schemes over Level 2 are their great operational simplicity due to the use of fixed and constant partial safety factors for a given class of design situations. The main disadvantage of Level 1 is the selection of partial safety factors for a given structural class in such a way that the efficiency of the method proposed is satisfactory. It must assure that the deviation of the reliability of a design made on the basis of the adopted coefficients from the desired reliability level laid down in the code is acceptable. Two Dimensional Space Example Let R and S be two random variables, where R defines strength and S the load. Then the limit state function z shown in Figure 4.1 is defined as 91 mi. mi. mi. mi. mi. Mode 3 Mode 4 Mode 5 mi. Mode 6 mi. mi. mi. Mode 7 Mode 8 Figure 5.4. System failure modes. 212 do 305 k=l,m c(k)=a(k,iflag) a(k,iflag)=a(k,i) a(k,i)=c(k) 305 continue do 310 k=l,nt c(k)=b(k,iflag) b(k,iflag)=b(k,i) b(k,i)=c(k) 310 continue do 280 j=i+l,nt if (abs(a(i,j)).gt.0.00001) then fact=-a(i,j)/a(i,i) do 290 kk=l,m a(kk,j)=a(kk,j)+a(kk,i)*fact 290 continue do 291 kk=l,nt b(kk,j)=b(kk,j)+b(kk,i)*fact 291 continue endif 280 continue 200 continue numec=nt-3*n c* ************************************************************* * c Forming bl c***************************************************** ********** lcount=nt-numec+l ki=l do 800 i=lcount,nt do 810 j=l,nt bl(j/ki)=b(j,i) 810 continue ki=ki+l 800 continue c************************************************************** c Creating Theta matrix c******************* ****************************************** do 156 j=l,numec do 157 i=l,2*n k=iqh+i theta(i,j)=bl(k,j) 157 continue 156 continue c*************************************************************** c Creating virtual displacements c***************************************************** ********* do 169 j=l,numec do 158 i=l,iqh r(i,j)=bl(i,j) 158 continue c************************************************************ c Adding joint mechanisms c************************************************************** do 161 i=3,iqh,3 if (abs(r(i,j)).gt.0.000001)then r(i,j)=0.0 29 XI area of beam X2 inertia of beam X3 horizontal tip displacement X4 vertical tip displacement X5 tip rotation VARIABLE INITIAL FINAL XI (in2) i.O 6.65 X2 (in4) 1.0 78625 X3 (in) 0.4 0.500 X4 (in) 0.4 0.353 X5 (rad) 0.4 0.006 Figure 2.2. Cantilever beam. 52 Spring Moment-Rotation Diagram Hu Ultimate moment My Yielding moment Kl 10e30 K2 (Hu My)/(Ou fly) Ksec Spring stiffness for M > My Figure 3.6. Secant spring stiffness. 223 al=a2 a2=a resl=res2 res2=res go to 100 endif arml=dd-a+2./3.*xl arm2=dd-gama*(a-xl) vmy=ccl*arail+epcs*es*(dd-co)*aste+cc2*arm2 phiy=epsy/(dd-a) c************************************************ ******** c LINEAR APPROXIMATION Q******************************************************** ro=aste/(b*d) vn=es/ec vk=sqrt(4.*ro*ro*vn*vn+2.*(ro+ro*co/d)*vn)-2.*ro*vn vkd=vk*d exc=epsy/(d-vkd)*vkd cc=0.5*vkd*b*ec*exc excs=exc/vkd*(vkd-co) fps=excs*es cs=fps*aste vmy=cc*(d-vkd/3.)+cs*(d-co) phiy=exc/vkd write(*,*)'cc=1,cc,'cs=',cs,'exc=',exc,'excs=',exes write(*,*)'linear approx','vmy=',vmy,'phiy=',phiy stop end 158 tvag=tvag+clag(k)*psi+psi*psi*rp 200 continue c*************************************************************** C VALUE OF LAGRANGIAN FUNCTION C*************************************************************** call valobf (n,ntot,vof,x,cl) avof = cv*vof vlag = avof+tvah+tvag return end subroutine modsti(area,ec,vki,vkj,cl,tinert,cosl,cos2) implicit double precision (a-h,o-z) common /esq/ u(6),ck(6,6),vksi(100),vksj(100) C* ************************************************************* * C FLEXIBILITY MATRIX (2x2) C********************************************************* ****** do 10 i=l,6 do 20 j=l,6 ck(i,j)=0.0 20 continue 10 continue x=cl/(3*ec*tinert)+l./vki y=cl/(3*ec*tinert)+l./vkj z=-cl/(6*ec*tinert) C*************************************************************** C INVERSION OF MATRIX C******* ********************* ******************** ******* * * * det=x*y-z*z a=y/det b=-z/det c=b d=x/det C************************************************************** C EXPANDED MATRIX (6x6) Q************************************************************ * ckll=ec*area/cl ckl4=-ckll ck41=-ckll ck44=ckll ck22=(a+b+c+d)/(cl*cl) ck25=-ck22 ck52=ck25 Ck55=ck22 ck23=(a+c)/cl Ck53=-Ck23 ck26=(b+d)/cl Ck56=-Ck26 ck33=a APPENDIX A AUGMENTED LAGRANGIAN SUBROUTINES 214 do 121 ikl=l,iqh vahk(ikl)=0.0 121 continue n=nel n3=n+n+n do 100 k = l,n sigma2=0.0 do 111 ipo=l,6 u(ipo)=0.0 111 continue do 200 i = 1,6 m=lm(i,k) if (m.eq.O) go to 200 if(cvload(m).gt.sigma2)sigma2=cvload(m) 1 = n3 + m u(i) = x(l) 200 continue cl=cosl(k) c2=cos2(k) d2=-c2*u(1)+cl*u(2) d3=u(3) d5=-c2*U(4)+Cl*U(5) d6=U(6) c*************************************************************** c element forces c*************************************************************** c = cl(k) base = x(3*k-2) height = x(3*k-l) aste = x(3*k) area = base height tinert = area*height*height/12.0 al = ec*tinert/(c*c*c) call eley (ec,tinert,c,vksi(k),vksj(k),d2,d3,d5,d6, * fo3,fo6) sigmal=cvmu(k) c*************************************************************** c ultimate and yield moments c*************************************************************** call newmum(k,x,sigmal,sigma2,fo3,fo6,vksi(k),vksj(k), * d2,d5,d3,d6) c*************************************************************** c global modified stiffness c*************************************************************** call modsti(k,tinert,area,vksi(k),vksj(k)) do 300 1=1,6 j=lm(l,k) if (j.eq.0) go to 300 do 400 11=1,6 m = lm(ll,k) if (m.eq.O) go to 400 jj = n3+m vahk(j)=vahk(j)+ck(l,11)*x(jj) 400 continue 300 continue 100 continue 122 i j O- plastic hinge location i.j element node location A Element AB Ks = GJ/L where Ks Secant spring stiffness; G Shear modulus; J Torsional moment of inertia; L Element length. Figure 7.1. Displacement verification. 71 average standard variation for 68 a good quality control testing at the construction site is 550 psi. Using a 3500 psi specified compressive strength of concrete, f'c, the required average compressive strength of concrete is the larger of the following (51): f,cr = f/C + 1 34*CJ = 3500 + 1.34 550 = 4237 psi or f'cr = f'c + 2.33*cr 500 = 3500 + 2.33*550 500 = 4282 psi The coefficient of variation of f'c for this range of characteristic compressive strength is then given by V = a/f'c = 550/4282 = 0.128 and consequently the coefficient of variation of the concrete compressive flexural strength was adopted to be 0.15. External loads have different coefficients of variation for the different types of loads (52-53). For most design and construction in the United States a good estimate for the coefficient of variation of dead loads is 0.10. For the live loads the coefficient of variations are very high and range from 0.39 to 1.04. For that reason and since the building codes prescribe large values for live loads that exceed the mean value a single coefficient of variation of 0.15 was adopted for the combination of dead loads plus live loads. 210 implicit double precision (a-h,o-z) dimension a(100,100),b(100,100),c(100),cm(100,100), * cost(n),qa(100,100),sint(n),cl(n),q(100,100), * lm(6,n),am(100,100),bl(100,100),theta(200,100), * r(iqh,100) Ã‚Â£* ************************************************************ * c Constraint matrix for the structure q**************************************************** ********** do 300 i=l,3*n do 400 j=l,6*n cm(i,j)=0.0 q(i,j)=0.0 400 continue 300 continue do 60 k=l,n i=3*k j=6*k at=l.0/cl(k) im2=i-2 iml=i-l jml=j-1 jm2=j-2 jm3=j-3 jm4=j-4 jm5=j-5 cm(im2,j m5)=-l.0 cm(im2,jm2)=1.0 cm(iml,jm3)=l.0 cm(i,j)=1.0 cm(iml,jm4)=-at cm(iml,jml)=at cm(i,jm4)=-at cm(i,jml)=at 60 continue a*************************************************************it c Coordinate trnsformation matrix c************************************************************** do 70 k=l,n co=cost(k) si=sint(k) j=6*k do 80 i=l, 2 jum=j-3*i+l jdois=j-3*i+2 jtres=j-3*i+3 q( jum,jum)=co q(jum, jdois)=si q(jdois,jum)=-si q(jdois,jdois)=co q(jtres,jtres)=1.0 80 continue 70 continue Q* ************************************************************ * c Compatibilibity matrix from LM matrix c************************************************************** n6=6*n 31 106Kin lOSKin ELEMENT INITIAL FINAL 1 Area (in2) 1.0 25.4 Inertia (in4) 1.0 120224 2 Area (in2) 1.0 179 Inertia (in4) 1.0 5912 3 Area Cin2) 1.0 35.1 Inertia (in4) 1.0 17058 Figure 2.3. One bay frame. 226 Parameter indicating alteration of the limit of number of iterations. Line 24 Maximum number of consecutive iterations without objective function improvement, maximum number of consecutive Newton iterations, maximum number of completed one dimensional searches. Line 25 Parameter that controls the quantity of information in the output file. Line 26 Number indicating minimum printed information. Line 27 Indication that tangent vector extrapolation should be used for estimating initial values of basic variables. Line 28 Number of design variables iniatially included in the basis. Line 29 Numbers of design variables of the initial basis. Line 30 Parameter that indicates if new data should be read. 179 * 8003 8004 if(invmec(lmi).eq.m) go to 8004 continue end if nk=nk+l numele=numele+l invmec(n)=m go to 8000 continue end if endif 8002 continue 8001 continue 8000 continue c************************************************************* c control of system reliability c************************************************************* write(8,*) write(8,*)BETA MINIMAL FOR THE SYSTEM = ';betmin write(8,*) jflag=0 if(betmin.It.relind)then delta=(relind-betmin)/relind do 3891 i=l,nel do 3892 j=l,numele if(invmec(j).eq.i)then elerel(i)=(1.+delta)*elerel(i) endif 3892 continue 3891 continue jflag=l endif return end subroutine mecsys(n,iqh,cl,cost,sint,lm,numec,r,theta) implicit double precision (a-h,o-z) dimension a(100,100),b(100,100),c(100),cm(100,100), * cost(n),qa(100,100),sint(n),cl(n),q(100,100), * lm(6,n),am(100,100),bl(100,100),theta(200,100), * r(iqh,100) c************************************************************** c Constraint matrix for the structure c************************************************************** do 300 i=l,3*n do 400 j=l,6*n cm(i,j)=0.0 q(i,j)=0.0 400 continue 38 reaches the yielding value there is a large strain and section curvature increase. Geometric nonlinearities are then created by extra rotations of flexural elements from the cracking and steel yielding. A basic assumption in nonlinear analysis of reinforced concrete frames is that the element rotations with relation to the line defined by the nodes, chord rotations, are small and the theory for straight elements may be applied with some adaptations. The most popular analysis techniques are based on incremental loadings of the structure and are known by the initial stiffness and tangent stiffness methods. A technique based on the application of the entire load at a single step is known by the secant stiffness method. This last technique was chosen for the analysis of the structure since it is more adequate to the optimization formulation. Element Modeling Survey In the last three decades there have been many attempts to create a simplified beam model of the inelastic reinforced concrete element (31-33). The main objective for this research has been to advance a solution providing precise results within reasonable computational and memory storage limits. The study has a significant importance for the analysis of reinforced concrete structures submitted to dynamic loads (34-35). In these examples the moments at the ends are close to the ultimate allowable values. This 79 Series systems with n elements have n failure modes. Parallel systems with n elements have more than n failure modes. These failure modes in parallel systems are dependent on whether the failure type of the elements is brittle or ductile (62-63). For redundant brittle systems the failure of an element and consequent redistribution of the loads will provoke the system failure. In these cases the system behavior may be considered to be generally identical to that of as a series system. Probability of failure of a series system can be considered as the union of the elements probability of failure Pfs = P(Ui(Zi< 0)|i=l,n) where U union of events; Pfs probability of system failure; Zi failure function of element i. If the element failure functions are not correlated then the evaluation of PfS is relatively easy and may be assumed as Pfs = 1 *i=in('l P(ei=0)) where 7T product; ej_ = 0 if element is in a failure state, 123 5.000 lb Coefficients of variation fc 0-15 loads -0.15 Materials fc = 3,000 psi fy = 40.000 psi Hinge location Figure 7.2. Debug frame. 83 Truss Bars ^ Load FAILURE GRAPH 2F F Bar Failure Figure 5.2. Failure graph. 18 Displacements allowed were based on serviceability requirements like cracking and relative story drift. The reliability indices were based on usual values of probability of failure used in design codes. Only the flexural behavior of the frames was analyzed since it is the most important for usual structures and the members were modeled as beam elements. Inelastic behavior of the structure due to the material nonlinearities imposes a change of the global stiffness terms independently of those dictated by the alterations of the dimensions during the optimization search. For that reason, the reinforced concrete element was modeled as a linear elastic beam with nonlinear rotational springs at each end. Rotational spring stiffness was considered infinite when the moment was below the yielding moment. Above that value the element stiffness was inverted to its flexibility and the inverse of the secant spring stiffness value was added to the corresponding diagonal terms. Spring stiffness was calculated using the secant value of the bilinear moment-rotation diagram corresponding to the current global rotation. Values of the yielding and ultimate moments were obtained by integrating the actual stress-strain diagram for the compressive force in the concrete. The corresponding rotation at a hinge was calculated by integrating the curvature diagram along the element. 84 adopted was conceived by Watwood (15). It is an automatic tool to generate all failure mechanisms with one degree of freedom, or elementary failure mechanisms, of a given frame. The set of these mechanisms and all their linear combinations constitute all possible collapse configurations (70). The technique is relatively simple to use since the input data for this method is the same for traditional elastic analysis like joint and element information. Elementary failure mechanisms are dictated by the geometry of the structure and potential hinge locations created by the external load configuration. Hinge locations are considered at the end of each member. In the case where there are loads in the middle of the element, they are also considered at the points of concentrated or discretized loads. The element axial collapse is not considered in this formulation although it was included in the original version. Element global displacements of a planar frame form a vector with six variables, {S}. Using a cartesian referential set of axes x and y the displacements, to Sg, may be represented as in Figure 5.3. Element deformation parameters may be defined by three independent quantities S'i displacement about node i; S'2 ~ rotation of node i; S'3 rotation of node j. 118 avoid possible precocious eliminations the fundamental mechanisms are ordered such that those that involve external work are placed before joint mechanisms. To facilitate the combination of the mechanisms all virtual displacements are scaled so that all hinge rotations are unitary. At the end, if the mechanism with lower reliability index is not acceptable, the elements with hinges that belong to this failure mechanism have their required element reliability indices increased. The optimization process is restarted with these indices modified by the same percentage of the system reliability violation. Table 7.2. Debug frame: Augmented Lagrangian version Element Section 1 2 3 4 Total Total Initial Final Initial Final Reliability Base Height Area Base Height Area Index (in) (in) (in2) (in) (in) (in2) 5.0 10.0 1.0 2.0* 6.0* 0.06* 4.4 5.0 10.0 1.0 2.0* 6.8 0.07* 4.3 5.0 10.0 1.0 2.1* 8.4 0.09* 5.2 5.0 10.0 1.0 2.0* 9.4 0.09* 5.7 Initial Cost.. * - lower bounds. Final Cost.... Global Displacements 1 2 3 4 5 6 7 8 9 (in) (in) (rad) (in) (in) (rad) (in) (in) (rad) .23 0.0 -.002 .23 -.041' 0.0 .23 -.003 0.0 1.1* 0.0 -.009 1.1* -.126 .004 1.1* -.008 -.006 Secant Spring Stiffness (lb.in/rad) 1, 3, 4, 5, 6, 7, 8 30 Hinge Number Spring Stiffness 10x10 2 6.7X108 225 User's Manual Generalized Reduced Gradient Method Example: Debug Frame Input File: DATA Line 1 Problem title. Line 2 Number of variables, number of constraints, number of equality constraints. Line 3 Number of variables with lower bounds. Line 4 to line 11 Variable number and respective lower bound. Line 12 Number of constraints with upper bounds. Line 13 to line 16 Constraint number and respective upper bound. Line 17 to line 19 Initial values of design variables. Line 20 Number of prescribed optimization parameters. Line 21 Constraint tolerance. Line 22 Convergence tolerance. Line 23 157 call jacequ(x,n,cl,lm,cosl,cos2,fc,ec,vn,co,epsy, * fy,ntot, iqh,vah,r,vahk,es,ecm,beta,cvmu,cvload,k, * vmu,vjac) call sol(iqh,vjac,r,vinv) do 5890 jgo=l,iqh x(jgo+n3)=vinv(jgo) 5890 continue 200 continue C********************************************************* C REINITIALIZE VALUES C******************************************************** k=0 call lagfun (vlag,tvah,r,x,cl,cosl,cos2,lm, * d,clah,vag,clag,vah,ch,eg,vo f,vahk,vn,co,epsy,fy,beta, * cvmu,cvload,k,vmu,rph) return end subroutine lagfun (vlag,tvah,r,x,cl,cosl,cos2,lm, * d,clah,vag,clag,vah,ch,cg,vof,vahk,vn,co,epsy,fy,beta, * cvmu,cvload,kl,vmu,rph) implicit double precision (a-h,o-z) dimension x(ntot),cl(n),cosl(n),cos2(n),lm(6,n),d(iqg), * clah(iqh), clag(iqgn), vag(iqgn), r(iqh), vah(iqh), * ch(iqh), cg(iqgn), vahk(iqh), beta(n), * cvmu(n), cvload(iqh),vmu(n) common /parr/ decfc,fcinc,cv,alpl,ec,rp,fc,es,ecm,relind common /pari/ iter,numey,niter,ga,iqh,iqg,n,ntot,iqgn common /esq/ u(6),ck(6,6),vksi(100),vksj(100) tvah =0.0 tvag = 0.0 rp2 = 2.0*rp C************************************************************** C EQUALITY CONSTRAINTS C************************************************************** call equeon (x,n,cl,lm,cosl,cos2,fc,ec,vn,co,epsy, * fy,ntot,iqh,vah,r,vahk,es,ecm,beta,cvmu,cvload,kl,vmu) do 100 k=l,iqh vc=vah(k)*ch(k) tvah=tvah+clah(k)*vc+rph*vc*vc 100 continue C*************************************************************** C DISPLACEMENT CONSTRAINTS C************************************************************** call inecon (iqg,n,ntot,vag,x,d,relind,beta) do 200 k=l,iqgn v=vag(k)*cg(k) z=-clag(k)/rp2 psi=max(v,z) 46 The force in the compressed steel is given by where cs ~ As fcs As steel area; fcs stress in compressed steel. The force in the steel under tension is determined by Ts As where fy yielding steel stress. For instance, the internal ultimate moment is given by the moments of these three internal forces about the top compressed fiber. For that reason a parameter Cl, that defines the centroid of the concrete compressive stress diagram, is introduced as Cl = 1 - ca ec fc dec /(Ã‚Â£ca eca fc dec) 0 These parameters, a and il, when the ultimate concrete strain is defined as ec = 0.004, become a = 2/3 (region AB) n = 3/8 (region AB) a = 0.9 (region BC) n = 0.51851 (region BC) where the regions AB and BC are defined in Figure 3.2. The section flexural strength, Mji, may be defined as 125 The displacements satisfied the global equilibrium equations and all the constraints were satisfied. The optimization problem with material nonlinear behavior was solved using two minimization techniques described in Chapter 6. Results of the version using Hooke and Jeeves method are listed in Table 7.2. However, evaluation of element moments did not correspond to the location of the hinges, i.e., the element moment values in some hinges were below the yielding moment value. There was no force equilibrium in some of the nodes. An extensive set of initial design points and optimization parameters were tested with negative results. For that reason the optimization technique was tentatively replaced by the Generalized Reduced Gradient. First attempt to optimize with the Generalized Reduced Gradient assuming nonlinear material behavior, was with a secant spring stiffness equal to the ratio between the yielding moment and the yielding rotation. The option of using this spring stiffness had the advantage of avoiding the oscillation of the spring stiffness between the rigid and lower values. The implementation resulting from this choice was named yielding stiffness. Although providing incorrect displacements as the spring stiffness values did not represent the true material behavior, the yielding stiffness would model a situation somewhere between the linear and the nonlinear material behavior. The results are presented in Table 7.3. After adequate analysis it was ACKNOWLEDGEMENTS I want to express my sincerest gratitude to Dr. Marc Hoit for monitoring my research, for the inventive ideas and for his constant support. I am specially thankful to Dr. Fernando Fagundo for the productive discussions, his friendship and his vigorous encouragement. I owe to these two my best recollections from the University of Florida. My sincerest appreciation is extended to Dr. Prabhat Hajela for the teachings and the careful reading of my dissertation. My indebtment goes also to Dr. Clifford Hays and Dr. John Lybas for the useful conversations and their activity in my committee. I am also grateful to Dr. David Bloomquist for his support to initiate my geotechnical reliability research and his continuous disposition to help. I would also like to acknowledge my appreciation to the Fulbright Comission and to the Department of Civil Engineering of the University of Florida for their financial aid and the opportunity to study the fascinating area of structural optimization. My thanks go also to the Department of Civil Engineering of the University of Porto, specially to Dr. Adao da Fonseca, for giving me the possibility to research and study in the USA. ii 175 if(lj.It.-0.1)then kmu=-l lj=abs(lj) endif do 400 kk=l,nel kkj=2*kk-l thesum(kk)=thesum(kk)+theta(lj ,kkj) *kmu thesul(kk)=thesul(kk)+theta(lj,kkj+l)*kmu 400 continue kmu=l 300 continue c************************************************************** c reliability of combined mechanisms c (external work) c************************************************************** do 372 l=l,nucome lcl=lc(j+l-l) if(lcl.It.-0.1)then kmu=-l lcl=abs(lcl) endif do 472 kk=l,ndof if(abs(p(kk)).lt.0.001)go to 472 dispsu(kk)=dispsu(kk)+rb(lcl,kk)*kmu 472 continue kmu=l 372 continue 5564 continue c************************************************************* c combination with fundamental mechanisms c (internal work) c************************************************************* do 100 k=l,numec vmeanr=0.0 vmeanl=0.0 stdevr=0.0 stdevl=0.0 vmanrm=0.0 vmanlm=0.0 stdvrm=0.0 stdvlm=0.0 do 499 lll=j,j+nucome-1 if(abs(lc(lll)).ge.k)go to 100 499 continue if(nucome.gt.1)then if (icontr.eq.1)then becomi=400. becopl=500. go to 5574 endif endif thesu=0.0 thesu2=0.0 thesui=0.0 theslm=0.0 do 600 kk=l,nel 98 during the implementation and testing are presented and discussed. Augmented Lagrangian Formulation The set of design variables is divided in two main groups. These are the dimensions and steel area, defining each element cross section, and the global displacements. The objective function is the cost of the structure as a function of the the volume of concrete and steel. Equality constraints are defined by the global equilibrium equations. Inequality constraints include the bounds on global displacements and the minimum element reliability. In a reinforced concrete portal frame the set of design variables x having n elements and m global degrees of freedom is partitioned as follows: x, i=l,4,...,3n-2 base of rectangular element section; Xj, j=2,5,...,3n-l height of rectangular element section; xjr, k=3,6,..., 3n area of steel on one side of section; xi, l=3n+l, ..., 3n+m global displacements. The objective function used in this formulation was defined using the average costs of cast in place concrete for reinforced concrete frames and main reinforcing steel (71). Combined relative cost function was obtained dividing 33 found to expedite calculations to obtain acceptable initial values. The method of steepest descent makes use of the gradient of the pseudo-objective function. The gradient vector represents the line along which there is the highest variation of the pseudo-objective function at the actual design point. Moving in the direction defined by the negative of the gradient vector is expected to decrease the value of the pseudo-objective function. This direction is called the steepest descent. A graphical representation of the method is displayed in Figure 2.4. Since the explicit formulation of the gradient of the pseudo-objective function was not practical to obtain, the gradient vector was obtained using a finite difference technique. To obtain the minimum point along the gradient direction another design point along that line is found such that it has a higher pseudo-objective function value. Then, the optimum value should lie in this interval and a line search is performed using the golden section method. The gradient vector was normalized to avoid numerical ill-conditioning. For the same reason, constraints and the design variables were also scaled. Numerical difficulties are predictable if just one of the constraint function, or the objective function, is of different magnitude than the rest of the terms or its rate of change is considerably different from the others. Scaling factors for each constraint were evaluated as the ratio between the gradient 168 end subroutine inecon (iqg,n,ntot,vag,x,d,relind,beta) implicit double precision (a-h,o-z) dimension vag(iqg+n),x(ntot),d(iqg),beta(n) nel3=3*n c* ******************************************************** * c displacements c********************************************************** do 100 k l,iqg j = nel3 + k vag(k) = abs(x(j)) / d(k) 1. 100 continue c************************************************* ********* c reliability c********************************************************** iqgpl=iqg+l iqgn=iqg+n do 200 k=iqgpl,iqgn vag(k) = relind/beta(k-iqg)-1. 200 continue return end subroutine inputd (cl,cosl,cos2,iqh,jdir,jm,lm,n,nj, * nol,no2,r,xc,yc) implicit double precision (a-h,o-z) dimension nol(n),no2(n),jdir(3),xc(nj),yc(nj), * jm(6,nj),lm(6,n),cl(n),cosl(n),cos2(n),r(iqh) c********************************************************** c filling lm matrix c********************************************************** do 600 i = l,n j = nol(i) k = no2(i) do 600 1 = 1,3 lm(l,i) = jm (1, j) lm(l+3,i) = jm(l,k) 600 continue c********************************************************** c geometric characteristics Q* ***************** ************************************** * do 800 ii = l,n 76 This option is not commonly used for inhabited structures since it is difficult to evaluate the economic value of a structural failure where human life losses are expected. A more popular alternative is to include an additional constraint representing the maximum probability of failure allowed for the structure (60). The constraint for the system reliability will be of the type Pf(X) < Pm where Pf probability of system failure; X vector of design variables; Pm allowable probability of system failure. When performing structural optimization one may consider serviceability and ultimate limit states. This possibility leads to another type of formulation where the objective function and constraints for these limit states are considered simultaneously (61). This type of problems are called reliability-based optimization and can be summarized as follows: subject to Minimize CD Gi(X) < 0, i=l,m Pu ^ puo ps ^ pso 182 do 290 kk=l,m a(kk,j)=a(kk,j)+a(kk,i)*fact 290 continue do 291 kk=l,nt b(kk,j)=b(kk,j)+b(kk,i)*fact 291 continue endif 280 continue 200 continue numec=nt-3*n c***************************************************** ********** c Forming bl c*************************************************************** lcount=nt-numec+l ki=l do 800 i=lcount,nt do 810 j=l,nt bl(j ,ki)=b(j ,i) 810 continue ki=ki+l 800 continue q************************************************************** c Creating Theta matrix c************************************************************* do 156 j=l,numec do 157 i=l,2*n k=iqh+i theta(i,j)=bl(k,j) 157 continue 156 continue c******** ********************************************** ********* c Creating virtual displacements c****************************************************** ******** do 169 j=l,numec do 158 i=l,iqh r(i,j)=bl(i,j) 158 continue c************************************************************ c Adding joint mechanisms c******************************************************* ******* do 161 i=3,iqh/3 if (abs(r(i,j)).gt.0.000001)then r(i,j)=0.0 do 162 k=l,n if(lm(3,k).eq.i)then lpo=2*k-l theta(lpo,j)=l endif if(lm(6,k).eq.i)then lpo=2*k theta(lpo,j)=1 endif 162 continue endif 161 continue 169 continue SERIES SYSTEM PARALLEL SYSTEM Y Figure 5.1. System models. 42 and related altered element stiffness are simulated by the linear element with nonlinear rotational springs at the extremities. Inelastic rotations of reinforced concrete hinges at the element ends are determined as a function of the respective moment-curvature relationship for each element. These curves are redefined every time any element sectional properties changes during the optimization process since the ultimate and yielding moments also change. A typical moment-curvature diagram for reinforced concrete elements is bilinear. It is obtained assuming material stress-strain curves that are parabolic-linear for the concrete and bilinear for the reinforcing steel as shown in Figure 3.2 (28). The stress in the concrete is designated by fc and the stress in the steel reinforcement is represented by fs. The algorithm used to compute the moment corresponding to a certain strain diagram is an iterative Newton based iteration that determines the depth of the neutral axis guaranteeing equilibrium of the internal forces. Then, after determining the internal coupled forces the related moment is computed. All reinforced concrete elements are doubly reinforced with equal areas of steel on both sides. This assumption is valid for columns and acceptable for beams since in continuous frames there are moments of different sign along the beams. Evaluation of the moments for each reinforced concrete section was based on the exact internal equilibrium equations as follows: 206 499 continue if(nucome.gt.1)then if (icontr.eq.l)then becomi=400. becopl=500. go to 5574 endif endif thesu=0.0 thesu2=0.0 thesui=0.0 theslm=0.0 do 600 kk=l,nel kkj=2*kk-l thesu=thesum(kk)+theta(k,kkj) thesu2=thesul(kk)+theta(k,kkj+1) term=(abs(thesu)+abs(thesu2))*cvmu(kk) * *vmu(kk) vmeanr=term/cvmu (kk) +vmeanr stdevr=stdevr+term*term thesui=thesum(kk)-theta(k,kkj) theslm=thesul(kk)-theta(k,kkj+1) termm=(abs(thesui)tabs(theslm))*cvmu(kk) * *vmu(kk) vmannu=tennm/cvmu (kk) +vmannn stdvrm=stdvrm+tennin*termm 600 continue c************************************************************* c combination with fundamental mechanisms c (external work) Q************************************************************* do 672 kk=l,ndof if(abs(p(kk)).lt.0.001)go to 672 dispkk=(dispsu(kk)+rb(k,kk))*p(kk) vmeanl=vmeanl+dispkk stdevl=stdevl+dispkk*cvload(kk)*dispkk * *cvload(kk) dispkm=(dispsu(kk)-rb(k,kk))*p(kk) vmanlm=vmanlm+dispkm stdvlm=stdvlm+dispkm*cvload(kk)*dispkm * *cvload(kk) 672 5574 138 continue becopl=(vmeanr-vmeanl)/sqrt(stdevr+stdevl) becomi=(vmanrm-vmanlm)/sqrt(stdvrm+stdvlm) continue if(becomi.lt.becopl) then do 138 lk=l,nucome ltemp=ltemp+l lc(ltemp)=lc(j+lk-1) continue 1temp=ltemp+1 lc(ltemp)=-k lctlp=-k becote=becomi else do 139 lk=l,nucome 154 c = cl(kel) 45 continue base = x(3*kel-2) height = x(3*kel-l) aste = x(3*kel) area = base height tinert = area*height*height/12.0 al = ec*tinert/(c*c*c) fo3 = al*(6.0*c*(d2-d5) + 2.0*c*c*(2.0*d3+d6)) fo6 = al*(6.0*c*(d2-d5) + 2.0*c*c*(d3+2.0*d6)) C**************************************************************** C ULTIMATE AND YIELD MOMENTS C**************************************************************** call mumy(base,height,aste,vn,co,epsy,ec,fc,fy, * vksi(kel),vksj(kel),c,fo3,fo6,es,ecm,betak, * sigmal,sigma2,tinert,kl,vmuk,vmy,ijflag) if(vksi(kel).It.10e20.or.vksj(kel).It.l0e20)then call eley(ec,tinert,c,vksi(kel),vksj(kel), * d2,d3,d5,d6,fo3,fo6) call mumy(base,height,aste,vn,co,epsy,ec,fc,fy, * vksi(kel),vksj(kel),c,fo3,fo6,es,ecm,betak, * sigmal,sigma2,tinert,kl,vmuk,vmy,ijflag) endif beta(kel)=betak vmu(kel)=vmuk C**************************************************************** C GLOBAL MODIFIED STIFFNESS Q******************************************************* ********* call modsti(area,ec,vksi(kel),vksj(kel),c,tinert,cl,c2) do 300 1 = 1,6 j = lm (l,kel) if (j.eq.O) go to 300 do 400 11 = 1,6 m = lm (11,kel) if (m.eq.0) go to 400 jj = n3+m vahk(j)=vahk(j)+ck(l,ll)*x(jj) 400 continue 300 continue 100 continue Q********************************************************** C SUBTRACTION OF EXTERNAL GLOBAL FORCES C*********************************************************** rmax=0.01 do 510 ilj=l,iqh if(abs(r(ilj)).gt.rmax)rmax=abs(r(ilj)) 510 continue do 500 kpj = l,iqh if(abs(r(kpj)).It.0.0001)then vah(kpj)=vahk(kpj)/rmax go to 500 endif vah(kpj) = (vahk(kpj) r(kpj))/rmax 500 continue return end 27 HOOKE and JEEVES l Initial Point 4/5 Pattern Move 6 Final Point Figure 2.1. Pattern Search. 82 To exemplify the determination of all possible failure modes the initial step is to build a directed network, or directed graph, with all possible events involving element failures that will lead to a collapse. Each node represents a stable configuration and each branch corresponds to a element failure. Each path is a set of branches connecting the initial and final nodes. A cut of the graph is a set of branches containing only one branch from every path. A simple example is presented in Figure 5.2. Methods based on the determination of fundamental failure mechanisms using practical simplifications from graph theory have been implemented (69). The Beta unzipping method and the branch and bound method are two examples. The principal advantages are that they are precise and easy to use. The Beta unzipping method finds the significant collapse mechanisms using combinations of fundamental mechanisms and rejecting those combinations that are outside a prescribed interval. The branch and bound method selects all failure paths that have high probabilities of occurrence. Although less exact, the Beta unzipping method was chosen due to its simplicity and performance. Generation of Failure Modes To define all failure mechanisms, the first step consists of determining the set through manipulation of elementary failure mechanisms. To obtain these, the method 227 Example: Debug Frame Input File: DATA1 4,5 1,2 2.3 3.4 4.5 1,1,1,0,0 0,0,0,0,100 0,0,0,50,100 0,0,0,100,100 1,1,1,100,0 2,1,5000 3,2,-5000 '0,0,0 3000,40000,1 29e6,0.004 2 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 149 c************************************************************** call datini (clah,clag,ch,eg) c*************************************************************** c subroutine optimization c*************************************************************** call optimi (vlag,r,x,cl,cosl,cos2,lm,d,clah,xol, * vag,toll,clag,vah,ch,eg,xo,vahk,grad,xu,xl,xl,x2, * delta,alpha,alp,vaho,vago,vn,co,epsy,fy,ast,beta,theta, * numec,vmu,cvmu,rv,cvload,become,1c,thesum,thesuml, * dispsum,ni,nucomb,lct,xll,xl2,romin,rph,grad,vahold, * vjac,vinv,bfal,epsilo) c*************************************************************** c write final data q***************************************************************** write (8,880) write (8,890) iter write (8,910) vlag write (8,840) do 1000 k = l,n na = 3 k ne = na 1 no = ne 1 write (8,851) k, x(no), x(ne), x(na) 1000 continue write (8,960) do 1151 k = n21, ntot i = k 3*n write (8,970) i, x(k) 1151 continue write (8,920) do 1101 k = l,iqh write (8,930) k, vah(k) 1101 continue write (8,940) do 1200 k = l,iqg write (8,950) k, vag(k) 1200 continue write (8,944) do 1211 k=iqg+l,iqgn write(8,946)k-iqh,beta(k-iqh) 1211 continue write(8,945) do 1241 k=l,n write(8,947)k,vksi(k),vksj(k) 124 continue write(8,*) write(8,*)' LAGRANGIAN MULTIPLIERS EQUALITIES' write(8,*) do 1277 i=l,iqh write(8,966)clah(i) 1277 continue write(8,*) write(8,*)' LAGRANGIAN MULTIPLIERS INEQUALITIES' write(8,*) do 1278 i=l,iqh+n write(8,966)clag(i) 142 and using a secant approach relies upon the fact the exact secant spring stiffness value is obtained. There are certain approximations in the determination of the yielding and ultimate rotations, that define the moment rotation diagram from which the secant stiffness is evaluated. All these instabilities and approximations may create the lack of convergence that the results have shown. Future Work A good approach to improve the adequacy of the formulation assuming linear material behavior would be the determination of the proper values for the mean and standard deviation values of the external loads and concrete strength. Presently, there is a lack of information to allow a practical choice of these parameters for each particular design situation. More research should be done to examine the influence of including other statistical parameters such as the cross section dimensions, position of reinforcing steel, steel strength and load characteristics. Addition of other element effects will transform this formulation into a more complete optimization package. The main element force to be considered is the axial force that is decisive for column design. This will transform the system reliability evaluation and the element reliability constraints. Fundamental failure mechanisms will include axial failures coupled with flexural failures and there will 120 Hooke and Jeeves and the Generalized Reduced Gradient methods. The relevant results of these examples are presented in Tables 7.1 through 7.8. Result Verification Validation of results from the three types of frames described above was accomplished with a common strategy implemented at three levels. These were element reliability, compatibility with element moment rotation diagram, and global structure equilibrium and compatibility. Control of results was extensively performed for the debug frame and carefully administered in the other two cases. To evaluate the element reliability and the compatibility of the moment rotation diagram at the end or during the optimization process, a group of two programs was used. These computer programs called YIEL and ELTES are listed in Appendix C. The input data is composed of the dimensions of the cross section, the steel area reinforcement, the values of the secant stiffness of the springs, the length of the element and the global displacements of the element nodes. The output includes the element moments at the ends, the yielding and ultimate moments, the yielding and ultimate rotations, and the element reliability. These values are compared with those reported by the program results. TABLE OF CONTENTS Page ACKNOWLEDGEMENTS LIST OF TABLES vi LIST OF FIGURES V ABSTRACT viii CHAPTERS 1 STRUCTURAL OPTIMIZATION 1 Introduction 1 Historical Background 3 Methods 6 Typical Applications 8 Study Objectives 15 Summary 17 2 INTEGRATED OPTIMIZATION OF LINEAR FRAMES 21 Original Research 21 Augmented Lagrangian Function 22 Unconstrained Minimization Techniques 25 Final Results 28 Further Improvements 32 3 NONLINEAR REINFORCED CONCRETE ELEMENT 37 Introduction 37 Element Modeling Survey 38 Beam Element with Inelastic Hinges 40 Beam Element Stiffness 49 4 STRUCTURAL ELEMENT RELIABILITY 54 Introduction 54 Two Dimensional Space Example 60 Reinforced Concrete Element Reliability 69 5 SYSTEM RELIABILITY. 74 Introduction 74 System Reliability and Optimization 75 Methods 77 Generation of Failure Modes 82 Beta Unzipping Method 90 iv 192 dh2/dX! dh2/dx2 dh2/dx3 dh2/dx4 dh2/dx5 dh3/dxi dh3/dx3 dh3/dx5 x23 (-0.15x3 + x4); XiX22(-0.45x3 + 3x4); -0.15x^x23; xix23; 0; dh3/dx2 = dh3/dx4 = 0; -0.45x3; 1; Step c Initial Design Point and Initial Values Dependent variables d^t = {xi,x2,x3}; Independent variables dj^ = {X4,xs}; ddfc = {1,10,-0.1333} difc = {-0.02,0.4} grad = {100,10,0} grad fit = {0,0} where grad f is gradient of f; H = [ J I C ] where H is Hessian matrix of the equalities J -l -0.3 30 0 0 -150 -0.6-0.12 0 C -150 0 1000 0 0 1 Step D Recurrence Formulas dkt = {ddk I dik}; 219 endif return end subroutine coxticon (aste, dd, b, vmy, phiy) implicit double precision (a-h,o-z) dimension x(l) common /vgeom/ cl(100),cosl(100),cos2(100),lm(6,100), * cvmu(lOO) common /vload/ r(100)/cvload(100),jm(6,100),nol(100), * no2(100) common /xcord/ xc(100),yc(100),d(100),jdir(3) common /parr/ cV/ec^pjfC/es^cm^elind^o^y^psy common /pari/ iqh,iqg,nel,ntot,iqgn common /inequa/ vag(100),beta(100),u(6),vahk(100), * ck(6,6) common /equal/ vah(100),vmu(100) node=0 epso=0.002 c************************************************ ************** c c exc concrete strain c epcs compressive steel strain c epsy yield strain c c************************************************************** c first value for a c************************************************************** al=dd/2. exc=al*epsy/(dd-al) epcs=exc*(al-co)/al t=fy*aste cs=epcs*es*aste eces=exc/epso alpha=eces-eces*eces/3. cc=alpha*fc*b*al resl=cc+cs-t c******************************************************** ****** c second value for a c* ************************************************************ * a2=0.25*dd exc=a2*epsy/(dd-a2) epcs=exc*(a2-co)/a2 cc=fc*a2*alpha*b eces=exc/epso cs=epcs*es*aste res2=cc+cs-t c*************************************************************** c newton iteration c*************************************************************** 100 a=a2-res2*(a2-al)/(res2-resl) Figure 6.2. Augmented Lagrangian version flowchart. 25 L(X,,V) = f(X) +UH+PHH+VG' +P G'G' where u, v lagrangian multipliers; P penalty factor; G' maximum of (G, -v/2P}. The optimization procedure consists of several cycles of unconstrained minimization of the pseudo-objective function. The values of the lagrangian multipliers are kept constant during each cycle of the unconstrained minimization. At the end of an unconstrained minimization cycle, the multipliers are updated using an appropriate rule (12). The procedure is repeated for successive cycles until there is no significant change of the objective function. At this point the primal and dual optima have been found and the algorithm stops. Unconstrained Minimization Techniques Initially the technique used for the unconstrained minimization of the augmented lagrangian function was a zero order method referred to as the Hooke and Jeeves method or Pattern Search. The classification as a zero order method means that it does not utilize any information about the form or shape of the function. After the phase when stress constraints were added, a first order method, Steepest Descent, was tested as an improvement in the algorithm's 209 C************************************************************* nk=l numele=0 do 8000 m=l,nel do 8001 k=l,numec do 8002 l=l,nummec if(abs(locmec(l)).eq.k)then * * 8003 8004 endif 8002 continue 8001 continue 8000 continue ml=2*m-l m2=2*m if(abs(theta(ml,k)).gt.0.0. or.abs(theta(m2,k)).gt.0.0)then if(nk.gt.1)then do 8003 lmi=l,nk if(invmec(lmi).eq.m) go to 8004 continue endif nk=nk+l nume1e=nume1e+1 invmec(n)=m go to 8000 continue endif c************************************************************* c control of system reliability c************************************************************* write(3,*) write(3,*)'BETA MINIMAL FOR THE SYSTEM = ',betmin write(3,*) jflag=0 if(betmin.It.relind)then delta=(relind-betmin)/relind do 3891 i=l,nel do 3892 j=l,numele if(invmec(j).eq.i)then elerel(i)=(1.+delta)*elerel(i) endif 3892 continue 3891 continue jflag=l endif return end subroutine mecsys(n,iqh,cl,cost,sint,lm,numec,r,theta) 89 release corresponds to an addition of an external global degree of freedom. Addition of external degrees of freedom is done by replacing a row in matrix [Q] [A] with zeros. The changed rows correspond to the element degrees of freedom S3 and Sg, the node rotations. For each row that is replaced, a unit column vector is added to the matrix [Q] [A] with a 1 in the row that has been replaced. The dimensionality of {r} is increased by the number of rows replaced in [Q] [A]. The total is a set of extra columns with a dimension that is twice the number of elements. The homogeneous system becomes [C] ([Q] [A])* {ra} = [B1] {ra} = 0 where ([Q] [A])* matrix with extra 2n columns; {ra} vector of increased global degrees of freedom. Matrix [B'] is not square and has a greater number of columns than the number of rows. The solution of the system of homogeneous equations above may be obtained using a technique similar to that when solving for redundant unknowns in the force method. Difference between number of rows and number of columns is the number of independent solutions, that coincides with the number of elementary mechanisms. Suppose the rank of [B'] is the number of columns, m, and that the number of columns is p. In this 102 order method because it does not rely on information about the shape of the function obtained from derivatives. Methods based on the second derivatives were discarded since the problem was highly nonlinear and inequality constraints had discontinuous second derivatives. Conjugate Gradient is based on obtaining consecutive directions that are linearly independent, thus accelerating the search. The algorithm for the method is summarized as follows Step 1: Calculate grad f(xk); Step 2: dk = -grad f(x); Step 3: Find ak so that f(xk + ak.dk) = rain; Step 4: xk+1 = xk + ak.xk Step 5: Check convergence. If converged, stop. Step 6: dk+1 = dk + [grad f(xk+1)2/grad f(xk)2].dk Go to 3. Conjugate Gradient method proved to be unsuitable for the type of function presented. Progress in the minimization was minimal due to the ridge-type shape of the function. Whenever the process started at any point where 112 As specified before, some subroutines were directly used from the Augmented Lagrangian formulation while others had to be adapted or created. Since the data transfer between subroutines in the Generalized Reduced Gradient program was made through common data blocks, the same methodology was used for most of the added subroutines. The flowchart of this package is presented in Figure 6.3. An example of the input data files and the listing of the new or modified subroutines is presented in Appendix C. The unmodified subroutines perform the same tasks as described before. Essential structure of this program is the same as presented by the authors of the optimization package. There is a program, OPTIMI, that calls the main subroutines PRINCI, DATAIN, GRG and OUTRES. PRINCI reads the initial data from file DATA1 that is not abridged by the typical input data of the optimization package, which is read in subroutine DATAIN. The subroutine GRG performs the problem optimization calling other subroutines. The only subroutine written for this implementation was GCOMP that computes the values of the equality constraints, the inequality constraints, and the objective function. The system reliability was evaluated and the process was restarted if the results were unsatisfactory. Subroutine OUTRES writes the final results of each optimization run to a file RESULT. This subroutine was also modified to include the relevant APPENDIX C GENERALIZED REDUCED GRADIENT SUBROUTINES LIST OF TABLES Table Page 7.1. Debug frame (GRG): linear version results 124 7.2. Debug frame: Augmented Lagrangian version 126 7.3. Debug frame (GRG): yielding stiffness results 127 7.4. Debug frame (GRG): secant stiffness results 129 7.5. Debug frame: element moments 130 7.6. Compared frame: initial steel area reinforcement 133 7.7. Compared frame results 135 7.8. Building frame results 138 vi 178 do 6981 klp=l,ndof dispsu(kip)=0.0 6981 continue 200 continue c************************************************************* c control of maximum number of tree rows c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * lpti=lpti+nucomb(nucome+1) nucome=nucome+1 ni(nucome)=ni(nucome-1)+(nucome-1)*nucomb(nucome-1) if(nucome.It.(numax))go to 111 c************************************************************* c find minimum beta value c************************************************************* betmin=100 do 7890 i=l,lnumbe if(become(i).It.betmin)then betmin=become(i) itab=i end if 7890 continue c************************************************************* c find mechanisms involved c************************************************************* mecoun=0 mcomb=l do 7891 j=l,nucome-1 do 7895 jk=l,nucomb(j) mecoun=l+mecoun if(itab.eq.mecoun)then nistar=mcomb do 7893 1=1,j locmec(l)=lct(nistar+l-l) 7893 continue nummec=j go to 7894 endif mcomb=j +mcomb 7895 continue 7891 continue 7894 continue c************************************************************* c find elements involved c************************************************************* nk=l numele=0 do 8000 m=l,nel do 8001 k=l,numec do 8002 l=l,nummec if(abs(locmec(l)).eq.k)then ml=2*m-l m2=2*m if(abs(theta(ml,k)).gt.0.0. * or.abs(theta(m2,k)).gt.0.0)then if(nk.gt.1)then do 8003 lmi=l,nk 19 Element reliability was evaluated using a Level 2 method, i.e., an approximation to the evaluation of the exact probability of failure. The statistical variables considered were those assumed to have greater influence on the final result. These were the compressive strength of concrete and the external loads, assumed as normal distributed variables. The corresponding reliability index was calculated for constraint evaluation using the ultimate moment obtained from the integration of the respective strain diagram. Optimization techniques tested were based on the Augmented Lagrangian and the Generalized Reduced Gradient methods. The optimization problem was run, and after termination, the structure probability of failure was compared with the assigned value. If the result was not satisfactory, the process was restarted with updated values of the element reliability indices for the members involved in the most probable collapse mechanism. Evaluation of the system reliability was divided in two phases. First phase consisted of the identification of the elementary collapse mechanisms. In the second phase these elementary mechanisms were linearly combined to generate all significant mechanisms. System reliability was calculated considering the frame as a series system where each element is one of these mechanisms with higher probability of failure. 231 (15) Watwood, V. B., Mechanism Generation for Limit Analysis of Frames, Journal of Structural Division. ASCE, Vol. 109, ST1, 1979, pg. 1-15. (16) Thoft-Christensen, P., and Murotsu, Y., Application of Structural Systems Reliability Theory. Springer-Verlag, Berlin, 1986. (17) Ditlevsen, 0., Narrow Reliability Bounds for Structural Systems, Journal of Structural Mechanics. Vol. 7, 1979, pg. 435-451. (18) Schmit, L. A., and Fox, R. L. Fox, An Integrated Approach to Structural Synthesis and Analysis, AIAA Journal. Vol. 3, No. 6, 1960, pg. 1104-1112. (19) Grierson, D. E., and Schmit, L. A., Synthesis under Service and Ultimate Performance Constraints, Computers and Structures. Vol. 15, No. 4, 1982, pg. 405-417. (20) Haftka, R. T., Simultaneous Analysis and Design, AIAA Journal. Vol. 23, No. 7, 1985, pg. 1099-1103. (21) Hughes, T. J. R., Winger, J., Levit, I., and Tezduar, T. E., New Alternating Direction Procedures in Finite Element Analysis Based on EBE Approximate Factorization, Computer Method for Nonlinear Solids and Structural Mechanics. AMD, Vol. 54, 1983, pg. 75-109. (22) Burns, S. A., Simultaneous Design and Analysis Using Geometric Programming and the Integrated Formulation. Proceedings of the Swanson Analysis Systems, Newport Beach, California, 1987. (23) Marc I. Hoit, Alfredo V. Soeiro and Fernando E. Fagundo, Integrated Structural Sizing Optimization, Engineering Optimization. Vol. 12, No.3, 1987, pg. 207-218. (24) Soeiro, A., and Hoit, M., Sizing Optimization. Proceedings of the ASCE Structural Conference, Orlando, Florida, 1987. (25) Hoit, M., and Soeiro, A., Integrated Structural Optimization. Proceedings of the International Conference on Computational Engineering Science, Atlanta, 1988. (26) Soeiro, A., Integrated Analysis and Optimal Design. Thesis for the degree of Master of Engineering, University of Florida, Gainesville, Florida, 1986. (27) Avriel, M., Nonlinear Programming: Analysis and Methods. Prentice-Hall, Englewood Cliffs, New Jersey, 1976. (28) Park, R., and Paulay, T., Reinforced Concrete Structures. John Wiley and Sons, New York, 1975. OPTIMIZATION OF REINFORCED CONCRETE FRAMES USING INTEGRATED ANALYSIS AND RELIABILITY By ALFREDO V. SOEIRO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1989 ACKNOWLEDGEMENTS I want to express my sincerest gratitude to Dr. Marc Hoit for monitoring my research, for the inventive ideas and for his constant support. I am specially thankful to Dr. Fernando Fagundo for the productive discussions, his friendship and his vigorous encouragement. I owe to these two my best recollections from the University of Florida. My sincerest appreciation is extended to Dr. Prabhat Hajela for the teachings and the careful reading of my dissertation. My indebtment goes also to Dr. Clifford Hays and Dr. John Lybas for the useful conversations and their activity in my committee. I am also grateful to Dr. David Bloomquist for his support to initiate my geotechnical reliability research and his continuous disposition to help. I would also like to acknowledge my appreciation to the Fulbright Comission and to the Department of Civil Engineering of the University of Florida for their financial aid and the opportunity to study the fascinating area of structural optimization. My thanks go also to the Department of Civil Engineering of the University of Porto, specially to Dr. Adao da Fonseca, for giving me the possibility to research and study in the USA. ii My sincere appreciation and best remembrances go to my friends in the Gainesville Portuguese community and to my colleagues Jose, Joon, Lin and Prasit that helped smoothe the life contours created by the research work. Finally, my gratitude goes to my wife, Paula, for her work, her patience and her support throughout the whole period during which this dissertation was completed. iii TABLE OF CONTENTS Page ACKNOWLEDGEMENTS LIST OF TABLES vi LIST OF FIGURES V ABSTRACT viii CHAPTERS 1 STRUCTURAL OPTIMIZATION 1 Introduction 1 Historical Background 3 Methods 6 Typical Applications 8 Study Objectives 15 Summary 17 2 INTEGRATED OPTIMIZATION OF LINEAR FRAMES 21 Original Research 21 Augmented Lagrangian Function 22 Unconstrained Minimization Techniques 25 Final Results 28 Further Improvements 32 3 NONLINEAR REINFORCED CONCRETE ELEMENT 37 Introduction 37 Element Modeling Survey 38 Beam Element with Inelastic Hinges 40 Beam Element Stiffness 49 4 STRUCTURAL ELEMENT RELIABILITY 54 Introduction 54 Two Dimensional Space Example 60 Reinforced Concrete Element Reliability 69 5 SYSTEM RELIABILITY. 74 Introduction 74 System Reliability and Optimization 75 Methods 77 Generation of Failure Modes 82 Beta Unzipping Method 90 iv Page 6 PROCEDURE IMPLEMENTATION 97 Introduction 97 Augmented Lagrangian Formulation 98 Generalized Reduced Gradient 108 Reliability 114 7 EXAMPLES 119 Introduction 119 Result Verification.. 120 Debug Frame 121 Compared Frame 131 Building Frame... 136 8 CONCLUSIONS AND RECOMMENDATIONS 139 Linear Material Behavior 139 Nonlinear Material Behavior 141 Future Work 142 APPENDICES A AUGMENTED LAGRANGIAN SUBROUTINES 145 B GENERALIZED REDUCED GRADIENT EXAMPLE 189 C GENERALIZED REDUCED GRADIENT SUBROUTINES 195 REFERENCES 230 BIOGRAPHICAL SKETCH 236 V LIST OF TABLES Table Page 7.1. Debug frame (GRG): linear version results 124 7.2. Debug frame: Augmented Lagrangian version 126 7.3. Debug frame (GRG): yielding stiffness results 127 7.4. Debug frame (GRG): secant stiffness results 129 7.5. Debug frame: element moments 130 7.6. Compared frame: initial steel area reinforcement 133 7.7. Compared frame results 135 7.8. Building frame results 138 vi LIST OF FIGURES Figure Page 1.1. Implicit optimization 5 1.2. Element optimization 10 1.3. Truss optimization 11 1.4. System optimization 13 1.5. Geometry optimization 14 2.1. Pattern Search 27 2.2. Cantilever beam 29 2.3. One bay frame 31 2.4. Gradient method 3 4 3.1. Element model 41 3.2. Material behavior 43 3.3. Reinforced concrete section 45 3.4. Element deformation diagrams 48 3.5. Curvature integration 50 3.6. Secant spring stiffness 52 4.1. Design safety region 61 4.2. Probabilistic functions..... 64 4.3. Safety checks 66 4.4. Reliability index 68 5.1. System models 78 5.2. Failure graph 83 5.3. Element displacements definition 85 5.4. System failure modes 91 5.5. Combinatorial tree 96 6.1. Augmented lagrangian function plot 104 6.2. Augmented Lagrangian version flowchart 106 6.3. Generalized Reduced Gradient version flowchart... 113 6.4. Bilinear elastic-plastic diagram 117 7.1. Displacement verification 122 7.2. Debug frame 123 7.3. Compared frame 132 7.4. Building frame 137 B.l. Integrated optimization example... 191 vii Abstract of the Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy OPTIMIZATION OF REINFORCED CONCRETE FRAMES USING INTEGRATED ANALYSIS AND RELIABILITY By ALFREDO V. SOEIRO August 1989 Chairman: Dr. Marc I. Hoit Cochairman: Dr. Fernando E. Fagundo Major Department: Civil Engineering Simultaneous analysis and design were considered in the optimization of reinforced concrete frames. Frame elements had rectangular cross sections with double steel reinforcement. Design variables were the section dimensions, the area of steel reinforcement and the structure global displacements. Equality constraints were the equilibrium equations and inequality constraints were generated by element reliability requirements, code reinforcement ratios and section dimension bounds. Optimization strategies were based on the Augmented Lagrangian formulation and on the Generalized Reduced Gradient method. Reliability of the frames was considered at the element and system level. An element failure function was defined using moment forces and flexural strength. The random viii variables considered were flexural strength of concrete and external loads. System reliability was evaluated at the mechanism level using combinations of the elementary failure mechanisms. Optimization of the frames considering material nonlinear behavior was also investigated. Inclusion of this property was performed using a one-component model for the reinforced concrete element. Inelastic rotational springs were added to the ends of the linear elastic element. The element matrix was obtained by condensation of element elastic stiffness and secant spring stiffness. Three frames were researched. Respective results using linear material behavior were discussed. In these three cases the optimal solutions were found. Element reliability constraints were active and system reliability was satisfied. The integrated formulation was validated in the linear behavior range. The nonlinear material behavior results were presented for the smaller frame. ix CHAPTER 1 STRUCTURAL OPTIMIZATION Introduction Optimization is a state of mind that is always implicitly present in the structural engineering process. From experience engineers learn to recognize good initial dimension ratios so that their preliminary designs demand small changes through the iterative process and that elements are not overdesigned. The motivation behind this attitude is to create a structure that for given purposes is simultaneously useful and economic. Structural optimization theory tries to rationalize this methodology for several reasons. The main one is to reduce the design time, specially for repetitive projects. It provides a systematized logical design procedure and yields some design improvement over conventional methods. It tries to avoid bias due to engineering intuition and experience. It also increases the possibility of obtaining improved designs and requires a minimal amount of human- machine interaction. 1 2 There are, however, some limitations and disadvantages when using design optimization techniques. The first one is the increase in computational time when the number of design variables becomes large. Another disadvantage is that the applicability of the specific analysis program that results from the optimization formulation is generally limited to the particular purpose to which it was developed. A common inconvenience is that conceptual errors and incomplete formulations are frequent. Another drawback is that most optimization algorithms have difficulty in dealing with nonlinear and discontinuous functions and, hence, caution must be exercised when formulating the design problem. Another factor of concern is that the optimization algorithm does not guarantee convergence to the global optimum design, yielding on most occasions local optimum points. These facts lead to the conclusion that optimization results may often be misleading and, therefore, should always be examined. Therefore, some authors suggest that the word "optimization" in structural design should be replaced by "design improvement" as a better expression to materialize the root and outcome of this structural design activity (1). Nevertheless, there is an increasing recognition that it is a convenient and valuable tool to improve structural designs has been increasing among the designers community. 3 Historical Background Throughout time there have been various attempts to address structural optimization. The earliest ideas of optimum design can be found in Galileo's works concerning the bending strength of beams. Other eminent scientists like Bernouilli, Lagrange, Young, worked on structural optimum design based on applied mechanics concepts (2). These pioneering attempts were based on a close relation to the thoughts and accomplishments of structural mechanics. They started with hypotheses of stress distribution in flexural elements and ended with material fatigue laws. The accepted first work in structural optimization discusses layout theory, or structural topology. The paper focused on the grouping of truss bars that creates the minimum weight structure for a given set of loads and materials. The author of this primary achievement was Maxwell, in 1854, and Michell developed and publicized these concepts in 1904 (3). The practical application of these theorems was never accomplished since significant constraints were not included in the original works. Some procedures widely used by structural designers are nothing more than techniques of structural optimization. A well known example is the so-called Magnel's diagram (4). It is used to find the optimal eccentricity of the cable that leads to the smallest prestressing force without exceeding the limits imposed on the stresses in prestressed 4 concrete beams with excess capacity. This is a typical maximization problem in a linear design space, where the design variables are the eccentricity and the inverse of the cable prestressing force. The objective function is the value of the inverse of the cable prestressing force, and is to be maximized. The constraints represent the allowable stresses in tension and compression at the top and bottom of the cross-section. The problem is solved using a graphic representation of the problem, as shown in Figure 1.1, but could be solved numerically using the Simplex method. Numerical optimization methods and techniques have been widely researched and used in the operations research area, commonly known as Mathematical Programming. The practical application of these theories has been carried out in several areas for some decades like management, economic analysis, warfare, and industrial production. Lucien Schmit was the first to use nonlinear programming techniques in structural engineering design (5). The main purpose of structural optimization methods was to supply an automated tool to help the designer distribute scanty resources. Presently, anyone who wants to consider optimum structural design must become familiar with recent synthesis approaches as well as with accepted analysis procedures. 5 Magne1's Diagram Optimum Pair P-e P Initial prestressing torce; e Eccentricity of cable; e*- maximum cable eccentricity; a).b) minimum 1/P; c).d) maximum 1/P. Figure 1.1 Implicit optimization. 6 Methods In the last twenty years researchers have made considerable advances in developing techniques of optimum design. Research and exploration of these methods were mainly developed in the aeronautical and mechanical industries, where the need for more economical and efficient final products was extremely important. More recently, with the availability of increasing computer capabilities, civil engineering researchers and designers have increased their participation in structural optimization following the lines defined by the other engineering disciplines. Optimization methods are, nevertheless, common to these different engineering design areas and are mainly divided in two groups. These are commonly known by the names Optimality Criteria and Mathematical Programming (6). Another area in structural optimization researched by a few scientists is based on duality theory concepts, and is an attempt to unify the two basic methodologies (7). Optimality Criteria methods are based on an iterative approach where the conditions for an optimum solution are previously defined. The concept can be used as the basis for the selection of a structure with minimum volume. This methodology derives from the extreme principles of structural mechanics and has been limited to simple structural forms and loading conditions. The formulation can be mathematically expressed as follows: 7 2k+l =
where x is the vector of design variables, uk+1 is an
recurrence relation. Estimation of the lagrangian where
distribution.
where x is the vector of design variables, uk+1 is an
recurrence relation. Estimation of the lagrangian where
distribution. |