Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo Niko Kristanič Limit State Design Using Exact Sensitivity Analysis and Shape Optimization DOCTORAL THESIS Sinteza konstrukcij z uporabo točne občutljivostne analize in optimizacije oblike v nelinearnem področju DOKTORSKA DISERTACIJA Ljubljana, 2008 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. VII Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. ERRATA Page Line Error Correction Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. IX Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. BIBLIOGRAPHIC-DOCUMENTALISTIC INFORMATION UDC: 624.04(043.3) Author: Niko Kristanič Supervisor: Assoc. prof. dr. Jože Korelc Title: Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Notes: 148 p., 5 tab., 72 fig., 42 eq. Keywords: shape optimization, limit load, worst initial imperfection, symbolic approach. Abstract: Optimization has become an important tool in engineering activities because it represents a systematic method to improve design with respect to certain criteria. Within the thesis a numeric-symbolic approach to limit load shape optimization is studied which enables the use of an optimization algorithm as an ultimate state design tool. Shape is parameterized symbolically using a general computer algebra system. Therefore the design velocity filed can be computed analytically and an exact sensitivity analysis can be carried out. Accurate sensitivity information is of crucial importance for proper gradient shape optimization. When analyzing imperfection sensitive structures it turns out that the choice of the shape and size of initial imperfections has a major influence on the response of the structure and its ultimate state. Further on, shape optimization applied on the perfect mathematical model can lead to non-optimal results, e.g. a very light structure but very sensitive to buckling. While imperfections are not known in advance, a method for direct determination of the most unfavorable imperfection of structures by means of ultimate limit states was developed. The method is implemented as an internal and separate optimization algorithm within the global shape optimization process. Full geometrical and material nonlinearity is considered throughout the global optimization process consistently, resulting in efficient and robust, ultimate limit load structure design algorithm. The numerical examples indicate that the use of a symbolic-numeric system for gradient shape optimization combined with the use of the most unfavorable imperfections can represent a superior alternative to conventional ultimate limit state design. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. XI Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. BIBILIOGRAFSKO-DOKUMENTACIJSKA STRAN IN RAZŠIRJEN IZVLEČEK UDK: 624.04(043.3) Avtor: Niko Kristanič Mentor: izr. prof. dr. Jože Korelc Naslov: Sinteza konstrukcij z uporabo točne občutljivostne analize in optimizacije oblike v nelinearnem področju. Obseg in oprema: 148 str., 5 pregl., 72 sl., 42 en. Ključne besede: oblika konstrukcij, optimizacija oblike, nosilnost konstrukcij, mejna obtežba, najbolj neugodne začetne nepopolnosti, vpliv na stabilnost, simbolni pristop. Izvleček. Optimizacija postaja vedno bolj pomembno orodje v inženirski praksi saj predstavlja sistematično metodo izboljšanja izdelkov glede na dane kriterije. V okviru disertacije je predstavljen simbolno-numerični pristop k optimizaciji oblike konstrukcij v mejnem stanju nosilnosti. Pristop omogoča uporabo optimizacijskega algoritma kot orodje za projektiranje konstrukcij. Oblika konstrukcije je parametrizirana simbolno s pomočjo sistema za splošno računalniško algebro, ki s pomočjo neposrednega odvajanja omogoča analitičen izračun polja začetnih občutljivosti. Posledično je možno izvesti natančen izračun občutljivosti odziva, kar je ključnega pomena, saj so točne občutljivosti pogoj za uspešno uporabo gradientnih metod optimizacije oblike. Kadar obravnavamo konstrukcije, občutljive na spremembo začetne geometrije, se izkaže, da ima izbira oblike in velikosti začetnih nepopolnosti velik vpliv na odziv konstrukcije in njeno mejno stanje. Poleg tega uporaba idealne oblike konstrukcije lahko privede do nestabilnosti optimizacijskih algoritmov ali do neoptimalnih rezultatov, na primer izjemno lahkih konstrukcij, ki so močno občutljive na nepopolnosti. Začetne nepopolnosti niso znane v naprej, zato je v okviru disertacije bila razvita metoda za določitev najbolj neugodne začetne nepopolnosti v smislu mejnega stanja konstrukcij. Metoda je implementirana kot ugnezden optimizacijski algoritem v okviru globalne optimizacije oblike. Skozi celoten proces optimizacije oblike je uporabljen polno nelinearen pristop, ki omogoča učinkovito in robustno sintezo konstrukcij. Prikazani primeri prikazujejo uporabnost metode in nakazujejo, da uporaba simbolno numeričnega okolja za gradientno optimizacijo oblike v povezavi z metodo določitve najbolj neugodnih začetnih nepopolnosti predstavlja napredno alternativo klasičnemu projektiranju konstrukcij. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. XIII Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Zahvala Delo, ki je bilo opravljeno v okviru disertacije, je bilo delno financirano s strani Ministrstva za visoko šolstvo, znanost in tehnologijo, v okviru programa mladih raziskovalcev. Iskreno se zahvaljujem svojemu mentorju izr. prof. dr. Jožetu Korelcu za vso pomoč in vodenje skozi znanstveno delo. Prof. dr. Darku Begu se zahvaljujem za strokovno pomoč in koristne nasvete. Zahvala gre tudi vsem sodelavcem, ki so poskrbeli za konstruktivno debato in prijetno vzdušje. Posebna zahvala je namenjena mojim staršem, ki so me vso dolgo pot podpirali in vzpodbujali ter moji Mateji, ki si mi vselej stala ob strani. Acknowledgement I am very grateful for the help received from Krzysztof Wisniewski, Institute of Fundamental Technological Research, Warszawa, Poland who kindly allowed his symbolic description of shell finite elements to be used within the presented research. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. XV Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Contents 1 Introduction ...................................................................................................... 1 1.1 Structural limit load design optimization ............................................................ 1 1.2 Background of work and state of the art ............................................................ 4 1.2.1 Limit Load Shape Optimization ................................................................................................ 4 1.2.2 Sensitivity analysis .................................................................................................................... 6 1.2.3 Imperfections ............................................................................................................................. 9 1.3 Motivation and Objectives ................................................................................ 13 1.4 Methodology ...................................................................................................... 13 1.5 Thesis Outline ................................................................................................... 15 2 Finite element modeling and symbolic approach ........................................ 17 2.1 Automation of FEM .......................................................................................... 18 2.2 Automatic differentiation .................................................................................. 20 2.2.1 Principles of automatic differentiation .................................................................................... 20 2.2.2 Automatic differentiation and FEM ........................................................................................ 21 2.3 Hybrid symbolic-numerical approach ................................................................ 25 2.3.1 Typical example of automatic code generation procedure ...................................................... 28 2.4 Abstract symbolic formulations in computational mechanics ............................ 29 2.5 Symbolic-numerical environment AceFEM ....................................................... 31 3 Direct and Sensitivity limit load analysis .................................................. 33 3.1 Definition of ultimate states .............................................................................. 33 3.2 Direct Analysis .................................................................................................. 34 3.3 Sensitivity Analysis ........................................................................................... 37 3.3.1 The Analytical Design Velocity Field ..................................................................................... 39 3.3.2 Exact sensitivity analysis ........................................................................................................ 44 3.4 Symbolical formulation of general finite strain plasticity .................................. 45 3.5 Finite element models ....................................................................................... 49 3.6 Example ............................................................................................................ 51 4 Initial imperfections ....................................................................................... 55 4.1 Introduction ...................................................................................................... 55 4.2 Optimization method for the determination of the most unfavorable initial imperfection ............................................................................................ 56 4.2.1 Representation of imperfections .............................................................................................. 56 XVI Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 4.2.2 Description of the algorithm .................................................................................................... 57 4.2.3 Formulation of the optimization problem ............................................................................... 60 4.2.4 Direct and sensitivity analysis ................................................................................................. 62 4.3 Numerical examples (Test Problems and Results) ............................................. 63 4.3.1 Elasto-plastic cantilever structure ........................................................................................... 63 4.3.2 Thin-walled T beam ................................................................................................................ 68 4.3.3 Thin-walled I beam ................................................................................................................. 73 4.3.4 Thin-walled Cylinder ............................................................................................................... 76 4.4 Partial conclusions ............................................................................................. 80 5 Limit load design Shape Optimization ........................................................... 83 5.1 Introduction ....................................................................................................... 83 5.2 Shape Optimization ........................................................................................... 84 5.2.1 Evolutionary and Optimality Criteria approaches .................................................................. 84 5.2.2 Mathematical programming .................................................................................................... 85 5.3 Gradient based Shape Optimization combined with imperfection analysis .............................................................................................................. 88 5.4 Numerical examples ........................................................................................... 92 5.4.1 2D cantilever shape optimization ............................................................................................ 92 5.4.2 3D H cross-section thin-walled cantilever structure ................................................................ 95 5.4.3 3D single storey steel building ............................................................................................... 111 6 Final Conclusions ........................................................................................... 121 7 Povzetek .......................................................................................................... 125 References ............................................................................................................. 143 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. XVII Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Index of figures Fig. 1: Flow charts of conventional structural design (a) and optimum structural design (b). . .......................... 2 Fig. 2: Schematic representation of the overall algorithm for limit load shape optimization. . .......................... 14 Fig. 3: Multi?language and multi?environment FE code generation .................................................................. 27 Fig. 4: Typical AceGen input. . ............................................................................................................................ 29 Fig. 5: Typical automatically generated subroutine in FORTRAN and C language. . ........................................... 29 Fig. 6: Definition of ultimate states (EN 1993 1?5 2004) ................................................................................... 34 Fig. 7: General formulation of direct analysis of transient coupled nonlinear problems. . ................................ 36 Fig. 8: General formulation of shape sensitivity analysis of transient coupled nonlinear problems ................... 38 Fig. 9: Geometry of cantilever with marked node numbers and shape parameters. . ........................................... 40 Fig. 10: AceFEM input data for the structure illustrated in Fig. 9. . .................................................................... 40 Fig. 11: AceFEM node coordinates for linear boundary shape approximation. . ................................................... 41 Fig. 12: AceFEM node coordinates for quadratic boundary shape interpolation. . .............................................. 41 Fig. 13: Design velocity field by symbolical derivation of FE node coordinates for the case of linear interpolation between shape parameters. . .............................................................................................. 42 Fig. 14: Design velocity field by symbolical derivation of FE node coordinates for the case of quadratic interpolation between shape parameters. . .............................................................................................. 42 Fig. 15: Graphical representation of the y component of the design velocity field for the case of linear interpolation between shape parameters. . .............................................................................................. 43 Fig. 16: Graphical representation of the y component of the design velocity field for the case of quadratic interpolation between shape parameters. . .............................................................................................. 44 Fig. 17: Flowchart of shape sensitivity analysis by dual symbolic?numeric FE environment .................................. 45 Fig. 18: Algorithm for the abstract symbolic description of elasto?plastic problems. . ........................................ 47 Fig. 19: Summary of the finite strain plasticity equations. . ................................................................................. 48 Fig. 20: Finite element model for the one story building. . ................................................................................... 52 Fig. 21: Example of a nodal point coordinate in symbolic form. . ......................................................................... 53 Fig. 22: Comparison between analytic and FD method. . ...................................................................................... 53 Fig. 23: Graphical sensitivity representation of the vertical displacement with respect to the roof angle. . ......... 54 Fig. 24: Flowchart of the method for the determination of the most unfavorable initial imperfection ............... 59 Fig. 25: Geometry and loading (a) and the logically most unfavorable shape without considering technological constraints (b) for the cantilever beam example. . ............................................................ 63 Fig. 26: The shape base for the cantilever beam example. . .................................................................................... 64 Fig. 27: Convergence to the most unfavorable shape by increasing the number of considered shapes in the shape base (scale factor fS = 10). . ............................................................................................................ 66 Fig. 28: Geometry, supporting and loading conditions for the thin?walled T beam example. . .............................. 68 Fig. 29: Maximal allowable deviation from the perfect geometry used in the evaluation. . .................................. 69 Fig. 30: The calculated ultimate load of the T?beam considering the evaluated most unfavorable initial imperfection varying the number of base shapes. . ................................................................................... 70 Fig. 31: The load displacement curves considering various imperfection shape bases for the T beam example. . ..... 71 Fig. 32: The most unfavorable initial imperfection shape of the T cross?section thin?walled girder (a) and the corresponding deformed shape at collapse (b). . ................................................................................ 71 Fig. 33: Convergence of the global iteration process of finding the most unfavorable imperfection shape for the example with 52 base shapes. . ........................................................................................................... 72 Fig. 34: Equilibrium paths for the T cross?section thin?walled girder. . ............................................................... 73 Fig. 35: Maximal allowable deviation from the perfect geometry used in the evaluation. . .................................. 74 Fig. 36: Ultimate load deformation curves for different initial imperfect geometries with the same amplitude for the I beam example. . ......................................................................................................... 75 Fig. 37: Initial imperfect geometries used in analyses according to Fig. 36. . ........................................................ 76 Fig. 38: Geometry of the axially and transversely loaded cylinder. . .................................................................... 77 XVIII Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 39: Most unfavorable initial imperfection for axially and transversely loaded cylinder considering the imperfection amplitude R/165. Scale factor fS=10. . ................................................................................ 78 Fig. 40: Load displacement curve considering the most unfavorable initial imperfection ..................................... 78 Fig. 41: Deformation state at limit load. . ........................................................................................................... 79 Fig. 42: Fold line of the axially and transversely loaded cylinder ........................................................................ 79 Fig. 43: Optimization using symbolic?numeric environment ................................................................................. 89 Fig. 44: Optimization loop using symbolic?numeric environment. . ....................................................................... 90 Fig. 45: Initial shape (a), Optimal shape (b). . ........................................................................................................ 93 Fig. 46: Convergence of the optimization process for two shape parameters. . ..................................................... 94 Fig. 47: Structural response of different shapes and the limit load optimal shape. . ............................................. 94 Fig. 48: Initial geometry of H cross?section cantilever ........................................................................................ 96 Fig. 49: H cross?section cantilever shape parameters. . ........................................................................................ 97 Fig. 50: The most unfavorable imperfection for the initial (a) and optimal shape (b) (Scale factor=30, Shell thickness A) ................................................................................................................................ ............ 99 Fig. 51: Limit load shape optimization process. . ................................................................................................. 101 Fig. 52: Initial and optimum shape of H cross?section cantilever geometry for case A shell thickness. . .............. 102 Fig. 53: Mises stress at limit state for optimal H cross?section cantilever shape (undeformed). . ....................... 103 Fig. 54: Deformation of H cross?section cantilever at limit state (Scale Factor = 1) with Mises stress plotted. . .............................................................................................................................. ................ 103 Fig. 55: Mises stress for the corresponding load?deformation curve plotted in Fig. 51 ..................................... 104 Fig. 56: The most unfavorable imperfection for the initial (a) and optimal shape (b) (Scale factor=30, Shell thickness B). . .............................................................................................................................. .......... 105 Fig. 57: Initial and optimum shape of H cross?section cantilever for case B shell thickness. . .............................. 106 Fig. 58: Mises stress at limit state for optimal H cross?section cantilever shape (undeformed). . ....................... 107 Fig. 59: Deformation of H cross?section cantilever at limit state (Scale Factor = 1) with Mises stress plotted. . .............................................................................................................................. ................ 107 Fig. 60: Mises stress for the corresponding load?deformation curve plotted in Fig. 61 ..................................... 108 Fig. 61: Load?deformation curve for optimized shape in Case A and Case B. . ...................................................... 109 Fig. 62: Initial geometry of single storey steel building. . ................................................................................... 111 Fig. 63: Loading conditions and shape parameterization for the inner frame ..................................................... 112 Fig. 64: Loading conditions and shape parameterization for the outer frame. . .................................................. 113 Fig. 65: Optimum shape of the inner frame of the steel structure with the initial geometry in the final optimization run. . .............................................................................................................................. .. 114 Fig. 66: Optimum shape of the outer frame of the steel structure with the initial geometry in the final optimization run. . .............................................................................................................................. .. 114 Fig. 67: Deformation of the inner frame at limit state (Scale Factor = 10). . ...................................................... 115 Fig. 68: Deformation of the outer frame at limit state (Scale Factor = 10). . ..................................................... 116 Fig. 69: Mises stress for inner frame for the corresponding load?deformation curve plotted in Fig. 71. . .......... 117 Fig. 70: Mises stress for outer frame for the corresponding load?deformation curve plotted in Fig. 71. . ......... 118 Fig. 71: Load?deformation curves for the inner and the outer frame with the optimized shape .......................... 119 Fig. 72: Optimized shape of single story steel building. . ..................................................................................... 120 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. XIX Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Index of Tables Table 1: Automatic differentiation exceptions ..................................................................................................... 25 Table 2: Residual form of equations for mechanical problems. . .......................................................................... 34 Table 3: Strain energies for the use in the Automatic differentiation exceptions. . ................................................ 50 Table 4: Ultimate load factors for various initial imperfection shapes for the cantilever beam example. . ............ 67 Table 5: Optimization parameter values for the optimized structure. . ............................................................... 120 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 1 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 1 Introduction It is in the nature of mankind and nature itself to strive after the achievement of a goal with the least possible amount of effort. Technical development has always required incessant improvement of solutions. This has been achieved by systematic improvement of initial design, by redesigning, implementing new knowledge and learning from past mistakes. In this manner the initial form of design optimization has been a sort of trial and error process which in the end has given better and better solutions. Although a very time consuming process, the basic idea is adopted by most modern design optimization strategies. Nowadays the technical development is accompanied with a growing competition which demands lower design and production cost, higher quality, less energy consumption, environment friendly and recycling ability design and design of products with the required aesthetic value. To meet the increasing necessity for competitive products, a scientific approach is needed. The latest knowledge from the field of computational mechanics, sensitivity analysis and optimization has to be used. Within the work presented by the thesis a limit load shape optimization method including worst imperfection evaluation was developed using the latest symbolic-numerical environment technology. It is difficult to fulfill all the requirements and to claim that one of the evaluated potential optimum designs was the best possible. The decision is left to the design engineer to define the performance criteria and restrictions. In the present work the optimum design is therefore defined as the one that maximizes the chosen performance criteria and satisfies al the constraints given by the engineer. 1.1 Structural limit load design optimization It is important to recognize the difference between structural analysis and structural design. The analysis problem is concerned with determining the behavior of an existing structural system with a known design, while structural design is a task where the design is varied to meet performance requirements. Until recently conventional structural design has mostly been used for designing structures. Conventional structural design is a trial and error procedure and depends on the designer’s intuition, experience and skill. It can be a difficult challenge for an 2 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. engineer to achieve efficiency, cost-effectiveness and an overall integrity of a designed structure. Sometimes this approach can lead to erroneous results and uneconomical structures when dealing with complex structural systems. The growing industrial competition is forcing engineers to consider different approaches, economical and better design. Design optimization is one of such approaches. The difference between conventional design and design optimization called synthesis is clearly illustrated in Fig. 1. Collection of data to describe structural system Definition of Design variables, Cost function, Constraints to satisfy T Collection of data to describe structural system Estimation of initial design Estimation of initial design Check of performance criteria Change of design based on heuristics/experience Evaluation of performance criteria and check for the constraints Design satisfiesYES convergence criteria ? Change of design using optimization method b) Fig. 1: Flow charts of conventional structural design (a) and optimum structural design (b). Slika 1: Običajni potek projektiranja (a) in projektiranje s pomočjo optimizacije (b) Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 3 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Design optimization process forces the designer to define a set of design variables, an objective function to be optimized, and the constraint functions to be taken into account. Despite the extra work the path to better design is certainly shorter than within conventional design. Although the optimization approach seems more formal, it can substantially benefit from the designer’s experience and intuition in formulating the problem, choosing proper parameterization and identifying the critical constraints. Thus, the best approach would be an optimum design process aided by the designer’s interaction. Proper mathematical formulation of the design optimization problem is a key to good solutions. The optimization loop within the optimum structural design algorithm (Fig. 1b) can be divided into three characteristic steps: - Evaluating the performance measure of the structure by using the current design variables (direct analysis). - Evaluating sensitivity of the design to changes of all design variables, where sensitivities are the gradients of the objective and constraint functions used by the optimization algorithm (design sensitivity analysis). - Updating the design variables using sensitivity information in a way that improves the objective function (optimization). Limit load structural optimization design demands a complex interaction of different approaches, used algorithms and methods. When designing structures for the ultimate state, the limit load of the designed structure has to be known and is usually given by technical standards and codes. In the optimization algorithm the limit load presents a constraint that has to be fulfilled. When dealing with imperfection sensible structures, e.g. thin walled structures, imperfections have to be taken into account. A method for automatic evaluation of the most unfavorable initial imperfection is presented. In the optimization process the design is changed by varying design variables represented by shape parameters. When using gradient based shape optimization approaches, the hardest problem is to evaluate accurate sensitivities with respect to shape parameters for the gradient information. With the use of a symbolic-numeric system, exact sensitivities can be evaluated by using an analytically calculated velocity field. Arbitrary symbolic shape parameterization can be used with no limitation. 4 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Direct and sensitivity analyses within shape optimization and the evaluation of the most unfavorable initial imperfections are performed to the ultimate limit state of the structure. Within conventional shape optimization methods usually the minimum weight is sought for a certain stress distribution state, the state at first yielding point or buckling load, but no real limit state is considered as it is difficult to define and to evaluate. Within the presented approach the optimal shape is sought for the ultimate limit state of the structure resulting in a robust optimal design which has to satisfy all design criteria. The combination of a symbolic–numeric system and an algorithm for the automatic determination of most unfavorable initial imperfections gives rise to an effective structural design shape limit load optimization method appropriate to be used as an efficient designing tool. 1.2 Background of work and state of the art 1.2.1 Limit Load Shape Optimization Historically, structural optimization used to be a global process of progress in structural design, mainly based on experience and experiments. The engineering skills have improved gradually and the results of designs have became more and more optimal. Many results of current basic structural optimization problems solved by contemporary optimization methods are in close relation to the results gained through the historical design development. The most important point in structural optimization is to define the relation between geometry and the internal flow of forces. Throughout the historical design development this was the focus of experiments and intuitive design. The first analytical works of structural optimization appeared in 17th and 18th century. Important contributors were Galilei (1638; Discorsi e Dimostrazioni Matematiche, intorno a due nuove scienze), Bernoulli (1687; The brachistochrone problem) and Lagrange (1770; Miscellanea Taurinensia). Mainly the first analytical works were interested in particular cases of optimal sections of beams and columns and no general application was developed. In the 19th and 20th century numerical methods constantly developed further. The key to the modern design was the development of computers, structural analysis methods and mathematical programming. The first to integrate numerical methods into optimization techniques was Schmit (1960) who proposed the concept of structural synthesis. The basic idea of structural synthesis is to integrate finite Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 5 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. element structural analysis into nonlinear mathematical programming methods, which results in an automated optimum design process. The field of structural synthesis developed further in the next decades. Several other disciplines had to be involved, i.e. structural sensitivity evaluation, CAD preprocessors, and control systems for managing different phases, which interact with each other in the solving process. The first optimization problems considered within structural synthesis were mainly concerned with size optimization of discrete structures. In size optimization, the geometry of the mathematical model is known (e.g. shape, topology) and the characteristics of the mathematical model elements have to be determined (e.g. cross section, material). An example may be a truss structure modeled with truss finite elements where the cross-sections of the trusses have to be determined so that the overall weight of the structure is the smallest. In shape optimization, the geometry of the mathematical model itself is the subject under consideration. The crucial difference between sizing and shape optimization is related to how design variables affect the analysis rather than to the physical optimization problem itself. In size optimization the design variables are related to the properties of the finite elements where in shape optimization the design variables are related to the positions of the finite element nodes and therefore directly affect the implementation within structural analysis. The first attempts of shape optimization were therefore performed on discrete modeled structural systems. Zienkiewicz and Campbell (1973) considered finite element nodes coordinates as design variables which later turned out not to be ideal for solving optimization algorithms. Nevertheless they were the first to set up a shape optimization problem in a general form. The use of optimization in design began to strengthen in the 1980’s together with the intensive development of numerical analysis methods and nonlinear mathematical programming and with expanding computer capabilities. The start of the development of new shape optimization methods gave a clear indication that numerous new problems would arise, which were not known in classical discrete optimization (Haftka, Grandhi 1986). Nowadays discrete optimization is well developed and already integrated in everyday structural analysis programs. Shape optimization, on the other hand, remains the subject of continuous scientific research (see e.g. Bletzinger, Ramm 2001, Camprubi, et al. 2004, Choi, Kim 2005a, Choi, Kim 2005b, Maute, et al. 1999, Ramm, Mehlhorn 6 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 1991, Schwarz, et al. 2001, Uysal, et al. 2007). In continuous shape optimization the design parameterization is not automatically given by the structural model and has to be explicitly defined by the designer in relation to the structural analysis model or to an underlying geometrical model. Within most modern approaches the need of general shape parameterization has evolved in an implicit relation of design variables with respect to the positions of the finite element nodes. As a consequence a complex interaction of the finite element method and sensitivity analysis has to be taken into account. An essential part of the procedure is to evaluate the design velocity field which is defined as the derivative of finite element node coordinates with respect to design variables. Within the presented work a numeric-symbolic approach to optimization is studied which enables the use of arbitrary symbolic shape parameterization and evaluation of an analytical design velocity field. As a consequence, an exact sensitivity analysis can be carried out. Accurate sensitivity information is of crucial importance for proper gradient shape optimization. The sensitivity analysis and the evaluation of the design velocity field will be further addressed in the next section. 1.2.2 Sensitivity analysis To overcome the difficulties of evaluating sensitivity within shape optimization, numerous methods have been developed in the past (see review by van Keulen, et al. 2005). The four broad categories of methods in common use for obtaining the derivatives of performance measures with respect to structural parameters are: a) Overall finite differences b) Discrete derivatives c) Continuum derivatives d) Computational or automatic differentiation The choice of the method is particularly important in gradient shape optimization where the shape design variables change the discretization of the discretizied problem, e.g. finite element mesh. All methods except for the finite differences can be implemented using direct or adjoint approach (called also the reverse mode of automatic differentiation explained in Section 2.2). In the direct approach, the derivatives of the entire structural response are evaluated. The sensitivities of performance measures can then be obtained from the chain rule of differentiation. In the adjoint approach an adjoint problem which depends on the performance measure is defined. The sensitivities of Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 7 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. performance measures can then be obtained using the structural and adjoint responses. This approach is of most advantage when the problem consists of many design variables and few performance measures while not all system response sensitivities are required. Overall finite differences (see e.g. Hörnlein H.R.E.M. 2000, Oral 1996) present the simplest method. The method, also called global finite difference, is based on repeated evaluation of structural analysis code and the use of a finite difference formula to obtain the derivative. Forward, backward and central differences can be used. Higher order difference formulae are very rare. Finite difference derivatives can suffer from truncation errors with large step sizes and also from errors when the step size is too small. The computational time is high due to the repeated code evaluation. Global finite differences become very useful when using commercial structural analysis programs where the analysis code is in a form of a black box with no ability to solve the sensitivity problem. Continuum derivatives are obtained by differentiating the governing continuum equations. Most commonly these consist of partial differential equations or an integral form, for example, derived from the principle of virtual work. The differentiation leads to a set of continuum sensitivity equations that are then solved numerically. The same discretization as for the original structural response can be used. For shape sensitivities, the two main approaches for continuum derivatives are the material derivative approach (see e.g. Saliba, et al. 2005) and the control volume approach (see e.g. Arora, et al. 1992). The advantage of these methods is the possibility of different meshes for response and sensitivity analysis. The most widely used methods are the discrete derivatives. While the continuum sensitivity equations are derived by differentiating the governing continuum equations with respect to the design variables and are then subsequently discretized, for discrete derivatives this order is reversed. The advances of these methods are low computational cost and high consistency. They can be separated into analytical and semi-analytical. The analytical methods use analytical derivatives on the global level as well on the finite element level (see e.g. Maute, et al. 2000). The analytical differentiation process may become tedious. This holds true especially for shape design variables, therefore symbolic computing software can be applied (Korelc 2002, Korelc 2007a, b) which often features the automatic generation of the source code. This code has to be integrated in the used software. Additional procedures must be implemented for each finite element used within the structural analysis. The procedure must account for all 8 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. possible design variables, i.e. size, material and shape design variables, as the actual code depends on the type of the design variable. Shape design variables are the most complex ones to implement. The semi-analytical methods (see e.g. Hörnlein H.R.E.M. 2000, Oral 1996) use finite difference schemes on the finite element level, while the dicretized governing equations on the structure level are differentiated analytically. The main reason for the simplification is the significant implementation effort needed for the evaluation of the finite element sensitivity pseudo-load vector. Therefore approximations in the form of finite differences are frequently accepted for the pseudo-load vector, which reduces the effort. Disadvantages of these methods are common for all finite difference methods, the dependency with respect to the perturbation size, which has a pronounced effect on both consistency and efficiency. The highest consistency is proved by the methods which use automatic differentiation. Analytical derivation is used on all levels. Even if the finite element program is composed of many complicated subroutines and functions, they are basically a collection of elementary functions. The automatic differentiation method defines the partial derivatives of these elementary functions, and then the derivatives of complicated subroutines and functions are computed using propagation of the partial derivatives and the chain rule of differentiation. Thus, no approximation is introduced. For the highest efficiency of automatically generated codes it is necessary to fulfill specific requirements in order to produce element source codes that are as efficient as manually written codes (see e.g.Korelc 2002). More reference is given in Chapter 2. A major step in performing shape design sensitivity analysis is the evaluation of the design velocity field. The purpose of design velocity field is to characterize the changes of the finite element nodal point coordinates with respect to the changes of arbitrary design parameters. While the design derivatives of the finite element quantities (residual, tangent matrix, etc.) can be constructed by automatic procedures (see e.g.Korelc 2002), this is not true for the design velocity field. Within standard approaches to finite element mesh generation, either with specialized preprocessors or with CAD tools, there exist no explicit relations between the positions of the finite element nodes and the shape design parameter as the choice of the shape parameters is an arbitrary decision of the designer. A number of methods have been proposed in the literature to compute the design velocity field (Choi, Chang 1994). Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 9 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Most frequently used methods for evaluating the design velocity field are the finite difference methods using mesh generators (a), isoparametric mapping methods (b), boundary displacement and fictitious load methods (c), and methods which combine isoparametric mapping and boundary displacement methods (d). All these methods are based on different numerical approximations. The most sophisticated methods in use are the design model concept methods basically fitting into the isoparametric approach section (see e.g. Kegl 2000, Samareh 1999). A great review was done by Haftka and Grandhi (Haftka, Grandhi 1986). The basic idea of the design element approach relies on the assumption that the geometrical data of the structure are not a simple set of constants defining directly the finite element mesh. Instead of that, the structural data are extracted from a set of geometrical objects called the design elements. The shape of the design elements is connected to the finite element mesh and varied with a few shape parameters using e.g. Bezier curves or polynomials. Although an analytical design velocity field can be evaluated, these methods are limited in the choice of design shape parameters. In complex structural systems the shape has to be composed of many design elements which limit the general applicability. For example, no shape parameters can be defined for global structure dimensions. To overcome the difficulties of the design element approach, an arbitrary symbolic parameterization is used in the current work using a symbolic-numeric system (Korelc 2007a, b) further addressed in Chapter 2. Most methods for evaluating design sensitivity are limited to the use of linear or simple nonlinear material models. The use of complicated nonlinear material models makes the evaluation of design sensitivity very difficult. The attempts to develop methods that address this field are therefore rare. Nonlinear limit load design sensitivity analysis used in the present work for the use in gradient based limit load shape optimization demands complex interaction between shape parameterization and finite element code. Within the development of integration of commercial finite element analysis programs and new design approaches using optimization, every effort was made (see e.g. Chang, et al. 1995, OptiStruct 2008) although to the author’s knowledge an efficient tool for nonlinear limit load design optimization is not available. 1.2.3 Imperfections Limit load design optimization requires an accurate evaluation of the limit load of a structure. In order to evaluate the limit load correctly, all the relevant phenomena 10 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. have to be considered e.g. geometrical and material imperfections, load imperfections, residual stresses and strains, damage, etc. It is now well known that geometrical, structural, material and load imperfections play a crucial role in the load carrying behavior, especially of thin walled structures. Most structural imperfections (residual stresses, geometrical and welding imperfections, etc.) are not known in advance. To include the imperfections in a structural analysis, they have to be assumed. Technical standards therefore suggest different approaches to include imperfections on an empirical basis. A convenient way to include all relevant imperfections (i.e. geometrical, structural and material imperfections) is to consider equivalent geometrical imperfections. In this way the geometrical imperfections are augmented by the influence of other relevant imperfections to produce the same effect on the load carrying behavior of a structure. The idea to find imperfections that will cause the structure to fail at the lowest possible load is as old as the ascertainment of the crucial role of imperfections itself. The known discrepancy between theoretical results and experimentally obtained values for ultimate loads of structures can be reduced by properly including imperfections in an analysis. In order to achieve this for a general structure, a series of full geometrical and material nonlinear analyses up to the ultimate limit state need to be performed for a large range of possible imperfections, varying both their shapes and amplitudes. The computational cost involved discourages this kind of direct approach and has been the motivation for the development of computationally less expensive methods. Numerous approaches for analyzing the effect of imperfections on the response of structures have been proposed. Among them one can basically distinguish between those which are derived from the hypothesis that it is possible to obtain a sufficiently accurate structural response for an imperfect structure from the properties of a perfect structure (“perturbation approach”) and those that obtain the structural response by analyzing the imperfect structure itself (“direct approach”). The review of different approaches accompanied with the impact on modern design procedures of engineering structures can be found in (Schmidt 2000). The origin of the perturbation approach goes to the pioneering work on stability of shells done by (Koiter 1945).The perturbation approach applies to structures showing bifurcation phenomena along their natural equilibrium path and is based on asymptotic descriptions of the initial post critical behavior. Although simple for implementation and numerically efficient, the original Koiters theory does not Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 11 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. account for several effects that can significantly lower the ultimate load of the structure, such as: a) nonlinear natural equilibrium path, b) nonlinear material behavior, c) buckling mode interactions, d) consideration of realistic technological constraints on the shape and amplitude of the imperfections, e) large postcritical deflections, f) large imperfections. The interaction of several buckling modes can be assessed by the “minimum path theory” of (Ho 1974, Lanzo 2000, Lanzo, Garcea 1996). An overview of some other methods was given by (Godoy 2000). There are numerous difficulties obtaining nonlinear post critical behavior. Even with modern methods such as arc-length schemes, branch switching procedures, direct computation of stability points and stabilization techniques, it is sometimes impossible to overcome problems like secondary bifurcations, coincidental or clustered singularities and post critical paths, that cross each other, to get all the possible hypothetic equilibrium paths or at least the minimum one. Another approach is the “minimum perturbation energy concept”, which has recently been applied also to dynamic stability problems (Dinkler, Pontow 2006, Ewert, et al. 2006). The basic idea of the method is to lower the buckling load of a perfect structure by introducing a certain amount of energy into the system, which causes a snap through to the post-buckling path, or to a secondary path in dynamic problems. Common to all perturbation methods is that they become exceedingly complicated for implementation and numerically inefficient when phenomena such as mode interaction or plasticity are included. This has been an inspiration for the research on the second branch of methods where the structure is analyzed using a full nonlinear analysis that by definition includes phenomena (a), (b), (e) and (f). The buckling or limit load of a theoretically perfect structure is lowered by introducing imperfections directly into the geometry of the structure. There exists a huge amount of uncertainties in the determination of the shape and amplitude of real imperfections because of the nature of the production processes. Therefore it would be natural to use probabilistic approaches where imperfections are introduced as random variables with a certain distribution (Elishakoff 2000, Papadopoulos, Papadrakakis 2005, Schenk, Schueller 2003). However, these methods rely upon very scarce data banks of 12 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. measured imperfections and are therefore useful for a certain assortment of applications. An alternative approach to include imperfections in an analysis is the concept of the “definitely worst” imperfection (Deml, Wunderlich 1997, El Damatty, Nassef 2001, Song, et al. 2004, Wunderlich, Albertin 2000, 2002). Within the “definitely worst” imperfection concept the shape of imperfections is searched that would lead to the lowest ultimate load of the structure. The shape of imperfections is additionally bounded by the given imperfection amplitude. Several variants of the procedure are possible and discussed (Schmidt 2000). Deml and Wunderlich (Deml, Wunderlich 1997) introduced an elaborate method to obtain the “definitely worst” imperfection. In their case imperfections are treated as additional degrees of freedom. The approach results in an extended system of equations composed of equilibrium equations, equations for direct computation of stability points (Wriggers, Simo 1990), condition equations for the worst imperfections shape and constraint equations that limit the amplitude of the imperfection. The system has to be solved simultaneously for the equilibrium state, the worst imperfection shape and the corresponding ultimate load. The result is the “definitely worst” imperfection within the considered amplitude. Such procedure is rather time consuming and is therefore useful for problems of small order (Wunderlich, Albertin 2000, 2002). Within the approach of Deml and Wunderlich (Deml, Wunderlich 1997) the condition equations for the worst imperfection are based on Koiters asymptotic theory, thus limiting the approach to small imperfections and linear fundamental paths. The determination of the “definitely worst” imperfection can also be formulated as a nonlinear optimization problem solved by one of the well-known nonlinear optimization methods. This is the most general approach that includes all the relevant phenomena, but is for the same reason also the most computationally expensive one. For example, a genetic optimization algorithm for obtaining the worst imperfection of shell structures was used by (El Damatty, Nassef 2001). Apart from the obvious disadvantage of the optimization approach, i.e. the need for a potentially large number of full nonlinear analyses, there are also advantages. Since the bifurcation-type instability always represents some type of a symmetry breakdown, bifurcations can be avoided if the symmetry is deliberately destroyed by introducing imperfections (Bažant, Cedolin 2003). By limiting the analysis to limit points, one avoids tedious procedures involved in proper determination, classification and sensitivity analysis of bifurcation points. Another certain advantage of the direct Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 13 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. approach is that only stable equilibrium states have to be considered and no hypothetic or unstable states need to be relied upon. It is doubtful that an approach requiring a large number of fully nonlinear analyses could be used in everyday engineering design in foreseeable future. In the present work a computationally less expensive optimization method is developed that would still retain the generality of the optimization based “definitely worst” imperfection approach. The method is further addressed in Chapter 4. 1.3 Motivation and Objectives The integration of optimization methods into engineering design is a complex task. As pointed out in the previous sections, there are numerous difficulties to achieve this for a general nonlinear case. To facilitate the use of synthesis to design engineering structures, an effective optimization method has to be used considering full nonlinearity with the use of automatic definition of proper initial imperfections. The general objective of this work is to develop a finite element based limit load shape optimization technique which can be used for the ultimate limit design of structures. The key aspects to be investigated in this work are: - The use of a symbolic-numeric system for limit load analysis and optimization purposes. Symbolic derivation and automatic code generation of finite elements for direct and sensitivity analysis. - The evaluation of most unfavorable initial imperfections of a structure by means of ultimate limit states. - Analytical evaluation of the design velocity field, using arbitrary shape parameterization, for the use in exact sensitivity analysis. - Development of an efficient gradient based limit load shape optimization method based on all the above components. 1.4 Methodology Limit load shape optimization is performed with state of the art optimization algorithms using the computer algebra software Mathematica (Wolfram 2008) and a symbolic-numeric approach using automatic code generator AceGen (Korelc 2007b) and finite element environment AceFEM (Korelc 2007a). Finite element method is used to model the structure subjected to optimization. Direct and sensitivity analysis 14 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. of structures by using the most unfavorable initial imperfect geometry is performed in AceFEM. All the necessary finite element codes are developed using abstract symbolic description with simultaneous optimization of expressions, automatic differentiation technique, theorem proving and automatic generation of finite element code (Korelc 2002) using AceGen. The use of the numeric-symbolic system offers the possibility one to use an analytical approach. Sensitivity calculation with the use of an analytically evaluated design velocity field is of crucial importance for the convergence of the gradient optimization algorithm, especially when dealing with highly geometrical and material nonlinear problems. The overall algorithm is presented in Fig. 2. Most unfavorable Imperfection evaluation Direct and Sensitivity analysis Optimization ( Mathematica ) Fig. 2: Schematic representation of the overall algorithm for limit load shape optimization. Slika 2: Prikaz celotnega postopka optimizacije oblike konstrukcij v mejnem stanju. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 15 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 1.5 Thesis Outline The thesis is organized in 4 main chapters. - Chapter 2 introduces the symbolic description of mechanical problems, the use of advanced automatic differentiation techniques and the hybrid symbolic-numeric environment. - Chapter 3 is devoted to the structural sensitivity analysis using an analytical design velocity field. The general expressions for the direct and sensitivity analysis which can be used for abstract symbolic description of transient nonlinear coupled systems are given. - Chapter 4 explains the use of imperfections in a limit load structural analysis. A method for the evaluation of the worst imperfect geometry is presented. - Chapter 5 incorporates all the relevant techniques presented in previous chapters to perform a limit load shape design optimization. Examples are given for typical civil engineering structures. - Chapter 6 summarizes the conclusions drawn from present work. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 17 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 2 Finite element modeling and symbolic approach The most effective and widely spread method for solving problems in solid mechanics nowadays is the finite element method (FEM). The method originates from the needs to solve complex, structural analysis problems in civil and aeronautical engineering. Today FE methods are far more advanced and can be used for highly nonlinear direct and sensitivity analyses, inverse modeling and optimization of Multi-field, Multi-scale, Multi-body, Multi-phase and Multi-objective problems. Te purpose of the present work is not to introduce fundamental knowledge of Continuum Mechanics and of Finite Element Methods, as many references can be found on these matters. An introduction to Continuum Mechanics is given by (see e.g. Marsden, Hughes 1994) and (Lemaitre, Chaboche 1990). An overture to the Finite Element Method can be found in (Zienkiewicz, Taylor 2000b) and (Crisfield 1996; Crisfield 1997). With the employment of FEM to increasingly complicated and bounded problems, the implementation of the method has become highly sophisticated. The struggle to automate some of the tasks in the overall processes which starts with the development of theories and ends with a working program is therefore highly appreciated. Some major achievements have been attained in this field of computational mechanics in the last decade. Automation of the finite element method has attracted attention of researches from the fields of mathematics, computer science and computational mechanics, resulting in a variety of approaches and available software tools. The use of advanced software technologies, especially symbolic and algebraic systems, problem solving environments and automatic differentiation tools influences directly how the mechanical problem and corresponding numerical model are postulated and solved, leading to the automation of the finite element method. In order to formulate nonlinear finite elements symbolically in a general but simple way, a clear mathematical formulation is needed at the highest abstract level possible. The finite element technology used in the present work is based on a symbolic-numeric approach with a high degree of automation implemented in software 18 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. packages AceGen and AceFEM (Korelc 2007a, b). The basic techniques and methods used within are described in next sections. 2.1 Automation of FEM Automation of FEM is a complex task because of the various transformations, differentiation, matrix operations, and a large number of degrees of freedom involved in the derivation of characteristic FE quantities, which often leads to exponential growth of expressions in space and time (Korelc 2002). The complete FE simulation can be, from the aspect of the automation level, decomposed into the following steps: - formulation of strong form of initial boundary-value problem; transformation of the strong form into weak form or variational functionals; - definition of the domain discretization and approximation of the unknowns and the virtual fields; - derivation and solution of additional algebraic equations defined at the element level (e.g. plastic evolution equations); - derivation of algebraic equations that describe the contribution of one element to the global internal force vector and to the global tangential stiffness matrix; - coding of the derived equations in required compiled language; - generation of finite element mesh and boundary conditions; - solution of global problem; - presentation and analysis of results. Alternatively, one can also start from the free Helmholtz energy of the problem and derive element equations directly as a gradient of the free energy. This approach is especially appealing for the automation due to the numerical efficiency of the solution when the gradient is obtained by the backward mode of automatic differentiation. There are almost countless ways of how a particular problem can be solved by the FE method. If the automation of all nine steps is chosen, then only very specific subset of possible formulations can be covered. On the other hand, the standard discretization is of little use for problems involving coarse mesh, locking phenomena and distorted element shapes where highly problem specific formulations have to be used. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 19 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. The following techniques, which are the result of rapid development in computer science in the last decades, are particularly relevant for the description of a nonlinear finite element model on a high abstract level, while preserving the numerical efficiency: - Symbolic and algebraic computational systems - Automatic differentiation tools - Problem solving environments - Hybrid approaches Computer algebra (CA) systems are tools for the manipulation of mathematical expressions in symbolic form. Widely used CA systems such as Mathematica or Maple have become an integrated computing environment that covers all aspects of computational process, including numerical analysis and graphical presentation of the results. In the case of complex mechanical models, the direct use of CA systems is not possible due to several reasons. For the numerical implementation, CA systems cannot keep up with the run-time efficiency of programming languages such as FORTRAN and C and by no means with highly problem-oriented and efficient numerical environments used for finite element analysis. However, CA systems can be used for the automatic derivation of appropriate formulae and generation of numerical codes. The FE method is within the general CA systems usually implemented as an additional package or toolbox such as AceFEM (Korelc 2007a) for Mathematica used in the work covered by this thesis. The major limitation of symbolic systems, when applied to complex engineering problems, is an uncontrollable growth of expressions and consequently redundant operations and inefficient codes. (see Korelc 1997, Korelc 2002) This is especially problematic when CA systems are used to derive formulae needed in numerical procedures such as finite element method where the numerical efficiency of the derived formulae and the generated code are of utmost priority. The automation level of FE method can be greatly increased by combining several approaches and tools. A hybrid symbolic numeric system is used within the thesis described in Section 2.3 20 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 2.2 Automatic differentiation Differentiation is an arithmetic operation that plays crucial role in the development of new numerical procedures. The exact analytical derivatives are difficult to derive, which is why the numerical differentiation is often used instead. The automatic differentiation (AD) represents an alternative solution to the numerical differentiation as well as to the symbolic differentiation performed either manually or by a computer algebra system. With the AD technique, one can avoid the problem of expression growth that is associated with the symbolic differentiation performed by the CA system. 2.2.1 Principles of automatic differentiation If one has a computer code which allows to evaluate a function f and needs to compute the gradient Vf of f with respect to arbitrary variables, then the automatic differentiation tools, see e.g. Griewank (2000), Bartholomew-Biggs et al. (2000), Bischof et al. (2002), can be applied to generate the appropriate program code. There are two approaches for the automatic differentiation of a computer program, often recalled as the forward and the backward mode of automatic differentiation. The procedure is illustrated on a simple example of function f defined by n f = bc,withb = ^ai 2 andc = Sin(b) (1) where a1,a2,...,an are n independent variables. The forward mode accumulates the derivatives of intermediate variables with respect to the independent variables as follows { 2 xi } i = 1,2,..., n {Cos(b)ybi} i = 1,2,...,n (2) { Vbic + bVci } i = 1,2,...,n - df In contrast to the forward mode, the backward mode propagates adjoins x = -x , which are the derivatives of the final values, with respect to intermediate variables: {dai} \dxi J [dai j Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 21 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. f = — = 1 1 df c df 9f f bf (3) Vf = {aij = jjabi bJ = { 2ai b } i = 1,2,...,n The numerical efficiency of the differentiation can be measured by numerical work ratio numerical cost(f(a1,a2,a3,...,an),Vf = -^) wratio(f) =----------------------------------------------------------i (4) numerical cost(f(a1,a2,a3,...,an)) The numerical work ratio is defined as the ratio between the numerical cost of the evaluation of function f together with its gradient Vf and the numerical cost of the evaluation of function f alone. The ratio is proportional to the number of independent variables O (n) in the case of forward mode and constant in the case of backward mode. The upper bound for the ratio in the case of backward mode is wratio(f) < 5 and is usually around 1.5 if care is taken in handling the quantities that are common to the function and gradient. Although numerically superior, the backward mode requires potential storage of a large amount of intermediate data during the evaluation of the function f that can be as high as the number of numerical operations performed. Additionally, a complete reversal of the program flow is required. This is because the intermediate variables are used in reverse order when related to their computation. There exist many strategies how the AD procedure can be implemented. The most efficient are source-to-source transformation strategies that transform the source code for computing a function into the source code for computing the derivatives of the function. The AD tools based on source-to-source transformation have been developed for most of the programming languages, e.g. ADIFOR for Fortran, ADOL-C for C, MAD for Matlab and AceGen for Mathematica. 2.2.2 Automatic differentiation and FEM The tools for automatic differentiation were primary developed for the evaluation of the gradient of objective function used within the gradient-based optimization procedures or the Hessian of objective function used within the Newton-type 22 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. optimization procedures. The objective function is often defined by a large, complex program composed of many subroutines. Thus, one can apply the AD tools directly on the complete FE environment to obtain the required derivatives when the evaluation of the objective function involves FE simulation. The AD technology can also be used for the evaluation of specific quantities that appear as part of a finite element simulation. It would be difficult and computationally inefficient to apply the AD tools on large FE systems to get e.g. the global stiffness matrix of large-scale problem directly. This is especially problematic when a fully implicit Newton type procedure is used to solve nonlinear, transient and coupled problems involving various types of elements, complicated continuation or arc-length methods and adaptive procedures. However, one can still use automatic differentiation at the single element level to evaluate element specific quantities in an efficient way, such as: - strain and stress tensors; - nonlinear coordinate transformations; - consistent stiffness matrix; - residual vector; - sensitivity pseudo-load vector. A direct use of automatic differentiation tools for the development of nonlinear finite elements turns out to be complex and not straightforward. Furthermore, the numerical efficiency of the resulting codes is poor. Another solution, followed mostly in hybrid object-oriented systems, is to use problem specific solutions to evaluate local tangent matrix in an optimal way. Another solution, followed in hybrid symbolic-numeric systems, see e.g. (Korelc 2002), is to combine a general computer algebra system and the AD technology. The implementation of the AD procedure has to fulfill specific requirements in order to get element source codes that are as efficient as manually written codes. Some basic requirements are: - The AD procedure can be initiated at any time and at any point of the derivation of the formulas and as many times as required (e.g. in the example at the end the AD is used 13 times during the generation of element subroutine). The recursive use of standard AD tools on the same code, if allowed at all, leads to numerically inefficient source code. This requirement limits the use of standard AD tools. An alternative approach is implemented in (Korelc 2007b) where the source-to-source Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 23 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. transformation strategy is replaced by the method that consistently enhances the existing code rather than produces a new one. - The storage of the intermediate variables is not a limitation when the backward differentiation method is used at the single element level. The finite element formulations involve, at the single element level, a relatively small set of independent and intermediate variables. - For the reasons of efficiency, the results of all previous uses of AD have to be accounted for when AD is used several times inside the same subroutine. - The user has to be able to use all the capabilities of the symbolic system on the final and the intermediate results of the AD procedure. - The AD procedure must offer a mechanism for the descriptions of various mathematical formalisms used within the FE formulation. The mathematical formalisms that are part of the traditional FE formulation are e.g. partial derivatives —f , total derivativesf or directional derivatives. They can d(•) D(•) all be represented by the AD procedure, if possible exceptions are treated in a proper way. However, the result of AD procedure may not automatically correspond to any of the above mathematical formalisms. Let us define a "conditional derivative" with the following formalism vf = 3f(a,b(a)) 0(a) (5) _- 9(a) 9(b) =M where function f depends on a set of mutually independent variables a and a set of mutually independent intermediate variables b. The above formalism has to be taken in an algorithmic way. It represents the automatic differentiation of function f with respect to variables a . During the AD procedure, the total derivatives of intermediate variables b with respect to independent variables are set to be equal to matrix M. Some situations that typically appear in the formulation of finite elements are presented in Table 1. In case A there exists an explicit algorithmic dependency of b with respect to a , hence the derivatives can be obtained in principle automatically, without user intervention, simply by the chain rule. However, there also exists a profound mathematical relationship that enables evaluation of derivatives in a more efficient way. This is often the case when the evaluation of b involves iterative loops, inverse matrices, etc.. 24 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Case B represents the situation when variables b are independent variables and variables a implicitly depend on b. This implicit dependency has to be considered for the differentiation. In this case, automatic differentiation would not provide the correct result without the user intervention. A typical example for this situation is a differentiation that involves a transformation of coordinates. Usually the numerical integration procedures as well as interpolation functions require additional reference coordinate. An exception for automatic differentiation of type B is then introduced to properly handle differentiation involving coordinate transformations from initial X to reference coordinates ? as follows: <9(.) d(.) dX dX Ö|_[ÖX (6) In case C there exists an explicit dependency between variables b and a that has to be neglected for differentiation. The status of dependent variable b is thus temporarily changed. For the duration of the AD procedure, it is changed into an independent variable. The situation frequently appears in the formulation of mechanical problems where instead of the total variation some arbitrary variation of a given quantity has to be evaluated. The exceptions of cases A, B and C are imposed within automatic differentiation only during the execution of the particular call of the AD procedure. Case D is equal to case A with an AD exception defined globally, thus valid for every call of the AD procedure during the derivation of the problem. When in collision, then exceptions of type A, B and C overrule the D type exception. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 25 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Table 1: Automatic differentiation exceptions. Tabela 1: Izjeme pri avtomatskem odvajanju. Type Formalism Schematic AceGen input A vf öf(a,b(a)) 3(a) 9(b) 9(a) =M B Vf öf(b) ö(a(b)) 9(b) 9(a) =M C vf öf(a,b(a)) 3(a) 9(b) =0 9(.) D ö(.) ö(.) 9(b) 9(a) =M 2.3 Hybrid symbolic-numerical approach The real power of the symbolic approach for testing and applying new, unconventional ideas is provided by general-purpose CA systems. However, their use is limited to problems that lead to large systems like finite element simulations. Furthermore, the use of large commercial finite element environments to analyze a variety of problems is an everyday engineering practice. The hybrid symbolic-numerical (HSN) approach is a way to combine both. Although large FE environments often offer a possibility to incorporate user defined elements and material modes, it is time consuming to develop and test these user defined new pieces of software. Practice shows that at the research stage of the derivation of a new numerical model, different languages and different platforms are the best means for the assessment of specific performances and, of course, failures of the numerical model. The basic tests, which are performed on a single finite element or on a small patch of elements, can be done most efficiently by using general CA system. Many design flaws, such as element instabilities or poor convergence properties, can be easily identified, if the element quantities are investigated on a symbolic level. Unfortunately a standalone CA system becomes very inefficient once there is a larger number of nonlinear finite elements to process or if iterative numerical procedures have to be executed. In order to assess element performances under real conditions, 26 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. the easiest way is to run the necessary test simulations on sequential machines with good debugging capabilities and with the open source FE environment designed for research purposes, e.g. FEAP, AceFEM or Diffpack. At the end, for real industrial simulations involving complex geometries, a large commercial FE environment has to be used. In order to meet all these demands in an optimal way, an approach is needed that would offer multi-language and multi-environment generation of numerical codes. The automatically generated code is then incorporated into the FE environment that is most suitable for the specific step of the research process. The structure of the hybrid symbolic-numerical system AceGen for multi-language and multi-environment code generation introduced by (Korelc 2002) is presented in Fig. 3. Using the classical approach, re-coding of the element in different languages would be time consuming and is rarely done. With the general CA systems, re-coding comes practically free, since the code can be automatically generated for several languages and for several platforms from the same basic symbolic description. An advantage of using a general CA system is also that it provides well known and defined description language for the derivation of FE equations, generation of FE code and possibly also for the complete FE analysis, as opposed to the hybrid object-oriented systems which introduce their own domain-specific language. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 27 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Computer algebra environment • symbolic input • problem formulation . derivation of formulas • code generation • interface code • initialization . numerical integration .... -----------------------------L---------------------------- J ABAQUS MATLAB Fig. 3: Multi-language and multi-environment FE code generation. Slika 3: Več jezično Več okoljsko generiranje kode končnega elementa. When the symbolic approach is used in a standard way to describe complex engineering problems then the common experience of computer algebra users is an uncontrollable swell of expression, as pointed out before, leading to inefficient or even unusable codes. Not many attempts have been undertaken to design a general FE code generator, where this key issue of the FE code generation would be treated within the automatic procedure. The classical way of optimizing expressions in CA systems is to search for common sub-expressions after all the formulae have been derived and before the generation of the numerical code. This seems to be insufficient for the general non-linear mechanical problems. An alternative approach for automatic code generation is employed in AceGen and called Simultaneous Stochastic Simplification of numerical code, see (Korelc 1997). This approach avoids the problem of expression swell by combining the following techniques: symbolic and algebraic capabilities of general computer algebra system Mathematica, automatic differentiation technique and simultaneous optimization of expressions with automatic selection and introduction of appropriate intermediate variables. Formulae are optimized, simplified and replaced 28 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. by the auxiliary variables simultaneously with the derivation of the problem. A stochastic evaluation of the formulae is applied for determining the equivalence of algebraic expressions instead of the conventional pattern matching technique. The simultaneous approach is also appropriate for problems where intermediate expressions can be subjected to the uncontrolled swell. 2.3.1 Typical example of automatic code generation procedure To illustrate the standard AceGen procedure, a simple example is considered. A typical numerical sub-program that returns a determinant of the Jacobean matrix of nonlinear transformation from the reference to initial configuration for quadrilateral element topology is derived. The syntax of the AceGen script language is the same as the syntax of the Mathematica script language with some additional functions. The input for AceGen can be divided into six characteristic steps: - At the beginning of the session the SMSInitialize function initializes the system. - The SMSModule function defines the input and output parameters of the subroutine ”DetJ”. - The SMSReal function assigns the input parameters X$$ and k$$ and e$$ of the subroutine to the standard Mathematica symbols. Double $ character indicates that the symbol is an input or output parameter of the generated subroutine. - During the description of the problem special operators (h,H,l=) are used to perform the simultaneous optimization of expressions to create of new intermediate variables. The SMSD function performs an automatic differentiation of one or several expressions with respect to the arbitrary variable or the vector of variables by simultaneously enhancing the already derived code. - The results of the derivation are assigned to the output parameter J$$ of the subroutine by the SMSExport function. - At the end of the session the SMSWrite function writes the contents of the vector of the generated formulae to the file in a prescribed language format. The generated subroutine in C and FORTRAN language are presented in Fig. 5. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 29 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 4: Typical AceGen input. Slika 4: Tipični vhodni podatki za AceGen. Fig. 5: Typical automatically generated subroutine in FORTRAN and C language. Slika 5: Tipična avtomatsko generirana subrutina za jezika FORTRUN in C. 2.4 Abstract symbolic formulations in computational mechanics The true benefit of using symbolic tools is not about the development of a theory which is normally done manually on a sheet of paper using a pencil, or if a computer shall be used a simple word processor is adequate for such task. The advantage of the symbolic approach in computational mechanics becomes apparent only when the description of the problem, which means that the basic equations are written down, is appropriate for the symbolic description. Unfortunately, some of the traditional descriptions used in computational mechanics are not appropriate for the symbolic description. The symbolic formulation of the computational mechanics problems often differs from the classical one and thus brings up the need for rethinking and reformulating known and traditional ways. Despite that, there exist strong arguments why at the end symbolic formulations are indeed beneficial, i.e.: 30 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. - A symbolic formulation is more compressed and thus gives fewer possibilities for an error. - Algebraic operations, such as differentiation, are done automatically. - Automatically generated codes are highly efficient and portable. - The multi-language and multi-environment capabilities of symbolic systems enable generation of numerical codes for various numerical environments from the same symbolic description. - An available collection of prepared symbolic inputs for a broad range of finite elements can be easily adjusted for the user specific problem leading to the on-demand numerical code generation. - The multi-field and multi-physic problems can be easily implemented. For example, the symbolic inputs for mechanical analysis and thermal analysis can be combined into a new symbolic input that would create a finite element for fully coupled and quadratically convergent thermo-mechanical analysis. For example, the standard formulation (see e.g. Crisfield 1996, 1997, Hughes 2000, Zienkiewicz, Taylor 2000a, Zienkiewicz, Taylor 2000b) of the tangential stiffness matrix BTDB can be easily repeated using the symbolic tools. Having in mind that element tangential stiffness matrix is either the jacobian of the resulting system of discrete algebraic equations or the hessian of the variational functional, then the automatic differentiation should be sufficient for obtaining the tangent matrix. The work of implementing BTDB formulation and the efficiency of the resulting code is inferior to the approach when tangent matrix is derived by the backward AD. The latter approach requires, regardless of the complexity of the topology and the material model, a single line of symbolic input. The standard BTDB formulation would require much more input for the same result. It should be pointed out that the symbolic differentiation is one of the algebraic operations prone to severe expression growth and it can results even for relatively simple nonlinear elements in hundreds of pages of code. Thus, the use of hybrid system (e.g. AceGen) that combines the symbolic tool with the automatic differentiation technique is essential for the high abstract symbolic formulation of FE models. To increase the numerical efficiency of the generated code and to limit the physical size of the generated code, it is essential to minimize the number of calls to automatic differentiation procedure. In backward mode of automatic differentiation the expression SMSD [a,c ] + SMSD [b,c ] can result in a code that is twice as large and twice slower than the code produced by the equivalent expression SMSD [a +b,c ]. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 31 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. In this section, an abstract symbolic formulation is described, which is needed to obtain the contribution of a single element to the internal force vector ? and to the tangential stiffness matrix K. The variational functional approach and the weak form approach are the two basic possibilities open for the derivation of variational formulation of equilibrium equations and their linearizations. A weak form approach is used for the derivation of a 2D quadrilateral finite element formulation in Section 3.4. 2.5 Symbolic-numerical environment AceFEM Within the work covered by the thesis the AceFEM package (Korelc 2007a) is used for direct and sensitivity analysis. AceFEM is a general finite element environment designed to solve multi-physics and multi-field problems. It explores advantages of symbolic capabilities of Mathematica while maintaining numerical efficiency of commercial finite element environment. The main part of the package includes procedures that are not numerically intensive such as processing of the user input data, mesh generation, control of the solution procedures, graphic post-processing of the results, etc. These procedures are written in Mathematica language and executed inside Mathematica. The second part includes numerically intensive operations,a such as evaluation and assembly of the finite element quantities (tangent matrix, residual, sensitivity vectors, etc.), solution of the linear system of equations, contact search procedures, etc. The numerical module exists in two versions. The basic version called CDriver is an independent executable written in C language and is connected with Mathematica via the MathLink protocol. The alternative version called MDriver is completely written in Mathematica's symbolic language. It has the advantage of using advanced capabilities of Mathematica, such as high precision arithmetic, interval arithmetic, or even symbolic evaluation of FE quantities to analyze various properties of the numerical procedures on relatively small examples. Both environments operate from Mathematica and they also have the same data structures, functions, command language and input data (for details of the environment see (Korelc 2007a, b)). Direct and sensitivity analysis using AceFEM is further explained in Chapter 3. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 33 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 3 Direct and Sensitivity limit load analysis In Chapter 2 the symbolic approach to computational mechanics was introduced and all advances were outlined. In the present Chapter direct and sensitivity limit load analyses is explained which are used later on within the procedure for the evaluation of the most unfavorable imperfections (Chapter 4) as well within the limit load optimization procedure (Chapter 5). In Chapter 1 an overview was given over the history and the development of different methods of sensitivity analysis. Although sensitivity analysis is used in many areas of science and is by itself a major field of research in structural engineering, the scope of the present work is mostly dedicated to shape optimization. For this reason application to gradient based shape optimization will be studied. Accurate sensitivity analysis with the use of symbolic approach, which is needed for correct gradient shape optimization, will be presented. 3.1 Definition of ultimate states For the limit load structural analysis, used for the limit load shape optimization, the criteria for the limit load have to be defined. An ultimate state of a structure is generally defined with the limit point of the equilibrium path. In real, imperfect structures, this criterion proves unreliable because of possible exceeding of permissible tolerances of displacements or deformations before reaching the limit point. It is therefore necessary to additionally define the ultimate state of a structure. The ultimate state can be defined as the lowest load factor obtained by the following criteria (see Fig. 6): a) The maximum load factor on the load-deformation curve (limit load). b) The bifurcation load factor, before reaching the limit point of the load-deformation curve (does not occur in the presented case). c) The largest tolerable deformation, where this occurs during loading path before reaching a bifurcation load or limit load. 34 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 6: Definition of ultimate states (EN 1993 1-5 2004). Slika 6: Definicija računskih mejnih stanj (EN 1993 1-5 2004). When the first criterion of all criteria defined in Fig. 6 is reached on the load-displacement curve, the equilibrium point is defined as the ultimate load. Throughout the thesis ultimate load analysis is used. 3.2 Direct Analysis Nonlinear mechanical problems can be in general classified into 4 categories shown in Table 2: Table 2: Residual form of equations for mechanical problems. Tabela 2: Ravnotežne enačbe za različne probleme v mehaniki. Steady-state non-linear systems Ü>(a) = 0 Transient non-linear systems *(a,ap) = 0 Steady-state coupled non-linear systems *(a,b) = 0 $(a,b) = 0 Transient coupled non-linear systems *(a,ap,b,bp) = 0 $(a,ap,b,bp) = 0 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 35 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. In the notation used, a represents a vector of global generalized displacement parameters (displacements, rotations, enhanced mode parameters, etc.), b a vector of unknown state variables defined for each integration point (plastic deformations, hardening variables, etc.), ap a vector of generalized displacement parameters at the end of previous time step, bp a vector of state variables at the end of previous time step,? a set of equilibrium equations, and ? a set of local plastic evolution equations. For an ultimate limit load structural analysis the consideration of geometrical and material nonlinearity is necessary. According to terminology introduced by Michaleris, Tortorelli and Vidal (Michaleris, et al. 1994), the formulation of the system which has to be solved presents a nonlinear transient coupled non-linear system. A standard “arc-length” type continuation method is used for structural analysis (see e.g. Crisfield 1996, 1997). Therefore, in structural analysis the equilibrium equations are extended with load factor ? as an additional variable and constraint gc as an additional equation imposed on the increments of generalized displacements. The complete set of equations which need to be solved for each integration point and for the whole structure can be written as: * p *(a,ap,b,bp <&(a,ap,b,bp) = 0 (7) a= {a,A} ap {ap,?p} The general formulation of the fully implicit quadratically convergent direct analysis is for our case presented in Fig. 7. 36 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Structural level *(a, b(a)) = 0 ob <93> #K • — = -—— => da da ob öa ^KAa + * = 0 a := a + Aa Integration point level $(b) = 0 #KAb+$= 0 b := b + Ab Fig. 7: General formulation of direct analysis of transient coupled nonlinear problems. Slika 7: Splošna formulacija za direktno analizo tranzientnih povezanih nelinearnih problemov. The actual form of the equations in Fig. 7 is for structural finite elements (trusses, beams, shells etc.) well known and presented elsewhere. In all examples the simplest form of continuation methods, called “displacement controlled at a specific variable”(Crisfield 1996, 1997), was used. In this case the following displacement increment constraint equation gc is used: gc = um -~fum (8) where um is the actual and um the prescribed m-th scalar component of a generalized displacement vector a. Parameter 7 is used to parameterize the direct analysis. As soon as one of the criteria for the ultimate load presented in section 3.1 is reached, the analysis is stopped and the ultimate load is determined. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 37 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 3.3 Sensitivity Analysis Design sensitivity analysis is used to compute a rate of performance measure change with respect to the system design parameters variation. In structural engineering problems the system performance measure can include any quantity that may be used to characterize system behavior, such as displacements, stress, strain, energy, buckling or limit load, frequency response, weight, etc. The dependence to design parameters such as material property, sizing, shape and configuration parameters is in general implicitly defined by the laws of mechanics. Rarely, in the case of simple problems, there exists an explicit relation (see e.g. Choi, Kim 2005a). Sensitivities are obtained by derivation. The level of derivation effort required differs drastically in dependence of the nature of the problem and the approach used to evaluate sensitivity. The four general categories of mechanical problems were presented in Table 2. The derivation of sensitivity terms is significantly more complex for nonlinear systems than linear which will not be addressed here. For transient systems an additional dependency on time has to be considered in the derivations. Further complexity is gained with the choice of design parameters. In the case of sizing or material optimization (e.g. beam cross-section, shell thickness, elastic modulus etc.) the design variables appear explicitly in the variational equations, where in the case of shape optimization the design variables gain an implicit relation to the variational equations. If using FEM for structural analysis, a change in shape design variables implies a change in the finite element model. The dependency of design variables with respect to the coordinates of finite element nodes presents the main difficulty in evaluating sensitivity analysis. The limit load structural shape optimization leads to a transient coupled nonlinear system of equations where geometrical and material nonlinearities are taken into account. The difficulty of sensitivity expression derivation was an encouragement to use a symbolic-numeric approach which is thoroughly explained in Chapter 2. The sensitivity analysis based on direct differentiation method (Michaleris, et al. 1994) is used to evaluate the sensitivity of the objective function f with respect to the shape parameters ?. Due to the transient nature of the problem, the sensitivity analysis has to be evaluated at the end of each time step and integrated through the whole analysis. The corresponding equations are presented in Fig. 8. 38 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Structural level ? (a ,ap,b ,bp, * K Da D* D(f> Dej) D^ D(f> D?DapD?DbpD? ++ Dap D(f> Dbp D(f> D D b UK)"1 D$ D$ Dap D& Dap D^ + Dap Dej) + Dap Dej) Integration point level $(a ,ap,b ,bp) = 0 Db D& $ K D D0 D0 D$Da Da D1 6 2 ^ /V Ji, L 12 ^3 15 11 10 14 13 T Fig. 9: Geometry of cantilever with marked node numbers and shape parameters. Slika 9: Geometrija konzolne konstrukcije z označenimi vozlišči in parametri oblike. In[63]:= << AceFEM*; L = 40; (* Length *) H = 20; (* Width *) T=l; (* Thickness *) fy = 23.5; (* yield stress *) Em=21000; (* elastic modulus *) intord=3; (* interpolation order of FE mesh parameterization*) Ni^> = 3; (*number of geometry parameters through the whole length*) 4>= Table [ToExpression [" ToString[i] ] , {i, Ntf>}] (* Definition of shape parameters *) SMTInputData["CDriver"]; domains= {{"Domainl", "ElastPlast2DSens", {T, Em, 0.3, fy}}, {"Load", "SurfaceLoadGConst", {0, -1, T}}, {"PrescDispl", "PrescribedDispl2DY", {1}}}; SMTAddDomain[domains]; SMTAddEssentialBoundary[{"X" == OS, 0, 0}] ; SMTMesh[" Domainl", "Ql", {4, 2}, {Table [{(i - 1) L/ (Ntf>- 1) , - (H/ 2 + H/2 4>[ [i] ]) } , {i, Ntf>}] , Table[{(i- 1) L/ (Ntf> - 1) , (H/2 + H/2 4>[ [i] ]) } , {i,Ntf>}]}, "InterpolationOrder"-> intord] ; (#/. srch &) ] ; Fig. 10: AceFEM input data for the structure illustrated in Fig. 9. Slika 10: AceFEM podatki za konstrukcijo iz slike 9 . The node coordinates can be kept in symbolic form. In the example the three chosen parameters ?1, ?2 and ?3 define the upper and the lower boundary line which is defined by an arbitrary spline function. In Fig. 11 and Fig. 12 the node coordinates in 9 2 4 6 8 2 5 8 1 3 5 7 4 7 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 41 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. symbolic form are shown with respect to the use of a linear and quadratic spline function. Out[88]= Node Number X coordinate Y coordinate 1 0 10. + 1. (10. 01 - 10. 02) +10. 02 + 1. (-20.-1. (10. 01-10. 02) -20. 02 + 1. (-10. 01 + 10. 02)) 2 0 10. + 1. (10. 01-10. 02) + 10. 02 + 3 0 0.5 (-20.-1. (10. 01-10. 02) -20. 02+1. (-10. 01 + 10. 02)) 10. + 1. (10. 01-10. 02) +10. 02 4 10. 10. + 0.5 (10. 01 - 10 . 02 ) + 10 . 02 + 1. (-20. -0.5 (10. 01-10. 02) -20. 02 + 0.5 (-10. 01 + 10. 02)) 5 10. 10. + 0.5 (10. 01- 10. 02) + 10. 02 + 6 7 8 9 10. 20. 20. 20. 0.5 (-20. -0.5 (10. 01-10. 02) -20. 02+0.5 (-10. 01 + 10. 02)) 10. + 0.5 (10. 01- 10. 02) + 10. 02 10. + 1. (-20. -20. 02) + 10. 02 10. + 0.5 (-20. -20. 02) + 10. 02 10. + 10. 02 10 30. 10. + 0.5 (10. 02 - 10 . 03 ) + 10 . 03 + 1. (-20. -0.5 (10. 02-10. 03) -20. 03 + 0.5 (-10. 02 + 10. 03)) 11 30. 1O.+O.5(1O.02-1O.03)+1O.03 + 12 13 14 15 30. 40. 40. 40. 0.5 (-20. -0.5 (10. 02-10. 03) -20. 03+0.5 (-10. 02 + 10. 03)) 10. + 0.5 (10. 02- 10. 03) + 10. 03 10. + 1. (-20. -20. 03) + 10. 03 10. + 0.5 (-20. -20. 03) + 10. 03 10. + 10. 03 Fig. 11: AceFEM node coordinates for linear boundary shape approximation. Slika 11: Koordinate vozlišč z linearno interpolacijo parametrizirane mreže v AceFEM. Out[74]= Node Number X coordinate Y coordinate 1 0 10. + 1. (10. 01 - 10. 02) +10. 02 + 1. (-20.-1. (10. 01-10. 02) -20. 02 + 1. (-10. 01 + 10. 02)) 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0 0 10. 10. 10. 20. 20. 20. 30. 30. 30. 40. 40. 40. 10. + 1. (10. 01-10. 02) + 10. 02 + 0.5 (-20.-1. (10. 01-10. 02) -20. 02+1. (-10. 01 + 10. 02)) 10. + 1. (10. 01-10. 02) +10. 02 10. + 10. 02-0.25 (-2. (10. 01-10. 02) +0.5 (2. (10 . 01 - 10 . 02 ) + 1. (-1O.01+1O.03))) + 1. (-20. -20. 02 - 0.25 (-2. (-10. 01 + 10. 02) + 0.5 (2. (-10 . 01 + 10 . 02) + 1. (10 . 01 - 10. 03) ) ) + 0.25 (-2. (10. 01 -10. 02) + 0.5 (2. (10 . 01 - 10 . 02 ) + 1. (-10 . 01 + 10 . 03 ) ) ) ) 10. + 10. 02-0.25 (-2. (10. 01 -10. 02) + 0.5 (2. (10 . 01 - 10 . 02 ) + 1. (-10 . 01 + 10 . 03 ) ) ) + 0. 5 (-20. -20. 02 - 0.25 (-2. (-10. 01 + 10. 02) + 0.5 (2. (-10 . 01 + 10 . 02) + 1. (10 . 01 - 10. 03) ) ) + 0.25 (-2. (10. 01 -10. 02) + 0.5 (2. (10 . 01 - 10 . 02 ) + 1. (-10 . 01 + 10 . 03 ) ) ) ) 10. + 10. 02- 0.25 (-2. (10. 01 - 10. 02) + 0.5 (2. (10 . 01 - 10. 02) + 1. (-10. 01 + 10. 03) ) ) 10. + 1. (-20. -20. 02) + 10. 02 10. + 0.5 (-20. -20. 02) + 10. 02 10. + 10. 02 10. -0.25 (-0.25 (-2. (10. 01 - 10. 02) + 2. (10 . 02 - 10. 03) ) - 2. (10 . 02 - 10 . 03) ) + 10 . 03 + 1. (-20. +0.25 (-0.25 (-2. (10. 01 - 10 . 02 ) + 2. (10 . 02 - 10 . 03 ) ) - 2. (10 . 02 - 10 . 03 ) ) - 20 . 03 - 0.25 (-2. (-10. 02 + 10. 03) - 0.25 (-2. (-10. 01 + 10 . 02) + 2. (-10 . 02 + 10. 03) ) ) ) 10. -0.25 (-0.25 (-2. (10. 01 - 10. 02) + 2. (10 . 02 - 10. 03) ) - 2. (10 . 02 - 10 . 03) ) + 10 . 03 + 0.5 (-20. +0.25 (-0.25 (-2. (10. 01 - 10 . 02 ) + 2. (10 . 02 - 10 . 03 ) ) - 2. (10 . 02 - 10 . 03 ) ) - 20 . 03 - 0.25 (-2. (-10. 02 + 10. 03) - 0.25 (-2. (-10. 01 + 10 . 02) + 2. (-10 . 02 + 10. 03) ) ) ) 10. - 0.25 (-0.25 (-2. (10. 01- 10. 02) + 2. (10 . 02 - 10 . 03 ) ) - 2. (10 . 02 - 10 . 03 ) ) + 10 . 03 10. + 1. (-20. -20. 03) + 10. 03 10. + 0.5 (-20. -20. 03) + 10. 03 10. + 10. 03 Fig. 12: AceFEM node coordinates for quadratic boundary shape interpolation. Slika 12: Koordinate vozlišč s kvadratno interpolacijo parametrizirane mreže v AceFEM. 42 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. The design velocity field can easily be computed with the use of symbolic derivation function (D[ ]) in Mathematica, as shown in Fig. 13 and Fig. 14. Out[110]= n[110]:= Map[D[NodeCoordinates, #] S, {01, 02, 03}] -10.}, 80, 0, 0.}, 80, 0, 10.}, 80, 0, -3.75}, 80, 0, 0.}, 80, 0, 3.75}, 80, 0, 0}, 0}, 80, 0, 0}, 80, 0, 1.25}, 80, 0, 0.}, 80, 0, -1.25}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}}, 0.}, 80, 0, 0.}, 80, 0, 0.}, 80, 0, -7.5}, 80, 0, 0.}, 80, 0, 7.5}, 80, 0, -10.}, 0.}, 80, 0, 10.}, 80, 0, -7.5}, 80, 0, 0.}, 80, 0, 7.5}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}}, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 1.25}, 80, 0, 0.}, 80, 0, -1.25}, 80, 0, 0}, 80, 0, 0}, 0}, 80, 0, -3.75}, 80, 0, 0.}, 80, 0, 3.75}, 80, 0, -10.}, 80, 0, 0.}, 80, 0, 10.}}} Fig. 13: Design velocity field by symbolical derivation of FE node coordinates for the case of linear interpolation between shape parameters. Slika 13: Polje začetnih občutljivosti izračunano s simboličnim odvajanjem koordinat vozlišč mreže končnih elementov za primer linearne interpolacije mreže med parametri oblike. Out[93]= 888 0 ln[93]:= Map[D [NodeCoordinates, #] S, {01, 02, 03}] 880 0, 80 0, 880 0, 80 0, 880 0, 80 0, -10.}, 80, 0, 0.}, 80, 0, 10.}, 80, 0, -5.}, 80, 0, 0.}, 80, 0, 5.}, 80, 0, 0}, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}}, 0.}, 80, 0, 0.}, 80, 0, 0.}, 80, 0, -5.}, 80, 0, 0.}, 80, 0, 5.}, 80, 0, -10.}, 0.}, 80, 0, 10.}, 80, 0, -5.}, 80, 0, 0.}, 80, 0, 5.}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}}, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 80, 0, 0}, 0}, 80, 0, -5.}, 80, 0, 0.}, 80, 0, 5.}, 80, 0, -10.}, 80, 0, 0.}, 80, 0, 10.}}} Fig. 14: Design velocity field by symbolical derivation of FE node coordinates for the case of quadratic interpolation between shape parameters. Slika 14: Polje začetnih občutljivosti izračunano s simboličnim odvajanjem koordinat vozlišč mreže končnih elementov za primer kvadratne interpolacije mreže med parametri oblike. The design velocity field can be graphically represented as a scalar function. The x-coordinates do not depend with respect to design variables in the present example as can be seen in Fig. 11 and Fig. 12. The y-coordinate dependence with respect to design variables is plotted in Fig. 15 and Fig. 16 for the case of linear and quadratic spline interpolation between shape parameters respectively. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 43 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. dY dY dY AceFEM AceFEM Fig. 15: Graphical representation of the y component of the design velocity field for the case of linear interpolation between shape parameters. Slika 15: Grafični prikaz y komponente polja začetnih občutljivosti v primeru linearne interpolacije mreže med parametri oblike. 44 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. dY dY djl Max . 10. Mi n . -10. AceFEM dY Ö03 Max . 10. Mi n . -10. AceFEM Fig. 16: Graphical representation of the y component of the design velocity field for the case of quadratic interpolation between shape parameters. Slika 16: Grafični prikaz y komponente polja začetnih občutljivosti v primeru kvadratne interpolacije mreže med parametri oblike. 3.3.2 Exact sensitivity analysis While the design velocity field can be defined and evaluated symbolically, the numerical analysis done by computer algebra systems cannot keep up with the runtime efficiency of programming languages such as FORTRAN and C. The key idea of the used approach is to use a dual symbolic-numeric finite element environment. Such environment (AceFEM) was introduced in Section 2.5. The whole procedure of Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 45 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. the evaluation of analytical sensitivities is presented in Fig. 17 and can be applied on problems with arbitrary complexity. The sensitivity analysis is done on the basis of automatically derived finite element code explained in Chapter 2. The general expressions for sensitivity analysis are given in Fig. 8. Numeric values for design parameters Design parameters left in symbolic form Numeric mesh generation by CDRIVER Symbolic mesh generation by MDRIVER Direct analysis & Sensitivity analysis Design velocity fielcP d NodeCoordinates d(j) Design sensitivity da Fig. 17: Flowchart of shape sensitivity analysis by dual symbolic-numeric FE environment. Slika 17: Potek občutljivostne analize s pomočjo simbolno numeričnega MKE okolja. 3.4 Symbolical formulation of general finite strain plasticity For the representation of automatic derivation of internal force vector ? and the tangential stiffness matrix K, the 4-node quadrilateral elastic-plastic finite element is derived next. Let a be a vector of generalized displacements parameters of the element, b a vector of unknowns at Gauss point level and bp a vector of history values at Gauss point level from the previous time step. The elasto-plastic problem is defined by a hyperelastic strain energy density function W , a yield condition f and a set of algebraic constraints to be fulfilled at Gauss point level ?(a,b,bp) that have to be solved for unknowns b when the material point is in plastic state. In general, vector 46 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. b is composed of an appropriate measure of plastic strains (or stresses in the small deformation case), the hardening variables and the consistency parameter ? where ? are composed of the corresponding set of discretized evolution equations that describe the evolution of plastic strains and hardening variables and the consistency condition f = 0 . The yield condition is evaluated for the trial state by freezing the state variables as follows ftr = f(a,bp ) (9) The general algorithm for the abstract symbolic description of elasto-plastic problems is presented in Fig. 18. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 47 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. summation over GA US S points x:=X(L) + u(a,0 use AD exception of type B for coordinate transformation F := Ttrial:=T(a,bp) f(Ttrial) < 0 {b := bp local NEWTON loop b:= bp repeat 6>$(a,b,bp) dx dX 9X 9X 9? f(Ttrial) > 0 ¦ A := ob Ab:=-A -1*(a,b,bp) b := b + Ab until |Ab| < TOL b := b define AD exception of type D for b <9(.) db_ ^^(a,b,b p) 9a 9a (10) use AD exception of type C 9(.) K := 6>* 3a endloop Fig. 18: Algorithm for the abstract symbolic description of elasto-plastic problems. Slika 18: Algoritem za abstrakten simbolni zapis elasto-plastičnega problema. Here b denotes a vector of the local unknowns at Gauss point level within the iterative loop. A is a matrix that follows from the linearization of the nonlinear equation set $. The "basic equation of the symbolic plasticity" is written as: *=; dW a 9a 9b W) dn (11) =0 48 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. An efficient and accurate numerical solution of the corresponding coupled non-linear system of algebraic equations requires a quadratically convergent numerical procedure. For this the linearization of (11) is needed, which leads to the tangent stiffness matrix. This matrix can be derived for a finite element by directly applying the automatic differentiation procedure leading to K <9* da (12) Tangent stiffness matrix derived in this way is already "consistent" with the algorithm used for plasticity. Hence no additional procedures to derive a consistent tangent modulus are required. The parts necessary for the abstract symbolic description are briefly summarized in Fig. 19. Fe = FF"l1 Ce = Fe Fe J2 = Det(Ce) W = 2 (tr(Ce)-3-ln(J 2 )) + 4 (J 2 - 1-ln(J2)) dW T t = 2Fe ——Fe dC tr(x) s = T------ 1 3 a = yj2/ 3A f = Vs • s - V2/ 3(Y0 + Ha) $ = Fe-exp (H^-App l ) )FFppl-1 = 0 if = 0 b= \^pl , \l\ w = {vpp lr1,%} Fig. 19: Summary of the finite strain plasticity equations. Slika 19: Povzetek enačb plastičnosti. where Ce is right Cauchy-Green tensor and Fe is the deformation gradient. µ and ? are the first and the second Lame's material constants and ?pl is the plastic Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 49 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. multiplier. b and bp are the vectors of state variables at the current step and at the end of previous time step, respectively. 3.5 Finite element models In the work covered by the thesis 5 types of finite elements were used: - 2D and 3D Point Load finite element - 2D and 3D Line Load finite element - 3D Truss finite element - Quadrilateral, 4-node, elastic-plastic finite element - 6-parameter, elastic-plastic, shell finite elements All finite elements were derived and coded with the help of AceGen (Korelc 2007b). The equations used in the general algorithm for the abstract symbolic description shown in Fig. 18 for the used elements are summarized in Table 3. The load finite elements are used only to apply load on the numerical model. 50 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Table 3: Strain energies for the use in the Automatic differentiation exceptions. Tabela 3: Izjeme pri avtomatskem odvajanju. Element Equilibrium Equations on FE level Point Load Z, w, Fz A Y, v, F P – point load X, u, F *= P ?u ?a Line Load Z, w, Fz A Y, v, F q - continous load X, u, F Truss Z, w, Fz A Y, v, F X, u, F /12 rcW W = -AL0Eme2 dL ; *= I ^— 2 J c/a dO 2d quadrilateral Y, v, Fy -------> X, u, F x See Section 3.4 3d shell Z, w, Fz Y, v, F X, u, F See references (Wisniewski, Turska 2000, 2001). 1 2 2 2 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 51 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 3.6 Example An example of sensitivity analysis of a single-storey steel building is presented. The finite element model of the structure is shown in Fig. 20. The model consists of the following parts: - The main structure consists of four portal frames modeled by the four node shell elements based on finite rotations, 6 parameter shell theory combined with ANS and two enhanced modes for improved performance (Wisniewski, Turska 2000, 2001) - The purlins and braced system are modeled by large displacement truss elements. - Special “load” elements were generated to apply wind and snow loads. The analytical shape sensitivity pseudo-load vector is derived for all elements by direct differentiation method and with the use of symbolic code generation explained in Chapter 2. 52 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 20: Finite element model for the one story building. Slika 20: Model konstrukcije enoetažne hale. In the example the angle of the roof (ß) is used for the design shape parameter. Fig. 21 presents the typical symbolic form of the nodal coordinate generated by the MDriver. Differentiation of the symbolically parameterized mesh with respect to ß results in a design velocity field that is used within the sensitivity analysis (Michaleris, et al. 1994). The sensitivity of the vertical displacement is presented in Fig. 23. The results of analytical sensitivity analysis are then compared with the results obtained by the finite difference method in Fig. 22. The finite differences are computed considering a relative perturbation size of 9.5·10-8, for which an optimal perturbation size study has to be done. With the evaluated optimal perturbation size finite differences coincide with the analytical method with an average relative error of 4·10-5. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 53 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 4950. + 4600. Tan[/3] + 0.3 I 50. - 4600. Tan[/3] 2.76x10 9950. + 4600. Tan[/3] 0.666667 1-50. +4600. Tan[ß] + 2.76 x10 6 9950. + 4600. Tan[/3] ^ 0.8 1-400. +0.3 (-50.-4 600. Tan[ß] + 0.666667 (50. +4600. Tan[ß] ) ) ( 2.76xlOb f 2.76xlOb }}} 0.3 I 50. - 4 600. Tan[/3]------------------------------------+ 0.666667 1-50. +4600. Tan[ß] + --------------------------------- I I I______________________9950. + 4600. Tan [/3]_____________________________________9950. +4600. Tan[/3] III Fig. 21: Example of a nodal point coordinate in symbolic form. Slika 21: Primer koordinat vozlišča mreže končnih elementov v simbolni obliki. ¦ Analitic (elastic) •X"» Finite differeces (elastic) ¦Analitic (plastic) •K-" Finite differeces (plastic) 10 20 30 40 50 Roof angle ß [°] Fig. 22: Comparison between analytic and FD method. Slika 22: Primerjava analitične občutljivosti in občutljivosti po metodi končnih diferenc. + 0 54 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 23: Graphical sensitivity representation of the vertical displacement with respect to the roof angle. Slika 23: Grafični način prikaza občutljivosti vertikalnega pomika glede na naklon strehe. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 55 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 4 Initial imperfections 4.1 Introduction Limit load design optimization requires the use of imperfect geometry of a structure in order to evaluate the correct limit load. While the imperfections are not known in advance, a method has to be used which is capable of proper involvement of imperfection effects. It is now well known that geometrical, structural, material and load imperfections play a crucial role in the load carrying behavior, especially of thin walled structures. The determination of the “definitely worst” imperfection can be formulated as a nonlinear optimization problem solved by one of the well-known nonlinear optimization methods as explained in Section 1.2.3. It is doubtful that an approach that requires a large number of fully nonlinear analyses could be used in everyday engineering design in foreseeable future. In the present work a computationally less expensive optimization method is developed that would still retain the generality of the optimization based “definitely worst” imperfection approach. Geometrical, structural and material imperfections are considered by means of equivalent geometrical imperfections. The basic idea of the approach is to replace the nonlinear optimization problem with an iterative procedure that would involve only linear optimization problems. Within the iteration the objective function for the minimum ultimate load is constructed by the means of a fully nonlinear direct and first order sensitivity analysis. Constraints on the shape and the amplitude of the imperfections have to be taken into account. When carefully constructed, they remain linear, thus enabling the use of efficient and readily available linear programming algorithms for the solution of the corresponding optimization problem. In the case where only the amplitude of the imperfections is constrained, numerical studies show that the procedure tends to lead to significantly lower ultimate loads than experimentally observed, as the result is associated with imperfection shapes that are not necessarily technologically feasible. Therefore it is essential to take into account realistic technological constraints on the shape and the amplitude of the imperfections. In the present approach the imperfections are represented by a linear 56 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. combination of base shapes with the base constructed from the sufficient number of buckling modes augmented by the eigenvectors of the structure subjected to “technological” boundary conditions and characteristic deformation modes. The construction of the shape base requires the solution of several generalized eigenvector problems and is done only once. The computational cost of each iteration is equivalent to one full nonlinear analysis up to the ultimate load accompanied by the nonlinear sensitivity analysis with respect to all base shapes and the solution of the linear optimization problem. Thus the total computational cost remains within the range that is acceptable for design procedures. 4.2 Optimization method for the determination of the most unfavorable initial imperfection 4.2.1 Representation of imperfections The applied initial imperfection shape with specified amplitudes has to represent a change in the geometry of a structure in the most unfavorable way so that the ultimate load of the imperfect structure is the smallest possible. The imperfections are represented as a linear combination of the chosen base shapes within amplitude e0 prescribed by the principle of equivalent geometrical imperfections. Equivalent geometrical imperfections include geometrical and structural imperfections. Geometrical imperfections represent a general deviation from the perfect geometry. Geometrical imperfections can be augmented to include structural imperfections that are not included into the finite element model directly. Structural imperfections arise from the manufacturing method, for example residual stresses produced by welding. The geometry of an imperfect structure X is defined by: N X = Xp + ^ctjTj , (13) where Xp is the initial perfect geometry, aj are the unknown shape parameters and Yj are the base shapes. The unknown shape parameters aj are obtained as a solution of the optimization problem. The base shapes can be chosen arbitrary, but they have to be linearly independent in order to have a well defined minimum of the corresponding optimization problem. The overall numerical efficiency of the procedure strongly depends on the number of base shapes (N). The obvious choice, well explored by other authors (see e.g.Song, et al. 2004), are buckling modes (TA) of the structure obtained by initial buckling analysis. Alternative and cheaper to Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 57 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. evaluate are the eigenvectors (?B ) of the initial elastic tangent matrix K0 . The kinematic boundary conditions for imperfections can be different than kinematic boundary conditions of the structure. This can be observed in rigid support connections where the member is usually considered clamped and it is not possible to describe the support imperfections with the eigenvectors of the original structure. For this purpose the base can be extended by eigenvectors (?C ) of the elastic tangent matrix K0 of the same structure but with different kinematic boundary conditions. In this way the technological imperfections can be added. Some authors have observed (Schneider 2006, Schneider, Brede 2005, Schneider, et al. 2005) that sometimes the most unfavorable imperfection resembles deformation shapes rather than buckling modes. In order to reduce the total number of the necessary considered base shapes, an additional set of deformation shapes (?D ) of the structure in elastic and plastic range can be added. And finally, the set of shapes which are empirically known to represent the worst imperfections for certain type of structures (?E ) can be added. The total base ? is then in general composed of: ?=?A??B??C??D??E (14) The optimized imperfection shape (most unfavorable initial imperfection shape) depends on the number of shapes included in the shape base and the density of the finite element mesh. For reasonable results it is necessary to increase them proportionally. With the increase of the considered shapes, the result converges to a final shape. For practical reasons it is necessary to include at least that much different shapes to allow including all local and global collapse mechanisms. The shape of the most unfavorable initial imperfection changes with different loading patterns, supporting conditions, changes in geometry or the amplitude of initial imperfections. In this sense, the shape of the most unfavorable initial imperfection in means of ultimate load of a structure has to be evaluated for every individual structure separately and can not be generalized. 4.2.2 Description of the algorithm In the presented approach a fully geometrically and materially nonlinear analysis is used. When dealing with thin-walled structures with moderate thickness, it is necessary to take geometrical and material nonlinearity into account. Since the algorithm starts from the beginning with the imperfect structure, bifurcation points usually do not occur prior reaching the limit point in a load-deformation curve. By limiting the analysis to limit points one avoids tedious procedures involved in proper determination, classification and sensitivity analysis of bifurcation points. Only stable 58 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. equilibrium states have to be considered and no hypothetic states need to be taken in account. Within this method the most unfavorable initial imperfection shape is sought, defined by the shape base ? and the shape parameters ? at which the ultimate load will be the lowest. Unknown shape parameters ? are evaluated iteratively by an optimization process. The iterative procedure for the k-th step can be written as: Xk — Xk-1 + AXk AXk N i=1 aik = a^1 + Aaik Xk = N i = 1 (15) where Xk is the imperfect geometry, Aaik the increment of the imperfection parameters, AXk the increment of the imperfection and Xk the total imperfection. The increment of the imperfection parameters in the k-th iteration Aaik is obtained as a solution of the corresponding optimization problem described in Section 4.2.2. The flowchart of the method is illustrated in Fig. 24. The algorithm starts with the first base shape T1, normalized by the amplitude e0 , as the initial guess X0 for the geometry of the imperfect structure: ai 0 = 0; Aai 0 = ¦ e 0 i = 1 maxTi 0 i^1 (16) X0 = Xp + A«10 r\ and then improves the solution by solving a sequence of optimization problems until the convergence condition | Aaik | < tolerance is reached. Within each step of the iterative procedure a fully nonlinear direct and sensitivity analysis of the structure with imperfect geometry Xk is performed followed by the formulation and solution of the optimization problem. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 59 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. INITIAL GEOMETRY Xp INITIALIZATION Shape BASE r=rAu rBu rCu rDu rE a? = 0 ; Aaf —— i = 1 max Ti ; k 0 i ^ 1 ITERATIVE PROCEDURE Shape parametrization with parameters A -------------- analysis dA Direct *lk —> \ analysis Minimization problem min Xlk ; Xlk = Xlk N E i=1 N dX <9Aa,k • Aa, -> Aak C(Xk,e0) < 0; Xk = ^ri(aiM + Aai i = 1 If Aaik > tol k = k + 1 If Aaik < tol MOST UNFAVORABLE GEOMETRY N X = X p+E La j=l ^ m=0 N rj = Xp + XXrj = Xp + X Fig. 24: Flowchart of the method for the determination of the most unfavorable initial imperfection. Slika 24: Potek metode določitve najbolj neugodne začetne nepopolnosti. 60 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. The presented approach gives the advantage of use of arbitrary state of the art optimization algorithms, as the optimization part is completely separated from the direct and sensitivity analysis. An alternative approach would be performing a fully coupled nonlinear optimization where the most unfavorable imperfection shape parameters would be determined simultaneously within the direct and sensitivity analysis. Such an evaluation is to authors experience not feasible for larger structural systems at this time. 4.2.3 Formulation of the optimization problem The ultimate load factor of the imperfect structure in k-th iteration ?lk represents the minimizing function, where the maximal amplitude of the total imperfection Xk has to be equal to or smaller than the amplitude of the prescribed equivalent geometrical imperfections e0 : min \lk C(Xk,e0) < 0 (17) where C(Xk,e0)is a constraint function. The decoupling of the direct analysis and optimization is achieved by expansion of the ultimate state load factor ?lk to a Taylor series around the ultimate state load factor of the current imperfect geometry. The ultimate load factor is then written as: \ k ¦ \ k M — M N Ak =0 + E d\k dAak Ak =0 Aq (18) where ?lk Ak =0 is the evaluated ultimate load factor of the structure and dXlk dAa Ak =0 the sensitivity of the ultimate load factor with respect to optimization parameters Aaik in the current step. The minimizing function Alk in this instance is a linear function. The constraining function, on the other hand, can be a highly nonlinear function or a simple set of linear constraints, depending on way it is defined. The employment of constraints (17) arises from the demand of the technical standards (see e.g. EN 1090/2 2007) which specify requirements for execution of structures or manufactured components of structures. The location of the point of Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 61 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. maximal amplitude is unpredictable, which makes it very difficult to choose the appropriate restraining function. For small order problems (i.e. up to 5000 nodes) where nonlinear optimization algorithms can be used, a simple norm of the total imperfection (Deml, Wunderlich 1997) gives satisfactory results. The norm that proved to be reliable for small order problems is the L2p vector norm with the exponent 2p. With the increase of the constant p the norm limits to the maximal value of components of the total imperfection vector Xk: lim X I = max X 2 p^oJ k2p k ; N , (19) Xk 2p 2p The minimization problem (17) with inequality constraints can be solved by an advanced penalty method or an extended Lagrange multiplier type method. For large problems the necessity of a high exponent p of the L2p norm in order to achieve the necessary accuracy makes the huge constraint function highly nonlinear and the minimization problem difficult to solve. In this case it is necessary to define a set of linear constraints for the maximal amplitude of the total imperfection vector: \Xkm\ =*----P ---------------- L= 2 m ------------! V y b) i ---------1 »» 29 ^S -_¦ ^^ ^^_^^^____^_ _ 0.196 Eg^^s»^™™ ' """""^^tes^^ -5f ^*""* 41 I5 0.194 -5[ ^^* Fig. 27: Convergence to the most unfavorable shape by increasing the number of considered shapes in the shape base (scale factor fs = 10). Slika 27: Konvergenca najbolj neugodne oblike z večanjem števila upoštevanih oblik (faktor povečave fs = 10) The shape and the amplitude of the equivalent geometric imperfections should be in general chosen in such a way that it has the same effect on the load bearing capacity of the structure as all relevant imperfections together. The shapes and amplitudes for geometric imperfections can be chosen in accordance with the manufacture tolerances (e.g. EN 1090/2 2007) although it is possible that taking the amplitude of the considered initial imperfections equal to the manufacture tolerances can lead to a too low characteristic resistance. This can be even more pronounced where several different imperfections interact (Johansson, et al. 2007). There is little information about equivalent geometrical imperfections for a general structure found in existing technical standards. In the present example a comparison has been made for the ultimate load factor of the structure, considering different manually defined combinations of base shapes and the calculated most unfavorable initial shape. In Table 4 there are ultimate load factors (Xu) shown for the cantilever structure with different initial imperfections considered according to the specified combination methods. The contribution of the chosen shapes from the shape base to the final result is represented by values of parameters ai. The smallest calculated ultimate load factor resulting from the use of various combinations of imperfection shapes is 0.219 for the combination method f = ri+rj where the most unfavorable combination of two base shapes is Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 67 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. considered. Other combination methods produce higher ultimate load factors. The ultimate load factor of the structure obtained by the proposed method, taking the calculated most unfavorable initial imperfect shape into account, is 0.194, which results in a 11% smaller ultimate load with respect to other combination methods. If the number of considered base shapes was increased, the corresponding most unfavorable imperfection shape would lead to an even lower ultimate load factor. To prevent the curvatures of the imperfection shape to exceed common values, the number of considered base shapes 41 was chosen in this example. Table 4: Ultimate load factors for various initial imperfection shapes for the cantilever beam example. Tabela 4: Mejni obtežni faktorji za različne kombinacije obtežb za primer konzole. Perfect geometry ^u= 0.366 Initial geometry according to recommended combinations of shapes: Shape combination method Combination at minimal Xu K r = ri+rj i = 41 i = 1, j = 41 0.268 0.219 N t = Ti+0.7^j j=1 j*i N f = ri+rj+o.7Xrk k=l k*i,j i = 41 i = 1, j = 41 0.342 0.343 r = max T Most unfavorable initial imperfection by proposed method: -2997 I\ + 106.6 T2 + 898 T3 + 33.5 T4 + 62.4 T5 - 4.4 T6 - 17.4 T7 + 3.3 T8 + 4.3 r9 + 1.23 r10 - 3.47 T11 + 0.41 T12 - 5.25 r13 + 1.30 T14 + 2.07 T15 - 0.04 T16 + 2.06 r17 - 1.28 r18 + 0.28 r19 + 0.04 r20 - 2310 r21 - 3722 r22 + 1858 r23 - 762 r24 + 0.71 r25 + 105 r26 + 10.0 r27 - 44.5 r28 - 18.9 r29 + 4.2 r30 + 11.8 r31 - 2.4 r32 - 7.7 r33 - 5.9 r34 - 3.5 r35 - 1.8 r36 + 1.3 r37 + 4.6 r38 - 0.3 r39 - 0.4 r40 + 0 J- / X,„ = 0.194 e o 68 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 4.3.2 Thin-walled T beam Structures composed of thin–walled components in general prove a high degree of imperfection sensitivity. The thin-walled girders in this section were modeled by elasto-plastic four node shell elements based on finite rotations, 6 parameter shell theory combined with assumed natural strain formulation and two enhanced strain modes for improved performance (Wisniewski, Turska 2000, 2001). The example refers to the ultimate load calculation of a simply supported thin-walled beam with a T cross-section, loaded with a concentrate force at the mid-length. The geometrical details and loads are presented in Fig. 28. The example was taken from (Lanzo, Garcea 1996) Cross section 38 h-----------------1 FEM model Constraints u[0,0,-32.5]=0 v[0,0, 32.5]=v[450,0, 32.5] w[0,y, z ]=w[450,y, z ] Load ^~ z,w X x,u 450 Fig. 28: Geometry, supporting and loading conditions for the thin-walled T beam example. Slika 28: Geometrija, podpore in obtežba tankostenskega T nosilca. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 69 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. B H 2 x B/200 2 x H/200 Fig. 29: Maximal allowable deviation from the perfect geometry used in the evaluation. Slika 29: Največje dopustno odstopanje od popolne oblike pri izračunu. Different constraint conditions (20) were used for the flange and the web of the girder. The maximal amplitude of the equivalent geometrical imperfection for the web was taken as H/200, where H is the height of the web and the amplitude for the flange was taken as B/200, where B is the width of the flange. Within the optimization problem (17) it was therefore necessary to define 3150 constraint equations for the maximal initial imperfection amplitude perpendicular to the web and 2025 constraint equations for the maximal imperfection amplitude perpendicular to the flange. In this case, only the y-direction web components and the z-direction flange components of the total imperfection vector Xk are constrained. First the structure is analyzed considering the shape base ? consisting of buckling modes. In Fig. 30 the calculated limit load of the T-beam with increasing number of base shapes is shown. The results show clear convergence of the calculated limit load. 70 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 3200 ?l 3000 2800 2600 2400 2200 2000 0 20 40 60 80 100 Number of considered shapes Fig. 30: The calculated ultimate load of the T-beam considering the evaluated most unfavorable initial imperfection varying the number of base shapes. Slika 30: Mejna obtežba T nosilca pri upoštevanju izračunane najbolj neugodne nepopolnosti pri različnem številu oblik v bazi. Furthermore the structure is analyzed by two different shape bases in order to assess the influence of the chosen shape base on the numerical efficiency of the proposed procedures. The first shape base (sbA) consists of 50 buckling modes (?A) and two deformation shapes ( ?D ) and the second (sbB) of 50 eigenvectors of the elastic tangent matrix K 0 (?B) and two deformation shapes ( ?D ). In Fig. 31 the ultimate load-deformation curves for the two cases are plotted. The difference in results between the two cases is small and decreases with increasing the number of considered shapes. In the present example the shape base composed of eigenvectors of the elastic tangent matrix (sbB) turned out to be more appropriate than the shape base composed of buckling modes (sbA), since a lower ultimate load was computed with the same number of considered base shapes. The shape base (sbB) is used for further analysis. The most unfavorable initial imperfection evaluated by the presented approach for the shape base (sbB) and the corresponding deformed state at collapse are presented in Fig. 32. 2 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 71 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 3500 3000 2500 2000 1500 1000 500 0 / / / / r = rBU td ----perfect geometry 9 J '¦ I ¦J 0 0.2 0.4 0.6 0.8 1 Displacement [cm] Fig. 31: The load displacement curves considering various imperfection shape bases for the T beam example. Slika 31: Graf sila-pomik za T nosilec z upoštevanjem različnih baz oblik. a) 8p*~ Scale factor fs=20 b) ÄSfe Scale factor fs=1 Fig. 32: The most unfavorable initial imperfection shape of the T cross-section thin-walled girder (a) and the corresponding deformed shape at collapse (b). Slika 32: Najbolj neugodna začetna nepopolnost za primer T nosilca (a) in pripadajoče stanje pri limitni obtežbi(b). Fig. 33 presents the convergence of the global iterative procedure for the considered shape base (sbB). The result of the first iteration, where the first eigenvector is taken for the initial imperfection, gives a good approximate to the final result. Only 9 72 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. iterations were necessary to achieve convergence within tolerances and to determine the most unfavorable initial imperfection. 3100 3000 2900 2800 2700 2600 2500 2400 2300 2200 23456789 Iteration Fig. 33: Convergence of the global iteration process of finding the most unfavorable imperfection shape for the example with 52 base shapes. Slika 33: Konvergenca globalnega iteracijskega procesa iskanja najbolj neugodne oblike z upoštevanjem 52 baznih oblik.. In Fig. 34 equilibrium paths calculated for various values of amplitudes of equivalent geometrical imperfections are plotted. The e0w is the maximal amplitude of the equivalent geometrical imperfections perpendicular to the surface of the web with height H and the e0f is the maximal amplitude of the equivalent geometrical imperfections perpendicular to the surface the flange with width B. The calculated ultimate load depends to a large extent on the amplitude of the used initial imperfection. The choice of the amplitude is therefore very important and a crucial part for determining the most unfavorable initial shape. For the purpose of comparison the equilibrium path calculated on a basis of the Koiters asymptotical approach to nonlinear instability for the same example was taken from (Lanzo, Garcea 1996). Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 73 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 3500 3000 2500 2000 1500 1000 500 0 e0w=H/200, e0f=B/200 ___ • Koiters appr. (Lanzo,1996) ¦ perfect initial geometry • u .................................»/... ^ e0w=H/100, e0f=B/100 tXSt e0w=H/300, e0f=B/300 e0w=H/400, e0f=B/400 e0w=H/500, e0f=B/500 i ^^^ ^^N*^^^ ¦ I ^^^^^^^^^^» ¦ 1 0.2 0.4 0.6 0.8 Displacement [cm] Fig. 34: Equilibrium paths for the T cross-section thin-walled girder. Slika 34: Ravnotežne poti za primer T nosilca. 4.3.3 Thin-walled I beam In the third example, the most unfavorable imperfection for a standard 8m long HEA400 structural steel I beam is computed. The beam is fully rigidly supported at the ends. The elastic modulus has been taken as 21000 kN/cm2 and the yield stress as 23.5 kN/cm2. The Poisson ratio was taken as 0.3. The structure is modeled by the same finite elements as in previous section. A vertical line load is applied along the upper flange center line. The considered shape base ? consists of 58 eigenvectors of the elastic tangent matrix K0 and 2 deformation shapes corresponding to the elastic deformed state and the plastic limit state. The optimization problem includes 3111 constraint equations. All components of the total imperfection vector are constrained with the amplitude e =H/200 , Xk 0 where the height of the cross-section H = 39cm (see Fig. 35). 0 1 74 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 35: Maximal allowable deviation from the perfect geometry used in the evaluation. Slika 35: Največje dopustno odstopanje od popolne oblike pri analizi. In Fig. 36 equilibrium paths are plotted for various combinations of imperfection shapes taken from T . Shapes and combinations of shapes considered are: a) initial imperfection in the shape of elastic deformed shape, b) initial imperfection in the shape of the plastic deformed shape, c) combination of two shapes ri + Tj, N d) combination of all shapes in form Ti+Tj+ 0.7 ^ Tk , k=1 ;k^i,j e) computed most unfavorable initial imperfection shape for shape base r. In the case (c) as in the case (d) the minimum limit load calculated was achieved for i = 34 and j = 37. All considered initial imperfection shapes were normalized by the value of the equivalent geometrical imperfection amplitude e0 for the purpose of comparison. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 75 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 1.4 1.2 0.8 0.6 0.4 0.2 01 234 Displacement [cm] Fig. 36: Ultimate load deformation curves for different initial imperfect geometries with the same amplitude for the I beam example. Slika 36: Mejne krivulje za različne oblike začetnih nepopolnosti z enako amplitudo za primer I nosilca. The corresponding imperfect geometries are plotted in Fig. 37. For case (e) two load curves are plotted belonging to the first iteration of the global iterative procedure and the final converged state. The various recommended combinations of shapes (a-d) result in significantly higher limit load when compared to the ultimate load obtained by the new approach. 1 0 76 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 37: Initial imperfect geometries used in analyses according to Fig. 36. Slika 37: Prikaz začetnih nepopolnosti. 4.3.4 Thin-walled Cylinder Among all thin-walled structures axisymetric structures (e.g. spheres and cylinders) prove to have the highest imperfection sensitivity. Several papers deal with the problem of finding the initial imperfection connected to the lowest ultimate load of cylindrical structures (Schmidt 2000, Schneider 2006, Schneider, Brede 2005, Schneider, et al. 2005, Schranz, et al. 2006, Song, et al. 2004). In the early stages of imperfection studies on cylinders, the analogy to column and plate buckling was considered appropriate and therefore the imperfection affine to the lowest eigenmode was taken as the worst initial imperfection. Recent studies dealing with the direct determination of the worst initial imperfections suggest that single dimple imperfections may be worse than eigenmode-affine patterns covering the whole structure (Wunderlich, Albertin 2000, 2002). In the present example an axially compressed cylinder is studied Fig. 38. The cylinder is fully rigidly supported at z=0 and free at z=H. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 77 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. t/R= 1/200 R=H=1000 cm E=21000 kN/cm2 fy=23.5 kN/cm2 v=0.3 1000 -500 X ¦ py= 0.01 A, • pz 1000 500 z 1000 500 1000-1000 Fig. 38: Geometry of the axially and transversely loaded cylinder. Slika 38: Geometrija osno in prečno obremenjenega cilindra. After the convergence study, the shape base was chosen that consists of 68 shapes including 59 eigenvectors (TB) of the elastic tangent matrix K0, 7 buckling modes (YA) and 2 deformation shapes (YD) corresponding to the elastic deformed state and the plastic ultimate state. Technical standards (EN 1993 1-6 2006) prescribe the maximal amplitudes for the equivalent geometrical imperfections for cylindrical structures. In the present example the maximal amplitude of the equivalent geometrical imperfections perpendicular to the cylinder wall was taken as e0 = 1.22t as for class A fabrication tolerance quality. In order for the constraints to remain linear and to preserve the possibility of using linear optimization methods the projection of the imperfection vector X" = {X",X",X"} in the n -th node in the radial direction was chosen to be bounded. The resulting optimization problem includes 2888 linear constrained equations of the form |Xrarra|<) ~d4> ~ dcf> Optimization environment Mathematical programming Convergence criteria reached Optimum result Fig. 44: Optimization loop using symbolic-numeric environment. Slika 44: Optimizacijska zanka s pomočjo simbolno-numeričnega okolja. Within the optimization process (see Fig. 44) the geometry of the structure is updated in two different loops. In the inner loop the geometry X(?) is changed due to the change of the design parameters done by the optimization algorithm. In the Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 91 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. outer loop the geometry is changed due to changed initial imperfections X. The method for evaluating the most unfavorable initial imperfections is presented in Fig. 24. The geometry can be written as: X = XP(0) + X (24) where XP() —— =^i = —-p^ (25) 8(f) 8(f) 8(f) The design velocity field does not depend on the total imperfection vector nor does it change by varying the design parameters. Therefore, it has to be evaluated only once. The most unfavorable imperfection is evaluated whenever the geometry changes to a certain extent. Together with the updated geometry the most unfavorable imperfection of the current geometry changes simultaneously. While the evaluation of the most unfavorable imperfections is a computationally demanding task, this is not done in every step of the optimization process. The initial imperfection is changed whenever the most unfavorable imperfection would cause a change in the limit load higher than a prescribed value AL or would evolve in a change of the limit load state. Conveniently, the basic shape of the most unfavorable imperfections generally does not change if minor updates of design variables are done by the optimization algorithm. Further on, the convergence of the optimization algorithm is better preserved when using the same initial imperfections and to update the initial imperfections sequentially on sets of optimal solution procedures. The practical application of the developed algorithms is shown in numerical examples in the next chapter. 92 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 5.4 Numerical examples 5.4.1 2D cantilever shape optimization To illustrate the limit load optimization procedure, a simple case of a 2D cantilever is studied first. The mathematical model consists of elastic-plastic finite strain, 2D, quadrilateral elements. An ideal elasto-plastic material is used. The geometry (see Fig. 45a) is parametrized with parameters

) has to be exactly 1. The goal is to minimize the volume of the structure: minf; f = V(4>) (27) <&kP < 0 The constraints include a limit load factor equality constraint XL = 1 which forces the calculated limit load factor XL((f)) to take the prescribed value 1 at which the calculated limit load q matches the prescribed limit load of the structure qu , and a set of inequality constraints <3kP < 0 which prescribe the minimum feasible values of parameters (f). The constrained problem is solved using an interior point method considering the merit function fB : fB = w1V + w2Log(XL - 1) + w3^2$kP (28) k where wi are the weights, V is the volume and <&P is barier function. The problem is solved with a standard quasi-newton algorithm using Mathematica (Wolfram 2008). Fig. 45 shows the initial shape (a) and the optimal shape (b). Mises stresses are plotted. The structure with the optimal shape shows a smooth distribution of yield stress over the whole length of the structure. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 93 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. load q = ? ?qu h0 MM j -4JM A > lllllllllill 1 L * va b h (/2,(x)= ——x , (29) PL Wpl(x) > P, p2)anaiyt. = 0) can be observed as the finite element mesh with the quad finite elements on the free end can not form an analytically sharp edge and there is no consideration of shear deformation in the analytical approach. 5.4.2 3D H cross-section thin-walled cantilever structure In the present example a thin-walled cantilever structure is studied. Structures composed of thin-walled components in general prove a high degree of imperfection sensitivity. Therefore the use of imperfections in an analysis is mandatory for correct optimization results and the flow of the optimization process itself, as possible bifurcation points in the analysis are avoided and a realistic lowest limit load can be calculated. The specific shape of the structure was chosen for representation purposes, as a variety of collapse mechanisms can be observed at the limit load. In conventional shape optimization dealing with bearing capacity of structures the shape is usually optimized for the state of the structure at the time when the first material point exceeds the elastic resistance or the first member buckles. No post-buckling or post-critical behavior is taken into account. In the present approach the collapse mechanism and the phenomena appearing at that time dictate the optimal shape of the structure. In this way, the shape of the structure is sought, which gives the maximal bearing resistance at the limit state, which is usually presented as a collapse of the structure. Varying the thickness of shell components by excluding it from the design variables specific collapse mechanisms can be enforced. The limit state of the structure is characterized by pure plastic limit state when the sheet components are thick enough. Lowering the thickness of sheet components results in a plastic buckling limit state and further on in an elastic buckling limit state. It has to be mentioned that within conventional shape optimization only one collapse mechanism is considered at the 96 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. same time and therefore the optimal shape is not computed considering more possible collapse mechanisms. The geometry of the studied cantilever structure is shown in Fig. 48. The thin-walled cantilever is modeled by elastic-plastic four-node shell elements based on finite rotations, 6 parameter shell theory combined with assumed natural strain formulation and two enhanced strain modes for improved performance (Wisniewski, Turska 2000, 2001). An ideal elasto-plastic material model has been used. case A case B E= 210000 MPa fy= 235 MPa P= 120kN A- P 20 0 -20 Z Y Fig. 48: Initial geometry of H cross-section cantilever. Slika 48: Začetna geometrija H konzolnega nosilca. Two fixed sets of wall thicknesses are considered: - Case A: thick sheets, where plastic behavior dominates. - Case B: thin sheets, where buckling behavior dominates. The shape of the structure is parameterized with parameters ? as shown in Fig. 49. The geometry is symmetrical with respect to the XZ plane. The boundaries of the vertical sheets are varied in both Y and Z axes. The vertical sheets have to remain vertical and the horizontal sheet horizontal. The minimum dimension of both, E Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 97 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. vertical and horizontal sheet is 10% of the initial dimension. Second order splines were chosen for interpolation of the boundary shape between the shape parameters. '9 20 Z y9™ -20 Fig. 49: H cross-section cantilever shape parameters. Slika 49: Parametri oblike H konzolnega nosilca. The goal is to minimize the volume of the structure: minf; f = V((p) \L(k - 4>b) (31) k where wi are the weights, V is the volume, XL is the calculated limit load factor and 4kB the minimum value of the k-th shape parameter. Within the optimization process the most unfavorable imperfection is evaluated according to the procedure described in Fig. 44. The shape optimization and the evaluation of the most unfavorable imperfection was done sequentially. The procedure stops when the evaluated most unfavorable imperfection does not evolve in a change of the calculated limit load factor XL of the optimal structure more than Ae, which was set to a value of 0.01. In Fig. 50 the most unfavorable imperfection shape is presented for the initial shape and the optimal shape considering case A shell thickness (thick sheets). Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 99 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. X 100 20 Z 0 -20 b) 0 -20-10 0 10 20 Fig. 50: The most unfavorable imperfection for the initial (a) and optimal shape (b) (Scale factor=30, Shell thickness A). Slika 50: Najbolj neugodna začetna nepopolnost za začetno in optimalno obliko konstrukcije (Faktor povečave=30, Debelina pločevin A). 100 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. The limit load optimization iteration procedure can be observed in Fig. 51. The limit load curves are plotted for different iterates within optimization. The load-displacement curve considering the optimal shape is plotted in red. Mises stresses are plotted in the illustrations. Red color represents yield stress. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 101 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 1,8 1,6 1,4 1,2 0,8 0,6 0,4 0,2 J__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I__I 10 12 14 16 1 20 22 24 26 28 Displacement [cm] Fig. 51: Limit load shape optimization process. Slika 51: Proces optimizacije oblike v mejnem stanju. 2 1 0 102 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. The evaluated optimal shape for case A shell thickness is illustrated in Fig. 52. The initial geometry at shape parameters being zero (transparent mesh) is sketched behind the optimal shape in order to illustrate the difference between the initial and optimal design. Fig. 52: Initial and optimum shape of H cross-section cantilever geometry for case A shell thickness. Slika 52: Začetna in optimalna oblika H konzolnega nosilca za primer debeline pločevin A. The material is not fully stressed to the yield point in the entire structure, as can be seen in Fig. 53, where Mises stresses are plotted for the optimal structure with a dense FE mesh. The explanation lies behind the way the constraints were chosen. The horizontal sheet has to remain horizontal and the vertical sheet under the horizontal sheet has a minimal height prescribed, which is at least 10% of the initial height of the vertical sheet. Using sheet thicknesses A, the sheets in compression stay compact throughout the optimization procedure and there is no risk of buckling. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 103 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 53: Mises stress at limit state for optimal H cross-section cantilever shape (undeformed). Slika 53: Misesove napetosti v mejnem stanju nosilnosti za H konzolnega nosilec z optimalno obliko (nedeformirana oblika). Fig. 54: Deformation of H cross-section cantilever at limit state (Scale Factor = 1) with Mises stress plotted. Slika 54: Deformacije H konzolnega nosilca pri mejni obtežbi (Faktor povečave = 1) z vrisanimi Misesovimi napetostmi. 104 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. There is only a small amount of lateral displacements in the limit state where extensive rotation occurs. The collapse mechanism can therefore be considered as full plastification with only small amount of plastic buckling. The plastification of the material with increasing load up to the limit state is shown in Fig. 55. Fig. 55: Mises stress for the corresponding load-deformation curve plotted in Fig. 51. Slika 55: Misesove napetosti za krivuljo sila-pomik prikazano na sliki 51. Next, smaller shell thicknesses (case B) were chosen to stimulate the possibility of buckling behavior in order to additionally optimize the shape and to lower the volume of the structure. The initial imperfections calculated according to the method Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 105 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. explained in Chapter 4 used in the initial run and the final run of optimization are drawn in Fig. 56. X 100 20-100 10 20 Fig. 56: The most unfavorable imperfection for the initial (a) and optimal shape (b) (Scale factor=30, Shell thickness B). Slika 56: Najbolj neugodna začetna nepopolnost za začetno in optimalno obliko konstrukcije (Faktor povečave=30, Debelina pločevin B). 106 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. The evaluated optimal shape for shell thickness B is illustrated in Fig. 57 together with initial geometries. The limit load Mises stresses are plotted on the unreformed and deformed mesh in Fig. 58 and Fig. 59, respectively. Fig. 57: Initial and optimum shape of H cross-section cantilever for case B shell thickness. Slika 57: Začetna in optimalna oblika H konzolnega nosilca za primer debeline pločevin B. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 107 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 58: Mises stress at limit state for optimal H cross-section cantilever shape (undeformed). Slika 58: Misesove napetosti v mejnem stanju nosilnosti za H konzolnega nosilec z optimalno obliko (nedeformirana oblika). Fig. 59: Deformation of H cross-section cantilever at limit state (Scale Factor = 1) with Mises stress plotted. Slika 59: Deformacije H konzolnega nosilca pri mejni obtežbi (Faktor povečave = 1) z izrisanimi Misesovimi napetostmi. 108 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. The plastification of the material with increasing load up to the limit state is shown in Fig. 60. Fig. 60: Mises stress for the corresponding load-deformation curve plotted in Fig. 61. Slika 60: Misesove napetosti za obtežno krivuljo prikazano na sliki 61. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 109 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 1,2 -------r --------1---------1----------- i i i i limit statje (A) i i i i i i -----------T----------- ---------T--------- --------------1 ------ -------------1-----------1--------------1 1 1 1 1 limit! state (B) i i i i i i ¦ f / M 1 ' i i i i .___]- ____i_______ —i------------ 1 I, "~r ¦ "h— -? , i : i : i i i j...... i f : 1 j 1 ' ' 1 ' 1 ' -----------J.. — — — Optima (thick s l shape in case a heets) il shape in case B eets) ----- ----- v—j----------- ™ uptime (thin sh — -j- — — "J — - h"-H ------ ----------1 - —i------------ k—-i- -"-h 1 ¦¦¦¦¦¦ i 1 ' ' ' ' ' ' i 1 ' ' ' ' ' ' 1 l-"-H ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ 0,8 0,6 0,4 0,2 0 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 Displacement [cm] Fig. 61: Load-deformation curve for optimized shape in Case A and Case B. Slika 61: Deformacijska krivulja za optimizirano obliko v primeru A in B. In Fig. 61 the load-displacement curve is plotted for both studied optimal structures. In case A the limit load deformation behavior can be described as mostly plastification, while in case B plastic buckling is more pronounced during the optimization procedure. However the optimized structure has been considerably improved regarding buckling behavior of initial designs. The horizontal sheet of the optimal structure was positioned by the optimization algorithm in the most favorable way to support the vertical sheets in the compressed part which is prone to buckle. The optimization algorithm searches for the minimum volume while the limit load must match the prescribed limit load. At the same time the most unfavorable initial imperfections are considered. This combination evolves in a search for a ductile, plastic structure behavior with a small sensitivity to buckling. The result is therefore 1 110 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. a robust structure with minimum weight and small sensitivity to buckling. The plastification zones are spread more widely through the structure which shows the full material usage. The whole optimization procedure can therefore be seen as an efficient tool for economical and safe design. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 111 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 5.4.3 3D single storey steel building In the present example a single storey steel building is being optimized according to the presented method. The initial shape is shown in Fig. 62. Loading conditions and the parameterization illustrated in Fig. 63 and Fig. 64 were used. The outer frames are loaded with the half load described in Fig. 63. Self weight is automatically added by the program. The structure’s basic initial geometry is symmetrical according to the YZ plane. The shell parts of the structure are modeled by elastic-plastic four-node shell elements (Wisniewski, Turska 2000, 2001). The truss parts are modeled by truss elements and have the function of lateral load transmission only. The entire load is added on the frames only. The distance between frames is e = 10m. Y 3000 Z 600 400 200 0 0 X 500 Fig. 62: Initial geometry of single storey steel building. Slika 62: Začetna geometrija enoetažne jeklene hale. 112 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. qr= 49 ql= 10 kN m Fig. 63: Loading conditions and shape parameterization for the inner frame. Slika 63: Obtežbe in prametrizacija oblike za notranji okvir. The shape parameters are used to change the height of the cross-section h in the way: h = h0 + (f)- h0 = h0(1 + (ß) (32) Z Y Second order splines were chosen for interpolation of the boundary shape between the shape parameters. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 113 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. q/= 5 kN m Fig. 64: Loading conditions and shape parameterization for the outer frame. Slika 64: Obtežbe in prametrizacija oblike za zunanji okvir. The procedure described in Fig. 44 is applied to find the optimum shape of the structure. Optimization of the entire structure is computationally demanding. Further on, only the frames are subjected to shape optimization. Therefore the two characteristic frames were chosen for investigation: the outer frame with only half external load applied and the inner frame fully loaded by external forces. The initial and the optimal shape are illustrated for the final run of the optimization algorithm for the inner frame and the outer frame in Fig. 65 and Fig. 66, respectively. Mises stresses are plotted at the limit load deformation state in Fig. 67 and Fig. 68. Z Y 114 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 65: Optimum shape of the inner frame of the steel structure with the initial geometry in the final optimization run. Slika 65: Optimalna oblika notranjega okvira jeklene konstrukcije v zadnjem krogu optimizacijskega procesa. Fig. 66: Optimum shape of the outer frame of the steel structure with the initial geometry in the final optimization run. Slika 66: Optimalna oblika zunanjega okvira jeklene konstrukcije v zadnjem krogu optimizacijskega procesa. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 115 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 67: Deformation of the inner frame at limit state (Scale Factor = 10). Slika 67: Deformacije notranjega okvira pri mejni obtežbi (Faktor povečave = 10). 116 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 68: Deformation of the outer frame at limit state (Scale Factor = 10). Slika 68: Deformacije zunanjega okvira pri mejni obtežbi (Faktor povečave = 10). Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 117 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 69: Mises stress for inner frame for the corresponding load-deformation curve plotted in Fig. 71. Slika 69: Misesove napetosti za notranji okvir za obtežno krivuljo prikazano na sliki 71. 118 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 70: Mises stress for outer frame for the corresponding load-deformation curve plotted in Fig. 71. Slika 70: Misesove napetosti za zunanji okvir za obtežno krivuljo prikazano na sliki 71. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 119 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 1,2 0,8 0,6 0,4 0,2 0 2 4 6 8 10 12 14 16 Displacement [cm] Fig. 71: Load-deformation curves for the inner and the outer frame with the optimized shape. Slika 71: Deformacijski krivulji za notranji in zunanji okvir z optimalno obliko. In Fig. 69 and Fig. 70 the Mises stress development with increasing the load up to the limit state is plotted. The corresponding load-displacement curves are plotted in Fig. 71. Because of the symmetrical shape parametrization and unsymmetrical loading conditions it is impossible for the entire structure to be in plastic state. The optimized structure has a shape which is optimal for the loading conditions considered in the example. Because of the symmetric parametrization the structure is optimized for opposite direction of loads in X direction also. If different ratios of horizontal and vertical loads had to be considered, this could be done with multi objective optimization procedures. Another way is to evaluate an optimum shape for every loading condition and then combine them in a way in which all constraints are still satisfied. 1 0 The optimized structure is shown in Fig. 72. The optimization parameters for the optimized structure are shown in Table 5. 120 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Fig. 72: Optimized shape of single story steel building. Slika 72: Optimizirana oblika enonadstropne jeklene hale. Table 5: Optimization parameter values for the optimized structure. Tabela 5: Vrednosti projektnih parametrov pri optimalni obliki konstrukcije. Frame ob da da da * = "öa ^KAä + * = 0 ä~ := ä + Aä~ Lokalni nivo $(b) = 0 #KAb+$= 0 b := b + Ab Slika 2: Splošna formulacija za direktno analizo tranzientnih, povezanih, nelinearnih problemov. 128 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. V uporabljenem zapisu a predstavlja vektor globalnih parametrov elementa, b je vektor neznanih lokalnih parametrov, definiranih za vsako integracijsko točko (plastične deformacije, spremenljivke utrjevanja, itd.), ap je vektor globalnih parametrov v prejšnjem koraku, bp je neznanih lokalnih parametrov v prejšnjem koraku, ? vektor globalnih enačb in ? vektor lokalnih enačb. 7.2.2 Analitična občutljivostna analiza in polje začetnih občutljivosti Občutljivostna analiza se uporablja za izračun spremembe odziva konstrukcije z ozirom na variacijo projektnih parametrov ? in predstavlja ključen del gradientnih metod optimizacije. Uporaba občutljivostne analize v optimizaciji oblike za mejno stanje konstrukcije zahteva rešitev tranzientnega, povezanega sistema enačb, z upoštevanjem geometrijske in materialne nelinearnosti. Zahtevnost izpeljave izrazov za občutljivostno analizo je bil ključen razlog za uporabo simbolno numeričnega sistema (Korelc 2007a, b). Uporabljena je metoda neposrednega odvajanja (Michaleris, et al. 1994). Zaradi tranzientne narave problema je potrebno občutljivosti izračunati na koncu vsakega obtežnega koraka skozi vso analizo. Ustrezne enačbe so predstavljene na sliki 3. Eden od ključnih problemov uporabe občutljivostne analize v gradientnih metodah optimizacije oblike, je izračun polja začetnih občutljivosti (Korelc, Kristanič 2005). Polje začetnih občutljivosti ( ?X/?? ) opiše spremembo koordinat vozlišč končnih elementov (X) glede na poljubno izbran projektni parameter ?. Medtem ko lahko odvode karakterističnih količin končnega elementa (reziduum, tangentne matrike, itd.) po projektnih parametrih izrazimo s pomočjo avtomatiziranih postopkov (Korelc 2007b), to ne velja za polje začetnih občutljivosti. Glavna ovira se pojavi pri povezavi projektnih parametrov s pozicijo vozlišč. Te povezave ni možno povsem splošno izraziti s standardnimi pristopi generacije mreže končnih elementov, niti s specializiranimi pred-procesorji ali CAD orodji, saj je izbira projektnih parametrov stvar svobodne izbire projektanta. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 129 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Globalni nivo ? (a ,ap,b ,bp, * K Da L>* D(f) Dej) D^ D(f> D?DapD?DbpD? ++ Dap D Dbp D D

$ D& Dap D& Dap D^+ Dap Dej) + Dap Dej) Lokalni nivo $ a ,ap,b ,bp) = 0 $K Db D® D® D(f> Dep D(j) D3>Da Da D

Projektna občutljivost da d4> Slika 4: Potek občutljivostne analize s pomočjo simbolno numeričnega MKE okolja. Opisane postopke je možno uporabiti na problemih poljubne kompleksnosti. Predstavljeni direktna in občutljivostna analiza, sta uporabljeni v določitvi najbolj neugodne začetne nepopolnosti kot tudi v optimizaciji oblike konstrukcij. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 131 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 7.3 Določitev najbolj neugodnih začetnih nepopolnosti Začetna nepopolnost konstrukcij je posledica napak pri izdelavi, ki se jim praktično ni mogoče izogniti. Rezultati nelinearnih numeričnih analiz konstrukcijskih elementov in konstrukcij so lahko v veliki meri odvisni od izbire oblike začetnih nepopolnosti, kar je še posebej izraženo pri obravnavi tankostenskih konstrukcij, občutljivih na spremembo začetne geometrije. Dobro znane razlike med mejno nosilnostjo konstrukcij, izračunane z računalniškimi analizami, ter izmerjene s preizkusi v laboratoriju, je možno zmanjšati z ustreznim upoštevanjem začetnih nepopolnosti. Pri določitvi mejnega stanja konstrukcij v okviru optimizacije oblike, je upoštevanje začetnih nepopolnosti zelo ugodno, saj pripomore k natančnejši določitvi mejne obtežbe in hkrati ugodno vpliva na proces optimizacije, saj se je na tak način mogoče izogniti bifurkacijskim točkam v ravnotežni poti idealnih modelov konstrukcij. Določitev najbolj neugodne začetne nepopolnosti predstavlja zahteven nelinearen optimizacij ski problem, ki je v splošnem v vsakdanji inženirski praksi, praktično nerešljiv v danem časovnem okviru. Z uporabo direktne in občutljivostne analize ter optimizacijskih algoritmov je možno neposredno določiti najbolj neugodno obliko geometrijske nepopolnosti, pri kateri konstrukcija izkaže najnižjo možno nosilnost v okviru obravnavanega problema (Kristanič, Korelc 2008). Pri tem so upoštevane geometrijske, konstrukcijske in materialne nepopolnosti, ki so zajete v obliki ekvivalentnih geometrijskih nepopolnostih ter predpisanih lastnostih materiala. Osnovna ideja predlaganega pristopa je zamenjava nelinearnega optimizacijskega problema z iteracijskim postopkom, v katerem je potrebno reševati le linearne optimizacijske probleme. V okviru metode je možno uporabiti tehnološke pogoje, katerih uporaba je ključnega pomena, saj se je izkazalo, da pri neupoštevanju pravilnih geometrijskih pogojev lahko privede do izračuna nerealistično majhnih mejnih obtežb. Oblika iskane najbolj neugodne začetne nepopolnosti je določena z linearno kombinacijo izbranih baznih oblik: N X = Xp + ^ctjTj , (3) kjer je Xp začetna, popolna geometrija, Nje število izbranih oblik v bazi, aj so neznani parametri oblike in Tj j-ta oblika iz baze oblik. Baza oblik T je lahko izbrana poljubno, vendar mora biti linearno neodvisna. Vsebuje lahko različne nabore 132 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. oblik. To so uklonske oblike (TA), lastni vektorji (TB) začetne togostne matrike K0, lastni vektorji (TC) togostne matrike K0 konstrukcije s spremenjenimi robnimi pogoji, deformacijske oblike (TD) in empirično znane neugodne oblike (TE). Končna baza oblik T je tako: r = rA u rB u rC u rD u rE (4) V okviru metode je iskana tista nepopolna oblika X, pri kateri je mejna nosilnost konstrukcije najnižja. Potek metode je prikazan na sliki 5. Neznani parametri a. , pri katerih bo mejna nosilnost najnižja, so računani iterativno v okviru optimizacijskega procesa. V k-ti iteraciji lahko zapišemo: Xk — Xk-1 + AXk AXk N i = 1 <*i = a^1 + Aaik Xk = N (5) i = 1 kjer je Xk nepopolna oblika, Aaik prirastek parametrov, AXk prirastek nepopolnosti in Xk skupna nepopolnost. Prirastek parametrov Aaik v k-ti iteraciji je dobljen z rešitvijo optimizacijskega problema. Za začetni približek je lahko izbrana kar prva bazna oblika V1: CK, i 0 =0; Aai 0 =- e0 i = 1 max Ti 0 i ^ 1 (6) X0 = Xp + Aa10 T1 Postopek je zaključen, ko je dosežen pogoj I Aaik II < toleranca. V vsaki iteraciji je izvedena nelinearna direktna in občutljivostna analiza konstrukcije z geometrijo Xk. Parametri oblike aik najbolj neugodne oblike Xk v trenutnem koraku, so izračunani s pomočjo optimizacijskega postopka, ki je popolnoma ločen od direktne in občutljivostne analize. Takšen pristop nam omogoča uporabo poljubnega naprednega optimizacijskega algoritma. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 133 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. ZAČETNA GEOMETRIJA Xp INICIALIZACIJA BAZA oblik Y=rAu rBu rCu rDu rE Občutljivostna analiza Analiza do mejnega stanja dxk 8Aak Direktna Ak —y \ analiza Minimizacijski problem k kk min Alk ; Alk = Alk N i=1 N C(Xk,e0)<0; Xk=J]ri i=1 ^lk dAak at1 + Aak •Aa -^ Aak If Aaik > tol k = k + 1 If Aaik < tol X = Xp + E EA«mj NAJBOLJ NEUGODNA OBLIKA Nk N j=1 j=1 ^ m=0 r\ X+X Slika 5: Potek metode določitve najbolj neugodne začetne nepopolnosti. Alternativen in bolj točen pristop bi lahko predstavljalo reševanje polno povezanega problema, vendar zaradi numerične prezahtevnosti trenutno za večje sisteme, še ni mogoč. Polno povezan problem je bil poenostavljen na ta način, da je bil z uporabo 134 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. občutljivosti mejni obtežni faktor nepopolne konstrukcije razvit v Taylorjevo vrsto okoli mejnega obtežnega faktorja nepopolne konstrukcije. Enačbo mejnega obtežnega faktorja lahko zapišemo na sledeči način: Alk \k N Ak =0 + E d\k dAa Agk =0 ¦Aa kjer je ?lk Ak =0 izračunani mejni obtežni faktor v k-ti iteraciji in d\k dAa (7) Ak =0 občutljivost mejnega obtežnega faktorja na optimizacijske parametre v trenutnem koraku. Uporabljen iterativni pristop nam omogoča, da najbolj neugodno nepopolno obliko konstrukcije iščemo na nepopolni konstrukciji. Pri tem se najbolj neugodna oblika iz prejšnje iteracije uporabi za začetno nepopolnost v trenutni iteraciji. S tem je zagotovljena natančnost tudi v primerih velikih nepopolnosti. V vsaki iteraciji je potrebno rešiti minimizacijski problem (8), kjer iščemo takšne Aaik, pri katerih bo Alk minimalen, pod pogojem, da je amplituda oblike, ki jo določajo parametri aik, v predpisanih mejah. Mejne amplitude e0 so določene s principom ekvivalentnih geometrijskih nepopolnosti, ki jih določajo tehnični predpisi (EN 1993 1-5 2004, EN 1993 1-6 2006). Optimizacijski problem lahko zapišemo: min \lk Ak C(Xk,e0) < 0 (8) kjer C(Xk,e0) predstavlja omejitveno funkcijo. Funkcija ?l je linearna, medtem ko je omejitvena funkcija C(Xk,e0) v odvisnosti od zasnove lahko samo ena, izrazito nelinearna funkcija, ali skupek več linearnih funkcij. V prvem primeru je potrebno problem reševati z razširjeno Lagrangevo metodo, v drugem primeru in predvsem pri obravnavi večjih problemov, pa lahko uporabimo metode linearnega programiranja. Z večanjem števila upoštevanih oblik v bazi oblik se napaka diskretizacije najbolj neugodne oblike približuje napaki diskretizacije mreže končnih elementov, kar zagotavlja konvergenco metode z gostenjem mreže končnih elementov in večanjem števila upoštevanih oblik. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 135 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 7.4 Optimizacija oblike konstrukcij V okviru disertacije je bil razvit algoritem za gradientno optimizacijo oblike konstrukcij za mejno stanje konstrukcije. Za izračun gradientov je bila uporabljena analitična občutljivostna analiza. V direktni in občutljivostni analizi so bile upoštevane začetne nepopolnosti, določene z metodo za določitev najbolj neugodne začetne nepopolnosti. Uporabljeno simbolno numerično okolje pri reševanju problemov z gradientnimi metodami optimizacije, omogoča uporabo kombinacije naprednih simbolnih zmožnosti ter hkratne numerične učinkovitosti okolja (Korelc 2002, Korelc 2007a, b). Postopek optimizacije v simbolno numeričnem okolju je prikazan na sliki 6. Z uporabo avtomatskega odvajanja, simultane optimizacije in preverjanja izrazov z uporabo AceGen–a (Korelc 2007b), je pridobljena učinkovita koda končnega elementa. Za simbolno obravnavo v okviru metode končnih elementov je uporabljeno okolje AceFEM (Korelc 2007a). V skladu s postopki, opisanimi v 7.2.2, je izračunano analitično polje začetnih občutljivosti z delom AceFEM-a, z zmožnostjo simbolnega obravnavanja problemov MDriver. Analitično polje začetnih občutljivosti je nato uporabljeno v izračunu občutljivosti z numeričnim delom AceFEM-a (CDriver). Odvodi namenske funkcije in pogojev po projektnih spremenljivkah so nato posredovani optimizacijskemu algoritmu v okviru okolja Mathematica (Wolfram 2008). Kot je bilo že poudarjeno, so natančne gradientne informacije, izračunane na podlagi analitičnega polja začetnih občutljivosti, odločilnega pomena za konvergenco optimizacijskega algoritma, še posebej, kadar obravnavamo probleme z upoštevanjem geometrijske in materialne nelinearnosti. AceGen in AceFEM delujeta znotraj okolja Mathematica, kar je zelo priročno, saj ni potrebno uporabljati vmesnikov za povezavo z drugimi okolji za optimizacijo. Poleg tega Mathematica nudi paleto modernih optimizacijskih algoritmov, ki se neprestano nadgrajujejo in jih je možno uporabiti direktno z AceFEM-ovo direktno in občutljivostno analizo. Splošen potek optimizacijske metode je predstavljen na sliki 6, splošna optimizacijska zanka na sliki 7. 136 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Okolje generiranja kode AceGen Izvorna koda za končni element CDRIVER C Numerična koda MDRIVER Mathematica Koda v simbolni obliki Opis modela z metodo končnih elementov Izračun najbolj neugodne začetne nepopolnosti (Poglavje 7.3) Direktna in občutljivostna analiza (Poglavje 7.2) Analitično polje začetnih nepopolnosti (Poglavje 7.2.2) Okolje za optimizacijo Matematično programiranje Optimalna oblika Slika 6: Optimizacija s pomočjo simbolno-numeričnega okolja. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 137 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Simbolno numerično okolje za KE AceFEM Opis modela z metodo končnih elementov I Preveri potrebo po ponovnem izračunu ¦+ —>-z------,—M začetnih nepopolnosti Izračun najbolj neugodne začetne nepopolnosti X Konvergenčni kriteriji NISO zadoščeni; Spremeni projektne spremenljivke Numerične vrednosti za parametre oblike I Numerična generacija mreže KE z CDRIVER X = XP(cj)) + X Direktna in občutljivostna analiza Parametri oblike v simbolni obliki Simbolna generacija mreže KE z MDRIVER XP(?) T Analitično polje začetnih občutljivosti ?X ?XP(?) = ?? dcf> Okolje za optimizacijo Matematično programiranje Konvergenčni kriteriji doseženi Optimalna oblika Slika 7: Optimizacijska zanka s pomočjo simbolno-numeričnega okolja. 138 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Znotraj optimizacijskega procesa (slika 7) je geometrija konstrukcije posodobljena zaporedoma v dveh iteracijskih zankah. V notranji zanki se v skladu z optimizacijskim algoritmom spreminjajo projektne spremenljivke (f), ki določajo geometrijo X((f)). V zunanji iteracijski zanki se geometrija konstrukcije spremeni zaradi spremenjenih najbolj neugodnih začetnih nepopolnosti X. Metoda za določitev najbolj neugodne nepopolnosti je predstavljena na sliki 5. Geometrijo konstrukcije lahko zapišemo kot: X = XP(0) + X (9) kjer je XP((j)) idealna geometrija konstrukcije in X celoten vektor nepopolnosti kot opisano v poglavju 7.3. Nepopolna oblika mora biti upoštevana v direktni in občutljivostni analizi, medtem ko se za izračun polja začetnih občutljivosti lahko uporabi idealna geometrija konstrukcije, saj velja: dX dXP(d>i) + X dXP(d>) = ---^---= —-p± (10) 8(f) 8(f) 8(f) Polje začetnih občutljivosti ni odvisno od celotnega vektorja začetnih nepopolnosti X in se ne spreminja pri spremembi projektnih spremenljivk (f). Zaradi tega ga je potrebno izračunati le enkrat. Najbolj neugodno začetno nepopolnost je potrebno izračunati vsakič, ko se geometrija konstrukcije v okviru optimizacije oblike spremeni do določene mere. Skupaj s spremenjeno geometrijo se spreminja tudi oblika najbolj neugodne nepopolnosti. Izračun najbolj neugodne nepopolnosti je računsko zahteven, zato je ugodno, če izračuna ni potrebno izvajati v vsaki iteraciji optimizacije oblike. Izračun je izveden samo takrat, ko bi spremenjena oblika začetnih nepopolnosti spremenila mejni obtežni faktor za več kot AL , ali bi se spremenilo mejno stanje. Ugotovljeno je bilo, da se oblika neugodnih začetnih nepopolnosti ne spreminja bistveno pri majhnih spremembah projektnih spremenljivk. Poleg tega je konvergenca optimizacijskega algoritma bolje ohranjena, če se oblika začetnih nepopolnosti ne spreminja. V predlaganem pristopu se izračun najbolj neugodnih nepopolnosti izvrši izmenoma z optimizacijo oblike, kar se je izkazalo kot uspešen pristop pri iskanju optimalne končne oblike konstrukcije. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 139 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. 7.5 Zaključek Cilj disertacije je bil razviti optimizacijsko metodo, katero bi bilo možno uporabiti v projektiranju konstrukcij. Večina modernih tehničnih standardov pri projektiranju konstrukcij zahteva, da je nosilnost konstrukcije takšna, da bo v svoji življenjski dobi z dano verjetnostjo prenesla vse predvidene obtežbe in da ne bo presegla projektiranih mejnih stanj. Standardni postopki optimizacije oblike konstrukcij v splošnem ne upoštevajo mejnih stanj konstrukcij. Običajno se konstrukcije optimizira z ozirom na njihov volumen, ceno, ali kakšno drugo lastnost, kot recimo lastna frekvenca ali togost. Izračun pravilne mejne obtežbe konstrukcije je zahtevna naloga, zato se pri klasični optimizaciji oblike konstrukcij uporabljajo kriteriji napetosti, prvega nastopa elastičnega uklona ali nastop bifurkacijske točke. Ti pristopi v splošnem ne omogočajo izračuna optimalne oblike za dejansko mejno stanje konstrukcije. V predstavljenem pristopu je prikazana optimizacija oblike za mejno stanje, v okviru katere je možno natančno določiti mejno obtežbo konstrukcije, s hkratnim upoštevanjem začetnih nepopolnosti ter materialnih in geometrijskih nelinearnosti ter jo uporabiti kot pogoj v procesu optimizacije. V metodo še ni implementirana možnost eksplicitnega upoštevanja zaostalih napetosti. Upoštevati jih je možno implicitno s pomočjo metode nadomestnih geometrijskih nepopolnosti. V skladu s tehničnimi standardi je na tak način ob predpostavki, da so zajeti vsi relevantni fenomeni, s pomočjo optimizacijskih algoritmov možno projektirati konstrukcije. V okviru disertacije je bil uporabljen simbolno numerični pristop k optimizaciji oblike. Ta omogoča poljubno simbolno parametrizacijo konstrukcije, kar ni mogoče pri klasičnih pristopih k optimizaciji oblike. Simbolna oblika omogoča analitičen izračun polja začetnih občutljivosti, s čimer je možno izvesti natančno občutljivostno analizo, ki je temeljnega pomena za uspeh uporabljenih gradientnih metod optimizacije. V primeru konstrukcij, ki so občutljive na spremembo začetne geometrije, določitev mejne obtežbe ter optimizacija oblike v mejnem stanju, zahtevata pravilno upoštevanje začetnih nepopolnosti. Velikost in oblika začetnih nepopolnosti lahko ima velik vpliv na odziv konstrukcije ter njeno mejno stanje. Nadalje je možno, da je ob neupoštevanju začetnih nepopolnosti rezultat optimizacije napačen, saj se kot rezultati lahko pojavijo zelo lahke konstrukcije, ki so zelo občutljive na pojav uklona. Kljub številnim raziskavam, eksperimentalnemu delu ter numeričnim analizam konstrukcij, občutljivih na spremembo geometrije, med strokovnjaki še vedno ni 140 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. enotnega mnenja, na kakšen način izračunati mejno stanje, kar je možno pripisati številnim težavam, ki se pojavijo. Nepopolnosti konstrukcije niso znane v naprej, zato je bila razvita metoda za določitev najbolj neugodne nepopolnosti nanašajoč se na mejno stanje konstrukcije. Metoda je implementirana kot interni ločen optimizacijski algoritem znotraj globalne optimizacije oblike konstrukcije. Prikazano je, da je s pomočjo geometrijsko in materialno nelinearne analize nepopolnih konstrukcij kombinirane z optimizacijo, možno direktno določiti začetno nepopolno obliko, pri kateri konstrukcija izkaže najmanjšo možno mejno nosilnost. Metoda ni omejena na linearno obtežno pot in ne na majhne nepopolnosti ter omogoča vpeljavo različnih tehnoloških pogojev glede same oblike. Običajno je potrebno za določitev neugodne kombinacije začetnih nepopolnosti opraviti številne analize. Kljub naporu je težko določiti stopnjo, do katere se zmanjša mejna obtežba konstrukcije z upoštevanjem tako dobljenih začetnih nepopolnosti in z izračunano najbolj neugodno nepopolnostjo. Na podlagi rezultatov predstavljenega pristopa je težko okarakterizirati določene konstrukcije z določenimi tipi nepopolnosti. Vsaka sprememba v debelini, obliki ali obtežbi lahko povzroči drastično spremembo najbolj neugodne oblike. Pri kompleksnih konstrukcijah, občutljivih na spremembo oblike, kjer je težko intuitivno določiti neugodne nepopolnosti in kjer empirično pridobljenih neugodnih oblik ni na voljo, je uporaba metode za določitev najbolj neugodnih nepopolnosti nujna. Zaradi nepredvidljivosti oblik nepopolnosti, tehnični standardi predlagajo uporabo empiričnih metod določanja začetnih nepopolnosti v numeričnih analizah. Mejne obtežbe konstrukcij, izračunane s pomočjo predstavljene metode, so se izkazale za manjše, kot bi jih dobili s pomočjo kombiniranja različnih oblik po priporočilih standardov ali s pristopi, ki temeljijo na Koiterjevi asimptotični teoriji ali parametričnimi študijami. Te metode lahko pripeljejo do preveč optimističnih rezultatov. Po drugi strani je verjetnost, da bi realna konstrukcija imela najbolj neugodno začetno obliko, zelo majhna. Kljub temu je informacija, katera nepopolnost je najbolj neugodna, zelo pomembna, kadar numerično simuliramo obnašanje konstrukcij. V tem smislu je mogoče povzeti, da je uporaba celostne metode določitve najbolj neugodnih nepopolnosti pri geometrijsko in materialno nelinearni analizi, nepogrešljiva. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 141 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Pomembnost možnosti izračuna najbolj neugodnih začetnih nepopolnosti se nadalje izkaže pri uporabi le teh v optimizaciji oblike v mejnem stanju nosilnosti. Skozi celoten proces optimizacije je konsistentno upoštevana polna geometrijska in materialna nelinearnost, kar omogoča efektiven in robusten način projektiranja konstrukcij. Pristop, ki upošteva dejansko mejno obtežbo ter uporabo najbolj neugodnih nepopolnosti, povzroči, da je v okviru optimizacije oblike iskana takšna konstrukcija, ki izkazuje plastično in duktilno obnašanje z zmanjšano nevarnostjo nastanka nestabilnosti ter po drugi strani minimalnim volumnom konstrukcije. Na tak način je rezultat predlagane metode optimizacije robustna konstrukcija z minimalno težo ter minimalno možnostjo uklona pri danih pogojih. Projektiranje konstrukcij z integracijo metod optimizacije oblike za mejno stanje ter določitve najbolj neugodne začetne nepopolnosti, predstavlja nov in napreden pristop k projektiranju konstrukcij. Uspešno kombinacijo analize mejne obtežbe in optimizacijskih metod omogoča uporaba simbolno numeričnega okolja za analizo konstrukcij. Z upoštevanjem vseh pomembnih fenomenov, lahko predlagan pristop predstavlja način projektiranja varnih in ekonomičnih konstrukcij, ter s tem boljšo alternativo klasičnemu projektiranju konstrukcij na mejna stanja. Pomembnejše izvirne prispevke k tehničnim znanostim predstavljajo naslednji razviti postopki: - Razvoj metode za določitev najbolj neugodne začetne geometrijske nepopolnosti. - Razvoj algoritmov za optimizacijo oblike konstrukcij v mejnem stanju z upoštevanjem najbolj neugodnih nepopolnosti. - Izračun analitičnega polja začetnih nepopolnosti in točna občutljivostna analiza s pomočjo simbolno-numeričnega okolja. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 143 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. References Arora, J. S. 2004. Introduction to Optimum Design. San Diego, Academic Press: 728 p. Arora, J. S., Lee, T. H., Cardoso, J. B. 1992. Structural Shape Sensitivity Analysis -Relationship between Material Derivative and Control Volume Approaches. Aiaa Journal, 30, 6: 1638-1648. Bažant, Z. P., Cedolin, L. 2003. Stability of structures : elastic, inelastic, fracture, and damage theories. Mineola, N.Y., Dover Publications: 1011 p. Bletzinger, K. U., Ramm, E. 2001. Structural optimization and form finding of light weight structures. Computers & Structures, 79, 22-25: 2053-2062. Bonnans, J. F., Gilbert, J. C., Lemaréchal, C., Sagastizábal, C. A. 2006. Numerical Optimization Theoretical and Practical Aspects. New York, Springer: 490 p. Camprubi, N., Bischoff, M., Bletzinger, K. U. 2004. Shape optimization of shells and locking. Computers & Structures, 82, 29-30: 2551-2561. Cappello, F., Mancuso, A. 2003. A genetic algorithm for combined topology and shape optimisations. Computer-Aided Design, 35, 8: 761-769. Chang, K. H., Choi, K. K., Tsai, C. S., Chen, C. J., Choi, B. S., Yu, X. M. 1995. Design Sensitivity Analysis and Optimization Tool (Dso) for Shape Design Applications. Computing Systems in Engineering, 6, 2: 151-175. Chang, K. H., Tang, P. S. 2001. Integration of design and manufacturing for structural shape optimization. Advances in Engineering Software, 32, 7: 555-567. Che, J., Tang, S. 2008. Research on integrated optimization design of hypersonic cruise vehicle. Aerospace Science and Technology, In Press, Corrected Proof. Choi, K. K., Chang, K.-H. 1994. A study of design velocity field computation for shape optimal design. Finite Elements in Analysis and Design, 15, 4: 317-341. Choi, K. K., Kim, N. H. 2005a. Structural sensitivity analysis and optimization 1, Linear systems. New York, Springer Science+Business Media: 446 p. Choi, K. K., Kim, N. H. 2005b. Structural sensitivity analysis and optimization 2, Nonlinear systems and applications. New York, Springer Science+Business Media: 447-772 p. Crisfield, M. A. 1996. Non-Linear Finite Element Analysis of Solids and Structures. New York, Wiley: 362 p. Crisfield, M. A. 1997. Non-Linear Finite Element Analysis of Solids and Structures, Volume 2, Advanced Topics. New York, Wiley: 508 p. 144 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Deml, M., Wunderlich, W. 1997. Direct evaluation of the 'worst' imperfection shape in shell buckling. Comput Method Appl M, 149, 1-4: 201-222. Dinkler, D., Pontow, J. 2006. A model to evaluate dynamic stability of imperfection sensitive shells. Comput Mech, 37, 6: 523-529. El Damatty, A. A., Nassef, A. O. 2001. A finite element optimization technique to determine critical imperfections of shell structures. Struct Multidiscip O, 23, 1: 75-87. Elishakoff, I. 2000. Uncertain buckling: its past, present and future. Int J Solids Struct, 37, 46-47: 6869-6889. EN 1090/2. 2007. Execution of steel structures and aluminium structures – Part 2: Technical requirements for steel structures. Brussels, European Committee for Standardization EN 1993 1-5. 2004. Eurocode 3: Design of steel structures, Part 1.5: Plated structural elements., Brussels, European Committee for Standardization EN 1993 1-6. 2006. Eurocode 3: Design of steel structures, Part 1.6: Strength and stability of Shell Structures., Brussels, European Committee for Standardization Ewert, E., Schweizerhof, K., Vielsack, P. 2006. Measures to judge the sensitivity of thin-walled shells concerning stability under different loading conditions. Comput Mech, 37, 6: 507-522. Garcia, M. J., Gonzalez, C. A. 2004. Shape optimisation of continuum structures via evolution strategies and fixed grid finite element analysis. Struct Multidiscip O, 26, 1: 92-98. Godoy, L. A. 2000. Theory of Elastic Stability: Analysis and Sensitivity. Philadelphia., Taylor and Francis. Haftka, R. T., Grandhi, R. V. 1986. Structural Shape Optimization - a Survey. Comput Method Appl M, 57, 1: 91-106. Hansen, J. S., Liu, Z. S., Olhoff, N. 2001. Shape sensitivity analysis using a fixed basis function finite element approach. Struct Multidiscip O, 21, 3: 177-195. Hardee, E., Chang, K.-H., Tu, J., Choi, K. K., Grindeanu, I., Yu, X. 1999. A CAD-based design parameterization for shape optimization of elastic solids. Advances in Engineering Software, 30, 3: 185-199. Ho, D. 1974. Buckling load of non-linear systems with multiple eigenvalues. International Journal of Non-Linear Mechanics, 6, 5: 649-662. Hörnlein H.R.E.M. 2000. Effiziente semi–analytische gradienten-berechnung in der strukturoptimierung. Z. Angew. Math. Mech., 81: 669-670. Hughes, T. J. R. 2000. Finite Element Method - Linear Static and Dynamic Finite Element Analysis Mineola, NY, Dover. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 145 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Jang, G. W., Kim, Y. Y. 2005. Sensitivity analysis for fixed-grid shape optimization by using oblique boundary curve approximation. Int J Solids Struct, 42, 11-12: 3591-3609. Johansson, B., Maquoi, R., Sedlacek, G., Muller, C., Beg, D. 2007. Commentary and worked examples to EN 1993-1-5 "Plated structural elements". Luxembourg, Office for Official Publications of the European Communities: 226 p. Kegl, M. 2000. Shape optimal design of structures: an efficient shape representation concept. Int J Numer Meth Eng, 49, 12: 1571-1588. Koiter, W. T. 1945. Over de stabiliteit van het elastisch evenwicht. doctor, Techische Hooge School, Delft. English translation : Edward Riks (1969), The stability of elastic equilibrium. Stanford University. Korelc, J. 1997. Automatic generation of finite-element code by simultaneous optimization of expressions. Theoretical Computer Science, 187: 231-248. Korelc, J. 2002. Multi-language and multi-environment generation of nonlinear finite element codes. Eng Comput-Germany, 18, 4: 312-327. Korelc, J. 2007a. AceFEM, Mathematica finite element environment. Ljubljana, University of Ljubljana, Faculty of Civil and Geodetic Engineering. http://www.fgg.uni-lj.si/Symech/. Korelc, J. 2007b. AceGEN, Multi-Language, Multi-Environment numerical code generation. Ljubljana, University of Ljubljana, Faculty of Civil and Geodetic Engineering. http://www.fgg.uni-lj.si/Symech/. Korelc, J., Kristanič, N. 2005. Evaluation of Design Velocity Field by Direct Differentiation of Symbolically Parameterized Mesh. In: Plasticity: fundamentals and applications: proceedings of the eighth international conference on computational plasticity, E. Onate, (ed.), CIMNE, Barcelona, Spain, 380-383. Kristanič, N., Korelc, J. 2008. Optimization method for the determination of the most unfavorable imperfection of structures. Comput. Mech., DOI:10.1007/s00466-008-0288-9. Lanzo, A. D. 2000. A Koiter's perturbation strategy for the imperfection sensitivity analysis of thin-walled structures with residual stresses. Thin Wall Struct, 37, 1: 77-95. Lanzo, A. D., Garcea, G. 1996. Koiter's analysis of thin-walled structures by a finite element approach. Int J Numer Meth Eng, 39, 17: 3007-3031. Lemaitre, J., Chaboche, J.-L. 1990. Mechanics of Solid Materials. Cambridge, Cambridge University Press: 556 p. Li, Q., Steven, G. P., Querin, O. M., Xie, Y. M. 1999. Evolutionary shape optimization for stress minimization. Mechanics Research Communications, 26, 6: 657-664. Li, W., Li, Q., Steven, G. P., Xie, Y. M. 2005. An evolutionary shape optimization for elastic contact problems subject to multiple load cases. Comput Method Appl M, 194, 30-33: 3394-3415. 146 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Marsden, J. E., Hughes, T. R. J. 1994. Mathematical foundations of elasticity. New York, Dover Publications. Maute, K., Nikbay, M., Farhat, C. 2000. Analytically based sensitivity analysis and optimization of nonlinear aeroelastic systems. 8th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, AIAA, Long Beach, CA, 2000–4825. Maute, K., Schwarz, S., Ramm, E. 1999. Structural optimization - The interaction between form and mechanics. Zeitschrift Fur Angewandte Mathematik Und Mechanik, 79, 10: 651-673. Meske, R., Lauber, B., Schnack, E. 2006. A new optimality criteria method for shape optimization of natural frequency problems. Struct Multidiscip O, 31, 4: 295-310. Michaleris, P., Tortorelli, D. A., Vidal, C. A. 1994. Tangent Operators and Design Sensitivity Formulations for Transient Nonlinear Coupled Problems with Applications to Elastoplasticity. Int J Numer Meth Eng, 37, 14: 2471-2499. Nocedal, J., Wright, S. 2006. Numerical Optimization. New York, Springer: 664 p. O.C. Zienkiewicz, J.S. Campbell. 1973. Shape optimization and sequential linear programming. In: R.H. Gallagher, O. C. Zienkiewicz, (eds.). Optimum Structural Design. New York, Wiley: 109-126. OptiStruct. 2008. OptiStruct manual. Altair Engineering Inc. www.altair.com. Oral, S. 1996. An improved semianalytical method for sensitivity analysis. Structural Optimization, 11, 1: 67-69. Papadopoulos, V., Papadrakakis, M. 2005. The effect of material and thickness variability on the buckling load of shells with random initial imperfections. Comput Method Appl M, 194, 12-16: 1405-1426. Ramm, E., Mehlhorn, G. 1991. On shape finding methods and ultimate load analyses of reinforced concrete shells. Engineering Structures, 13, 2: 178-198. Rong, J. H., Jiang, J. S., Xie, Y. M. 2007. Evolutionary structural topology optimization for continuum structures with structural size and topology variables. Advances in Structural Engineering, 10, 6: 681-695. Ryu, C. H., Lee, Y. S. 2007. A study on Ranked Bidirectional Evolutionary Structural Optimization (R-BESO) method for fully stressed structure design based on displacement sensitivity. Journal of Mechanical Science and Technology, 21, 12: 1994-2004. Saliba, R., Padra, C., Venere, M. J., Taroco, E., Feijoo, R. A. 2005. Adaptivity in linear elastic fracture mechanics based on shape sensitivity analysis. Comput Method Appl M, 194, 34-35: 3582-3606. Samareh, J. A. 1999. A novel shape parametrization approach. In: NASA/TM-1999-209116. Hampton, Virginia. Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. 147 Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Schenk, C. A., Schueller, G. I. 2003. Buckling analysis of cylindrical shells with random geometric imperfections. International Journal of Non-Linear Mechanics, 38, 7: 1119-1132. Schmidt, H. 2000. Stability of steel shell structures - General Report. J Constr Steel Res, 55, 1-3: 159-181. Schmit, L. A. 1960. Structural Design by Systematic Synthesis. In. Proc. 2nd Conference on Electronic Computation. New York, American Society of Civil Engineers: 105-132. Schneider, W. 2006. Stimulating equivalent geometric imperfections for the numerical buckling strength verification of axially compressed cylindrical steel shells. Comput Mech, 37, 6: 530-536. Schneider, W., Brede, A. 2005. Consistent equivalent geometric imperfections for the numerical buckling strength verification of cylindrical shells under uniform external pressure. Thin Wall Struct, 43, 2: 175-188. Schneider, W., Tirnmel, I., Hohn, K. 2005. The conception of quasi-collapse-affine imperfections: A new approach to unfavourable imperfections of thin-walled shell structures. Thin Wall Struct, 43, 8: 1202-1224. Schranz, C., Krenn, B., Mang, H. A. 2006. Conversion from imperfection-sensitive into imperfection-insensitive elastic structures. II: Numerical investigation. Comput Method Appl M, 195, 13-16: 1458-1479. Schwarz, S., Maute, K., Ramm, E. 2001. Topology and shape optimization for elastoplastic structural response. Comput Method Appl M, 190, 15-17: 2135-2155. Song, C. Y., Teng, J. G., Rotter, J. M. 2004. Imperfection sensitivity of thin elastic cylindrical shells subject to partial axial compression. Int J Solids Struct, 41, 24-25: 7155-7180. Steven, G. P., Li, Q., Xie, Y. M. 2002. Multicriteria optimization that minimizes maximum stress and maximizes stiffness. Computers & Structures, 80, 27-30: 2433-2448. Tang, P. S., Chang, K. H. 2001. Integration of topology and shape optimization for design of structural components. Struct Multidiscip O, 22, 1: 65-82. Uysal, H., Gul, R., Uzman, U. 2007. Optimum shape design of shell structures. Engineering Structures, 29, 1: 80-87. van Keulen, F., Haftka, R. T., Kim, N. H. 2005. Review of options for structural design sensitivity analysis. Part 1: Linear systems. Comput Method Appl M, 194, 30-33: 3213-3243. Wisniewski, K., Turska, E. 2000. Kinematics of finite rotation shells with in-plane twist parameter. Comput Method Appl M, 190, 8-10: 1117-1135. Wisniewski, K., Turska, E. 2001. Warping and in-plane twist parameters in kinematics of finite rotation shells. Comput Method Appl M, 190, 43-44: 5739-5758. 148 Kristanič, N. 2008. Limit State Design Using Exact Sensitivity Analysis and Shape Optimization. Doctoral thesis. University of Ljubljana, Faculty of Civil and Geodetic Engineering. Wolfram. 2008. Mathematica 6.0, The user guide. Wolfram Research, Inc. Wriggers, P., Simo, J. C. 1990. A General Procedure for the Direct Computation of Turning and Bifurcation Points. Int J Numer Meth Eng, 30, 1: 155-176. Wunderlich, W., Albertin, U. 2000. Analysis and load carrying behaviour of imperfection sensitive shells. Int J Numer Meth Eng, 47, 1-3: 255-273. Wunderlich, W., Albertin, U. 2002. Buckling behaviour of imperfect spherical shells. International Journal of Non-Linear Mechanics, 37, 4-5: 589-604. Zienkiewicz, O. C., Taylor, R. L. 2000a. The finite element method. Vol. 2, Solid mechanics. Oxford Butterworth-Heinemann: 459 p. Zienkiewicz, O. C., Taylor, R. L. 2000b. The fnite element method. Volume 1. The Basis. Oxford Butterworth-Heinemann: 689 p.