# Full text of "Boundary element methods in structural shape synthesis"

## See other formats

BOUNDARY ELEMENT METHODS IN STRUCTURAL SHAPE SYNTHESIS By JUNHAUR JIH A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1989 To Mei-Ming and my parents ACKNOWLEDGEMENTS I would like to express my deep appreciation to Dr. Hajela, chairman of my supervisory committee, for his encouragement, guidance, and valuable suggestions throughout my graduate program and the preparation of this manuscript. He has always treated me as a friend rather than a student, for which I am deeply indebted to him. My appreciation is extended as well to the other members of the committee, Dr. Malvern, Dr. Sun, Dr. Sankar, and Dr. Hoit. Thanks also go to Ms. Bloebaum, Ph.D. student in this department, who helped me with the writing of this work. Computational support received under a grant from the Pittsburgh Supercomputer Center is gratefully acknowledged. Finally, I wish to express my heartfelt gratitude to my wife Mei-Ming for her unselfish love and support during the completion of this study. iii TABLE OF CONTENTS PAGE ACKNOWLEDGEMENTS iii ABSTRACT vi CHAPTER 1 INTRODUCTION 1 1.1 Motivation and Objectives 1 1.2 Scope of Present Work 5 2 OVERVIEW OF OPTIMUM SHAPE DESIGN 7 2 . 1 Literature Survey 7 2.1.1 Finite Element Method for Optimum Shape Design 7 2.1.2 Boundary Element Method for Optimum Shape Design 12 2.2 Formulation of Shape Optimization Problem 2.3 Overview of Mathematical Programming Method 3 BOUNDARY ELEMENT METHODS FOR STRUCTURAL ANALYSIS 28 3.1 Introduction 28 3.2 Application in Torsional Shaft Problems 2 9 3.3 Application in Elasticity Problems 37 3.4 Numerical Results 42 4 BOUNDARY ELEMENT METHOD FOR OPTIMUM SHAPE SYNTHESIS 52 4.1 Selection of Design Variables 52 4.2 Optimum Shape Design for Torsional Shaft 4.3 Optimum Shape Design in 2-D Elasticity Problems 4.4 Discussions 69 IV 5 SENSITIVITY ANALYSIS FOR BOUNDARY ELEMENT APPLICATIONS 71 5.1 An Overview of Sensitivity Analysis Methods 71 5.2 Semi-analytical Approach 75 6 GRID REFINEMENT AND GRID ADAPTATION IN THE BOUNDARY ELEMENT METHOD 83 6.1 Introduction 83 6.2 Grid Refinement Technique 84 6.3 Adaptive Scheme Based on Variational Principles 85 6.4 Grid Adaptation Obtained by Nonlinear Programming Methods 91 6.5 Examples of Grid Refinement and Adaptation 95 7 AN INTEGRATED APPROACH FOR SHAPE OPTIMIZATION USING THE BOUNDARY ELEMENT METHOD 113 7.1 Introduction 113 7.2 Integrated Approach Using the Boundary Element Method 115 8 CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER RESEARCH 125 APPENDIX 129 REFERENCES 133 BIOGRAPHICAL SKETCH 139 v Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy BOUNDARY ELEMENT METHODS IN STRUCTURAL SHAPE SYNTHESIS By Junhaur Jih August, 1989 Chairman: Prabhat Hajela Major Department: Aerospace Engineering, Mechanics and Engineering Science The purpose of this study is to present an efficient implementation for optimal structural shape synthesis which is based on a combination of the boundary element method (BEM) with nonlinear programming (NLP) optimization algorithms. Shape optimization involving two-dimensional elastic, torsional shaft problems, and small-strain elastostatics problems, is introduced in this study. An approach for grid refinement and grid adaptation for the BEM is developed for application in plane elasticity problems. To prevent distortion of the original domain, master nodes are introduced to control the geometry of the structure. Two strategies for the solution of the adaptive schemes are investigated, one based on a variational approach and another employing nonlinear programming methods. The use of nonlinear programming methods in grid distribution allows vi greater flexibility in selecting the form and relative importance of criteria employed for grid adaptation and refinement. By utilizing the characteristics of both the BEM and the LU-decomposition methods, an efficient use of a semi- analytical sensitivity analysis for optimum design is obtained. This study also examines the use of BEM in optimal shape synthesis from the standpoint of an integrated formulation. A set of linear eguations obtained from the boundary element analysis is treated as equality constraints in the mathematical programming optimization problem. Several numerical examples are presented for illustrating the effectiveness and convenience of using BEM with these techniques. vii CHAPTER 1 INTRODUCTION 1 . 1 Motivation and Objectives In most structural design problems, the geometric configuration of the structure is usually fixed, and structural member dimensions, such as cross-sectional areas of bar elements, cross-sectional properties of beam elements, and thickness of membrane and plate elements, are changed to obtain the optimum design. It is obvious, however, that an additional mass reduction can be obtained by changing the structural geometry in addition to changing structural member dimensions. The problem of optimum shape design includes design variables that characterize either the shape of an exterior boundary of the structure, or the location and shape of an interior cutout. An extensive survey of attempts in shape optimization is presented in Haftka and Grandhi [1]. There are two major components that must be suitably integrated when attempting optimum structural design. The first relates to choosing a mathematical programming algorithm to achieve an optimal objective function and simultaneous satisfaction of a set of design constraints. 1 2 The second is the selection of an analytical tool for determining the structure behavioral response needed in the computational procedure. In shape design problems where design variables control the geometry of the structure, additional degrees of difficulty are introduced. At the very outset, there is generally an increase in computational resource requirement, as it is more difficult to obtain good sensitivity information with respect to shape design variables, than with respect to member sizing variables. This problem is further compounded if the finite element method is used as the analytical tool. Since the finite element domain is changed continuously in the redesign, undesirable element distortion is often introduced, making it difficult to ensure that the accuracy of the finite element analysis remains adequate throughout the design processes. Although adaptive mesh generation procedures such as r- and ^“finite element methods have been developed to reduce the error due to element distortion [2], it is cumbersome to perform a new mesh generation for each iteration of the shape design. The boundary element method (BEM) has been successfully employed in various engineering disciplines [3], and offers several advantages over the finite element method (FEM) [4], Among these, the more significant include a replacement of the FEM mesh generation for the entire domain by a boundary discretization in the BEM, a higher accuracy in sensitivity 3 evaluation along the boundary, and simpler preparation of input data for the problem. Furthermore, the level of experience necessary in using BEM is generally not as high as is reguired in the use of FEM for mesh generation and selection of optimum element sizes. However, some judgement is still reguired in using the BEM for optimum shape design. In order to reduce the number of grid points in the boundary element analysis and to obtain a more accurate numerical solution, a technigue using a combination of grid refinement and grid adaptation is explored in this study. In the conventional application of adaptive technigues, the variational approach is usually used to redistribute the analytical nodes. In the present study, a more general nonlinear mathematical programming approach is developed for this purpose. Such a strategy allows consideration of multiple criteria in the grid distribution problem. The sensitivity of structural response to changes in the design variables provides valuable information for design, and is extensively used in most efficient nonlinear programming algorithms. Such sensitivity may be obtained by either finite difference methods or an explicit analytical differentiation of the response eguations. The former is easier to implement but has an associated computational cost that may be undesirable in large problems. The analytical approach is sometimes difficult to implement due to the 4 implicit nature of the response equations. In the present work, a semi-analytical approach based on the characteristics of the boundary integral equation is investigated. The concept of integrated approach in structural design is also studied in this work. This approach was first presented by Schmit and Fox [5,6] in a FEM based structural design. The basic idea behind the effort was to treat the behavioral response equations (finite element equilibrium equations) as additional equality constraints in the optimization problem. The behavior response variables were added to the list of design variables. The equality constraints and the inequality design constraints were subsequently handled by a penalty function approach in an integrated design problem. The advantage of the approach resided in the fact that no explicit solution of the response equations was necessary, and the optimum solution simultaneously met the desired objective, satisfaction of the response constraints, and a solution to the analysis problem. Ineffective methods of solving equality constrained problems did not allow this method to be efficient till very recently [7]. In the present work, a generalized reduced gradient approach is used for the solution of the equality constrained problem in such an integrated formulation. 5 1 . 2 Scope of Present Work The literature regarding optimum shape design problems is presented in Chapter 2. An overview of mathematical programming methods and the corresponding mathematical formulations are also presented in the same chapter. In Chapter 3, the boundary integral equations of two- dimensional elasticity and torsional problems are derived. The applications and numerical accuracy of the BEM for both cases are illustrated. In Chapter 4, the selection of shape design variables is outlined. Several shape optimization problems are attempted by using the BEM and a nonlinear programming technique for analysis and optimization, respectively. The subject of sensitivity analysis is explored in Chapter 5. The form of the boundary integral equation is used advantageously in deriving a semi-analytical approach for sensitivity analysis. The grid refinement and grid adaptation techniques are investigated in Chapter 6. Both the variational approach and the nonlinear programming method are employed to achieve grid adaptation. Numerical examples of shape optimization, combined with these techniques, are also presented. In Chapter 7, the concept of the integrated approach is extended to shape optimization problems, wherein the BEM analysis equations are represented as a set of equality constraints in the optimization problems. 6 Finally, conclusions drawn from this study and recommendations for future research are presented in Chapter 8 . CHAPTER 2 OVERVIEW OF OPTIMUM SHAPE DESIGN 2 . 1 Literature Survey With the increased availability of high-speed digital computers since the 1960s, structural optimization has occupied an important place in the field of engineering design. Recent studies in structural optimization have extended the traditional structural member sizing problem to the determination of optimal shapes for desired structural response. The technique of structural shape optimization has been demonstrated successfully in various publications. A general review of papers published during the last decade can be found in reference [l]. In this chapter, the review of work related to optimum shape design is classified according to the two analytical tools known as the finite element method and the boundary element method. 2.1.1 Finite Element Method for Optimum Shape Design The first numerical treatment of shape optimization was proposed by Zienkiewicz and Campbell [8] in 1973. The importance and applications of shape optimal design was 7 8 emphasized in this article. A finite element model was formulated to analyze the shape optimization problem, and the coordinates of boundary nodal points were treated as design variables. The physical problem in this study was the design of dam structures using isoparametric elements. Sequential linear programming techniques were employed for optimization purposes, with the derivatives of structural responses obtained by an implicit finite difference method. Ramakrishnan and Francavilla [9] characterized this shape optimization problem in more detail. A hybrid penalty function method was employed to obtain optimum design for dam structures. The structural domain was discretized into two-dimensional isoparametric elements, and the coefficients of the element polynomials were chosen as design variables. In 1975, Francavilla, Ramakrishnan and Zienkiewicz [10] first applied the finite element method in a fillet shape optimization problem with the objective of minimizing the stress concentrations in the structural domain. Kristensen and Madsen [11] formulated a class of shape optimal design problems for planar solids. They used orthogonal polynomials to locate the boundary of the structural body and treated the coefficients of these polynomials as design variables. A finite element approach was employed to obtain the structural response and the derivatives of stress with respect to design variables. Chung and Haug [12] used an approach based on the 9 calculus of variations, combined with adjoint differential equations, to obtain the sensitivities necessary to perform the two-dimensional shape optimal design. The design objective was to minimize structural weight with constraints expressed by von Mises yield stress along the boundary. A generalized steepest descent algorithm was implemented for determination of the optimum shape. To simplify the optimization procedure, Dems and Mroz [13] modified the design variables to a set of prescribed shape functions and a set of shape parameters; an optimality criterion was derived and satisfied by using an iterative numerical algorithm. The two-dimensional torque arm design problem was first solved by Botkin [14], resulting in the basic concept for shape optimization of plate and shell structures. Later, Queau and Trompette [15] used an automated mesh generator for several two-dimensional problems. Shape optimization of three-dimensional structural components was first proposed by Imam [16]. In order to reduce the number of design variables in the three- dimensional shape optimization problem, Imam suggested the use of low order curves and surfaces to model the domain. The isoparametric representation of the surfaces and the numerical superposition of shapes to simplify the shape representation was also proposed. A feasible directions algorithm of nonlinear programming was used to obtain the 10 minimum mass. The dilemma of finite element distortion in each new shape design during the optimization procedure was also discussed in this work. The initial configuration of the finite element mesh may result in very distorted elements for the new shape. Therefore, even the best mesh design for the initial shape is usually not appropriate for intermediate or final designs. As a result, the analysis may be inaccurate in the computations for stresses and displacements. A solution- based adaptive mesh refinement scheme was studied by Bennett and Botkin [17], in which the strain density gradients were used to identify the regions which require mesh refinement. This concept was also successfully extended to three- dimensional structural shape optimization [18]. Kikuchi et al. [19] treated the automated mesh generation as an integral part of shape optimization. Two general grid adaptive methods, the r- and h- methods, applicable in shape optimum design are discussed in detail in their study. They emphasize the strong dependence of the final shape of design boundaries on the mesh generation and conclude that a good adaptive mesh refinement strategy can avoid the jagged shape produced by using the coordinates of boundary nodes as design variables. Furthermore, an accurate analysis for the discrete model may improve design convergence and reduce the total computational time. Using a mixed variational formulation, Rodrigues [20] 11 derived the optimality condition for which the maximum value of the von Mises equivalent stress in a two-dimensional linear elastic body is minimized. Choi and Seong [21] proposed domain integrals for shape design sensitivity analysis of built-up structures. To avoid the inaccuracies of finite element results on the boundary, the sensitivity information was expressed as domain integrals rather than boundary integrals. In such an approach, design sensitivity formulas for each standard component type (rod, beam, plate, etc.) are obtained, and one can define a library of structural components that may be assembled to form a built- up structure. Another important problem in shape optimization is the determination of the optimal cross section of a torsional shaft. Banichuk and Karihaloo [22] present a solution for an optimal cross section of an isotropic-solid bar for minimum weight and subject to constraints on the torsional rigidity and bending stiffness. Dems [23] used a domain perturbation technique to derive the optimality criterion for an isotropic, hollow bar which was designed for maximum torsional rigidity with a constant area and a fixed inner boundary. Adali [24,25] concluded that the cross section of an anisotropic hollow bar with a high torsional rigidity was given by a domain bounded by two geometrically similar ellipses . The torsional shape optimization problem was also 12 solved by using an approach in which the material derivatives were employed to compute the design sensitivity [26] . A finite element formulation was used for the analysis. An extension of this investigation [27] involved consideration of internal-node movement associated with the regridding process on the shape design sensitivity. 2.1.2 Boundary Element Method for Optimum Shape Design The first numerical treatment using the BEM in structural shape optimization was proposed by Mota Soares, Rodroigues, Oliveira Faria and Haug [28]. The optimization of the geometry of solid and hollow shafts for torsional requirements was formulated in terms of the BEM. Various analysis models, based on constant, linear and quadratic boundary elements were used requiring only a discretization along the boundary. The design objective was to obtain the optimal cross section of a shaft with a given area, so as to maximize its torsional stiffness. The discrete boundary element models are much more efficient and robust than the corresponding finite element discretizations. This is largely attributable to the higher accuracy of analysis response, and hence the more accurate sensitivity information obtained. Furthermore, there is no need to compute the state variables in the domain, because the maximum stress is at the boundary. 13 Mota Soares, Rodrogues and Choi [29,30] have also developed the shape optimal design of two-dimensional elastic structures based on linear and quadratic boundary elements. The corresponding nonlinear programming problem is solved by a linearization method. The design objective in this problem requires a minimization of the compliance of the structure, subject to an area constraint. All degrees of freedom of the model are defined only along the boundary. Thus, there is no need for calculating displacements and stresses in the domain. The determination of the optimal shape for two- and three-dimensional elasticity problems, based on the boundary element method, has recently been developed by Burczynski and Adamczyk [31]. A maximal stiffness criterion, together with the limitation of the volume of a body, have been used. The optimality conditions were derived for an optimal boundary. An iterative process, based on finite differences and the Newton-Raphson method, is used to solve a set of nonlinear algebraic equations. Choi and Kwak have developed a general method for the shape design sensitivity analysis as applicable to elasticity problems using the boundary integral equation formulation [32]. The material derivative concept and adjoint variable method were applied to obtain sensitivity formulae for the stress constraint functions. Accuracy of the design sensitivity analysis obtained from the material 14 derivative approach was compared with the conventional finite difference method. The design of a fillet and an elastic ring were used to illustrate the effectiveness of the material derivative concept. A smooth characteristic function was found to be better than a plateau function for the localization of the stress constraint. Sandgren and Wu [33] used the BEM with substructuring concepts for two- and three-dimensional structural shape optimization. Substructuring in the boundary element method allows the shape changes of each substructure to be isolated for separate analysis. The control points on B-spline curves and surfaces, which were introduced to describe the shape of the domain were selected as design variables. The generalized reduced gradient method was used to solve the nonlinear programming problem. 2 . 2 Formulation of Shape Optimization Problem A general formulation of the shape optimization problem usually involves minimizing either the structural mass or stress by selecting shape design variables and simultaneously satisfying prescribed geometry and response constraints (or behavior constraints) . The mathematical formulation of a structural shape optimization problem can be stated as follows: 15 Find xi i=l , 2 , . . . ,n Minimize F(xjJ Subject to gj (Xi) < 0 j =1 , 2 , ,m h k( x i) = 0 k=l,2, . . . ,q Xi 1 < xi < xi u i=l , 2 , . . • , n ( 2 . 1 ) where xi is a vector of shape design variables; F is the objective function; gj and h^ are inequality and equality constraints on the response quantities; x^ and xi u are the lower bound and upper bound of each design variable. The structural geometry is mainly described by the definition of shape design variables. The dimensionality of the design space is dictated by the number of design variables. From the standpoint of the efficiency of the optimization process and from manufacturing considerations, design variables are usually chosen as control parameters, most typically, the coefficients of polynomials or spline curves in the shape representation. The side constraints on each design variable serve as the geometry constraints which represent technological limitations on the design variables and reflect fabrication or analysis validity considerations. The equality constraints govern the structural response associated with 16 the loading conditions, which are usually obtained from the analysis tool. The inequality constraints represent the restrictions on structural response (eg, stresses, displacements, or natural vibration frequencies) . In the shape optimization problem, the von Mises stress function is usually used to calculate the equivalent stress so that each boundary point has one stress constraint. Once the optimal shape design problem is formulated, there are two tasks that must be considered. The first relates to choosing the optimizer to achieve an optimal objective function and satisfy a set of constraints. The second is in the selection of the tool to implement the structural analysis. Between the optimizer and analyzer, a grid generator is used to facilitate data preparation for the structure analysis. The grid generator, a preprocessor of the analytical tool, converts the design variables obtained from the optimizer to structural geometry variables and generates the nodes for the analyzer. During the optimization process, the optimizer proposes trial design variables to the analyzer through the grid generator, and the analyzer calculates the response constraints and feeds all analysis results back to the optimizer. Figure 2.1 illustrates the relation of the optimizer, the analyzer and grid generator. 17 Optimizer Figure 2 . 1 The relation among optimizer, grid generator, and analyzer. 18 Sensitivity analysis deals with the calculation of the derivatives of response with respect to the design variables. During a typical optimization seguence, the sensitivity analysis accounts for most of the CPU time. It is precisely for this reason that the issue of design sensitivity analysis in optimal shape design has been studied so extensively. The choice of proper sensitivity analysis method depends on the characteristics of both the design model and the analytical tool. A detailed introduction of sensitivity analysis will be discussed in Chapter 5 . 2 . 3 Overview of Mathematical Programming Method The purpose of mathematical programming techniques is to provide a computer tool to aid the designer in performing various optimal design tasks. The use of these techniques in engineering design problems offers a logical approach to design automation. The category of the mathematical programming techniques of optimization can be classified into two distinct types: (A) linear mathematical programming, (B) nonlinear mathematical programming. In the linear programming (LP) problem, all constraints and the objective function are linear functions of the design variables. The simplex method is available for solving the linear programming problem, as is the recently introduced 19 and still controversial approach proposed by Karmarkar [34]. In these problems, the optimum always occurs on the boundary of the feasible domain and the exact global optimum is reached in a finite number of steps. Although most engineering problems are not of the LP form, some nonlinear programming (NLP) problems can be approximated by linear expressions as in the linear sequential programming (LSP) method. Also the search direction in the NLP method often requires the solution of an auxiliary LP problem, as in the usable-feasible search direction method. In most of the structural optimization problems, both the objective function and the constraints might be highly nonlinear functions or even be implicit functions of the design variables. The following discussions are focused on describing the characteristics of the nonlinear optimization algorithms used in the present work. In the application of NLP method in engineering optimization problems, convergence to the global optimum is seldom assured. The design spaces in such problems typically nonconvex and several local optima exist in the region of search. Another pitfall is that a poorly scaled design space obtained for a particular choice of design variable may result in extremely slow convergence of the optimization algorithm. One way to increase the probability of locating the global optimum is to start the optimization process from several initial points. 20 Nonlinear programming algorithms can be classified into distinct categories of unconstrained and constrained optimization methods. The necessary conditions for a local minimum for unconstrained problems is that the gradient of the objective function must vanish at that point, and furthermore, the Hessian matrix must be positive definite. Here the Hessian matrix is defined as the matrix of second partial derivatives of the objective function with respect to the design variables. In a constrained optimization problem, the necessary conditions for a design point (x*) to be considered a local optimum are defined in terms of the Kuhn-Tucker conditions of optimality. These are written as follows : x* is feasible ^ j 9j ( x ) = 0 / A j > 0 m 1 VF (x*) + S Aj vg-j (X*) + s A Vh k (x*) = 0 j=l k=l K m Aj > 0 , A m+k unrestricted in sign where the symbols, j, k, m, gj and h k , are defined in equations (2.1); the A' s are the Lagrange multipliers which are uniquely determined from the gradients of the objective function and active constraints; VF, Vgj , and Vh k are the 21 first derivatives of objective function, inequality constraints, and equality constraints, respectively, with respect to the design variables. The available algorithms for unconstrained NLP can be categorized on the basis of the function and function derivative information required in computing the search direction. Examples of this category is as follows, (1) Zero-order methods: Random search method [35]. Powell's method [36]. (2) First-order methods: Steepest descent method [35]. Conjugate direction method [37]. Variable metric methods [38]. (3) Second-order methods: Newton's method [39,40]. The above methods each have their own advantages and drawbacks. Generally, the more sophisticated methods converge more rapidly and provide more accurate results. However, the efforts involved in computing the searching direction, particularly those related to obtaining gradient information, must be taken into consideration when determining the effectiveness of a given method. The available algorithms for constrained NLP are divided into two groups. The first, indirect methods, converts a constrained problem to an unconstrained problem by appending the violated constraints to the objective by means of a 22 penalty function. In such an approach the optimum is achieved by sequentially solving an unconstrained minimization problem with several values of the penalty parameters. This approach is generally separated into the exterior penalty function method, the interior penalty function method, and the extended penalty function method. These methods are discussed in more detail by Fiacco and McCormick [41]. The augmented Lagrange multiplier method [42] utilizes the Kuhn-Tucker conditions and includes the Lagrange multipliers into the exterior penalty function, improving the convergence of the optimization process. The second approach is a class of direct methods which solve the constrained problem directly. The methods of sequential linear programming (SLP) method [43], feasible direction method [44], generalized reduced gradient (GRG) method [45], and sequential quadratic programming method [46] belong to this category. When the evaluations of the objective function and sensitivity are costly, these approaches are preferable. The following discussion focuses on those methods which are used in the present research. The sequential linear programming method linearizes the nonlinear problem with a first-order Taylor expansion as follows . 23 Minimize F(x) = F(x°) + VF(x°) • fix S.T. gj (x) = gj (x°) + vgj (x°) • fix < 0 h k (x) = h k (x°) + Vh k (x°) • fix = 0 Xi 1 < Xi+6Xi < Xi U where fix = x-x° (2.2) The objective function and the constraints are linearized at the initial design x°. By constructing move limits for each design variable, the optimum of the nonlinear problem is obtained within a few linear programming cycles. The method converges well for fully constrained problems. However, if the optimum is located at a point where only a few constraints are active, the convergence pattern sometimes exhibits serious oscillations. The usable feasible search direction method belongs to the category of direct methods. A FORTRAN program, CONMIN [47], based on this method is widely applied in engineering problems. In this approach, a direction vector {S} is determined which satisfies both feasible and usable conditions, expressed mathematically as {S) T {Vgj} < 0 and {S} t {VF} < 0, respectively. A physical interpretation of 24 such a search direction is that it must point in a direction which allows no constraint violation and also result in a decrease of the objective function (for nonlinear problems) . The search direction vector { S } is obtained as a solution of the linear programming problem stated as follows: Max. p Subject to {S} T (Vgj) + p <0 {S} T {VF} + p < 0 -{1} < { S > < (1) (2.3) Here p is a positive scalar quantity, larger value of which makes (S) T {VF) more negative, and forces the search direction (S) to point in a direction in which the objection function decreases more rapidly; 6 j is a positive push-off factor which drives the search direction away from the constrained boundary and into the feasible domain. The advantages of this method are that only a linear programming problem is solved in determining the search direction, and furthermore, only the active constraints need to be considered in calculating the gradient information, usually a time-consuming process. However, this method is only well suited for fundamentally inequality constrained 25 problems. If equality constraints are introduced into the problems, this method will not be applicable in the form presented above. The generalized reduced gradient method is developed to solve nonlinear problems in the sense that the simplex algorithm exists to solve the linear problem. This is accomplished by adding slack variables to inequality constraints. The general form of eqn. (2.1) is thus stated as , z n-q independent variables Find x = { } y m+q dependent variables Minimize F(x) = F(z,y) Subject to V>k( x ) = 0 k=l , 2 , . . . , m+q xi 1 < xi < xi u i=l , 2 , . . . , n x j +n * 0 j=l,2, . . . ,m (2.4) Differentiating the objective and constraint functions in eqn. (2.4), the following forms are obtained. dF(x) = V Z F (x) • dz + VyF (x) • dy (2.5) d ^k( x ) = V z V>j c (x)-dz + VyV>k(x)dy = 0 k=l,m+q ( 2 . 6 ) 26 Here, the subscripts z and y identify the gradients with respect to the independent and dependent variables, respectively. By solving for dy in eqn. (2.6) and substituting dy into eqn. (2.5), the generalized reduced gradient, V r , is obtained as follows. dF (x) v r F T (x) = = V z F t (x) - v y F T (x) • [V y V>f 1 • V Z V> dz (2.7) The reduced gradient defines the rate of change of the objective function with respect to the decision variable with the state variable adjusted to maintain feasibility. The problem of eqn. (2.4) is then converted to, Min. f(z) ; z = [ z ± , z 2 , . . . , z Q ] T S.T. z (1) < z < z (u) (2.8) where the state variables y have been eliminated from the original problem by using the constraints ^ m (z,y)=0 to solve for y in terms of z. A necessary condition for the existence of a minimum of an unconstrained nonlinear function is that the elements of the gradient vanish. 27 V r f(x)i = > 0 if zi =zj/ 1 ^ < 0 if zi =zj/ u ^ = 0 otherwise i 1,2,. ..,Q (2.9) The use of GRG method effectively solves the problems with a few inequality constraints and several equality constraints. If there are a large number of equality constraints, the computer storage becomes prohibitive. CHAPTER 3 BOUNDARY ELEMENT METHODS FOR STRUCTURAL ANALYSIS 3 . 1 Introduction The boundary element method is an alternative approach to obtaining a solution to a given set of governing equations in continuum mechanics. The technique consists of the transformation of the partial differential equation describing the behavior of the unknowns inside and on the boundary of the domain into an integral equation relating only boundary values. The geometry of the boundary is then discretized into several elements to formulate a set of linear system equations. Numerical solutions of these equations with prescribed boundary values allow a computation of the unknown responses along the boundary. If values at the internal points are required, they are calculated from the boundary information. Since all numerical approximations take place only at the boundaries, the dimensionality of the resulting linear system of equations is considerably smaller than that obtained in the finite element method. Hence, the amount of computation required is substantially reduced, as is the 28 29 amount of labor involved in defining the geometry of the problem. These favorable characteristics contribute towards the recognition of the boundary element method as an efficient analysis tool for the optimization process. The interest in engineering applications of the BEM increased significantly in the latter part of the last decade, and there have been several recent contributions to the literature in this field [48]. In this chapter, a brief introduction and application of the boundary element method is introduced. Two-dimensional elastic, torsional shaft problems, including singly- and doubly-connected prismatic bars, and a small-strain elastostatics problem is solved by the BEM, and these implementations are presented here for completeness . The BEM for the potential problem is based on Green's formula, eqn.(3.1), which allows the formulation of certain boundary value problems as integral equations [48]. 3 • 2 Application in Torsional Shaft Problems (w- du • u) dr (3.1) Since the governing equation of the torsional problem is an 30 e llip"tical partial differential equation, it may be treated as a potential problem. Let u denote the stress function of a torsional elastic problem. This stress function should satisfy the Poisson's equation as follows: Lu = v z u = -b in 0 (3.2) Here L is the differential operator, V 2 is the Laplace operator, and b is a prescribed constant ( b=2 for the problem of torsion of an elastic bar) . The boundary conditions corresponding to this torsional problem are of two types, as shown in Figure 3.1. These types are (a) essential condition, u = u on r^, where u is a prescribed value of stress function; (b) natural condition, q = — = q on T 2 , and q is a prescribed value of normal derivative over the boundary. The total boundary is denoted by = r 2 + T 3 « Assume that w is a fundamental solution of Lw = -6 (x,x Q ) (3.3) where S(x,x 0 ) is the Dirac delta function which has the property 31 Figure 3.1 Definition of the boundary element method for potential type problem. 32 Jn u(x) 8 (x,x Q ) dx u(x Q ) , x D e o 1/2 u(x Q ) , x Q e r (3.4) The fundamental solution of the above problem is of the form, w = i 1 4 r -1 2n for 3-D for 2-D (3.5) where r= | x - x D | . Applying Green's identity, eqn. (3.1), and making use of the fundamental solution in eqn. (3.4) allows one to write an expression for u at a point Xp as follows. f 3u 3w f u ( x p) = Jr (w--^ — -u) dr + b-w dQ (3.6) It is observed that the fundamental solution, eqn. (3.5), is valid for any point in the domain n. However, if x Q approaches the boundary r , then the solution at x Q becomes singular # In order to formulate the problem by the BEM approach, the singularity must be specially handled by mathematical analysis. For instance, in a 3-D problem, if the boundary is assumed to be smooth, one may cut out a small hemispherical region of the boundary to analyze the 33 singular point (see Figure 3.2). The boundary integration around the singular point Xp is calculated as r aw f i iim ( J e U-— dr ) = lira ( -L u — — - dr ) = -1/2 u e >0 ° n e >0 4,re lira ( L € >0 J i au dr ) 4?rr an 0 where dr = r sin 6 ■ dc/> ■ d?r Thus, if Xp is on the boundary, equation (3.6) becomes f aw r au r 1 / 2 • u (x p ) + J r u dr = r -w dr + r b-w do an an (3.7) To implement the BEM on a digital computer, the boundary of the structure geometry is discretized into n segments. Equation (3.7) can then be written in a discrete form as, n n l/2-Ui(Xp) + 2 Uj • Hij = 2 qj • Gij + Bi (3.8) j=l j=l where, f aw H ij = Jr dr an G^j — J* r w dr 34 Figure 3.2 Boundary surface assumed to be hemispherical for integration purpose. 35 B i = Jn b-w do The first two terms of the left hand side in equation (3.8) are combined in one term, still represented by H^j . The system equation can now be written as, n n S Hi j • Uj = S Gij • qj + Bi (3.9) j=l j=l or can be expressed in matrix form as follows. [H] • (U) = [G] • (Q) + { B) (3.10) The boundary conditions described above are only applicable to the singly-connected torsional bar. When the problem involves torsion of a prismatic bar with a doubly- connected domain (Figure 3.3), the boundary conditions will change as follows [23], u = 0 on r u = u Q on r Q 3u ds = 2 Do on r Q (3.11) 3 n where D 0 is the area enclosed within the interior boundary, and K, the torsional rigidity, is expressed as 36 Figure 3 . 3 A multiply-connected domain. 37 K = 2 I n u dft + 2 u Q • D 0 (3.12) For the multiply-connected domain, one unknown constant u D and one additional equation (3.11) are introduced into the linear equation (3.10). The dimensions of the matrices, H and Q, are then modified as (n+l)*(n+l), and the dimensions of the matrices U, Q, and B are changed to (n+1) *1. The boundary element method for elasticity problems is based on Somigliana's identity [3], and is given by. where uj (£ ) is the displacement at point £ in j direction. The quantities pj (x) and bj (x) are the actual state of the traction and the body force at point x, respectively; G ij(£/X) and F^(£,x) represent the Kelvin solution for displacement and traction at £ in the j direction due to a 3 . 3 Application in Elasticity Problems (3.13) 38 un it force applied at x in the i direction. These functions can be expressed as follows, 1 8ttG(1-j/) t (3-4./) ln(— £) -fiij + dr dr ax. dXjJ (3.14) F ij<*' x > = 4?r ( 1 -v ) r L an r ar ,,,,,, , . ar dr t — ((i- 2 -)' ‘ij + 2 ) ,, _ . ,3r „ ar , , (1-2l/) ' ( ax. ' n j ax. ' n i } ] 1 D (3.15) where, r = | x - e I - If the field point £ is close to the boundary, one arrives at an expression involving only the boundary unknowns. Since the singularity point is introduced during the boundary integration, the integral equation (3.13) assumes the form, uj U) Jr [ P(x) Gij (e,x) - F ij (£,x) uj (x) ] dr (x) b-Gij (£,x) do (x) (3.16) where C. . ID can be expressed as follows, 39 C. . = S . . - I . . 13 ID ID (3.17) Here, 1^^ is defined in terms of 9 ]_ and 9 2 (Figure 3.4), and is written as, 1 8 tt ( l~i/ ) 4 (1-u) (n + 9 2 -9 i ) + sin 20 x - sin20 z COS20 Z - 0032^ cos 20 z “ COS 20 -L 4 (1-u) (n+9 2 ~9 1 ) +sin20 z “ sin 20 x where C-^j is the coefficient that depends on the geometry of the boundary at point £, i.e., equation (3.16), or may be determined from rigid body motion [48]. The numerical solution of eqn. (3.13) is facilitated by discretizing the boundary into N elements. The values of displacement u and traction p over each element are approximated by using the interpolation function <t > , expressed in terms of the nodal values u n and p n as follows. u = <t> T u n (3.18) P = n P (3.19) Substitution of expressions for u and p into equations 40 Figure 3.4 Definition of domain boundary orientation as defined by angles 8 ± and 82 - 41 (3.16) and summing over all the N segments, results in an alternative matrix notation as follows. where [H] and [G] contain the results of all the integrations. These integrals may be evaluated by any suitable method such as the Gauss quadrature integration scheme. The system of equations (3.20) can be rearranged in such a way that all unknowns are written on the left hand side in a vector of unknowns (x), such that where [A] is an N*N matrix, and (C) is a vector of prescribed boundary values. Once the boundary values are calculated, one can obtain stresses at any internal point by using the stress and displacement relation. [H] • (U> = [G] • (P) + (B) (3.20) [A] (x) = (C) (3.21) (3.22) where A and G are Lame' constants. Differentiating equation (3.13) and substitution into the equation (3.22), the stresses at the internal points are obtained as follows. 42 2Gi/ " ij = Jr ‘ TX ^lk 3G ik 3G ik 6.. — + G ( — + — 1 — )} Pk dr J a xi sx-i axi r 2 Gj/ { JQ 1 - 2 is aG lk aG ik aG ik S-M — ~ + G ( — + — 1 — ) } b k dQ J a xi axj axi I 2Gv aF 3F. aF. { 6., + G( — + — 1 — ) } U k dr c i- 2 i^ J a xi axj axj (3.23) In general, the stress tensor on the boundary is also important for structural shape design problem. The derivation of the stress tensor along the boundary is shown in the Appendix. 3 . 4 Numerical Results Numerical implementations of the boundary element method are illustrated by an elliptic torsional shaft and a cylindrical shell in a state of plane strain. Since the performance of structural analysis may affect convergence during the shape optimization, both numerical results are compared with analytical solutions. Computer programs for the solution of the isotropic torsional problem and two- 43 dimensional elastostatics problem are generated in this study. These programs are available for both constant and linear elements along the boundary. As mentioned earlier, the program requires that only the boundary of the domain be discretized. There is no requirement of a special assembler subroutine as in the finite element method. Furthermore, not only is the data preparation much simpler but also the order of the resulting linear system equations is considerably lower than that in the finite element approach. However, the sparseness and the symmetry of the influence matrix [A] are generally lost. Figure 3.5 (a) shows the solution of the stress function for a torsional prismatic shaft of an elliptic cross section. Figure 3.5 (b) displays the stress distribution along the elliptic boundary. The comparison of boundary stresses between computational and analytical solutions are listed in table 3.1. Table 3 . 1 Comparison of stresses along the boundary X y r/G/3 (analytical) r/G/3 ( computational ) error % 5 0 2.6471 2.6493 0.083 4 . 05 1.76 3 .3632 3 . 3749 0.347 2 . 5 2 . 6 4 . 0435 4.0583 0.366 0.52 2 . 98 4 .3964 4 .4116 0.346 44 Figure 3.5 (a) Solution of stress function for elliptic torsional shaft. stress Figure 3.5 (b) Shear stress distribution along the boundary. 45 The numerical results of a prismatic bar under torsion with a doubly-connected domain, is also examined in this section. Assume that A and B define the major and minor lengths of the outer elliptic boundary, respectively, and the inner elliptic cutout is represented by the major axis a and a minor axis b. Two examples, figure 3.6, for different major and minor lengths are used to illustrate the numerical results. In the tables, u Q represents the inner unknown constant in equation (3.11), K is torsional rigidity as expressed in equation (3.12), and the maximum and minimum stresses along the boundary are also listed. The torsional rigidities for this case are used to compare the numerical solutions with analytical solutions. Figure 3.7 shows a cylindrical shell subjected to the internal pressure of P a = 100 units and external pressure of Pb = “25 units acting on the surfaces at r a = 3 and rb = 8 , respectively. The problem can be assumed to be one of plane strain. Variations in both the computational and analytical solutions of and t along the radial direction are shown in figure 3.8 (a) and (b) , respectively. It is apparent that the accuracy of the numerically computed stress at internal points near the boundary is lost. Since the internal stress is determined from the derivatives of the integral representation for displacement, (equation (3.23)), the term containing r’ 2 obtained from the derivation of Equation (3.23) in the integral stress stress 46 Figure 3.6 Stress distributation along the multiply- connected elliptic domain of the torsioal shaft. 47 Figure 3 . 6 Continued 48 Figure 3. P b Description of cylindrical shell subjected to internal pressure P a and external pressure Pfc,. 49 Figure 3.8 (a) Stress distribution r „ versus cylindrical shell. xx 50 Figure 3.8 (b) Stress distribution r versus cylindrical shell. ^ 51 results in a singularity for small values of r. Due to such a strong singularity, integration points cannot be allowed to approach the boundary. Stress at such a point can be obtained from the stress computed either by using the interpolation method or by decreasing the singularity by one order [ 49 ] . CHAPTER 4 BOUNDARY ELEMENT METHOD FOR OPTIMUM SHAPE SYNTHESIS 4 . 1 Selection of Shape Design Variables The choice of shape design variables has a rather critical influence in determining the character of the boundary in a shape optimization problem. A proper selection of shape design variables not only substantially influences the characteristics of the design space in the nonlinear programming, but also avoids undesired or impractical shape design. Figure 4.1 shows the final shape obtained in a the fillet design problem involving the use of a boundary element method. The jagged shape mainly results from improper grid distribution along the boundary. This often happens when coordinates of the boundary nodes are defined as design variables. Thus, either from the standpoint of manufacturing and practicality of a design, or from the view of computational savings, the selection of design shape variables plays a very important role in the shape optimization problem. Several representations of shape design variables have been proposed to describe the structural boundary in a 52 53 Design generated when node locations are used as design variables Figure 4 . 1 54 design domain. Except for the use of node coordinates to represent the boundary, various shape representations of structural boundaries have been summarized in reference [1]. The polynomial representation has been commonly used to represent boundaries in optimum shape design. After taking the coefficients of the polynomial as design variables, the design shape is then characterized by a set of polynomials [50,51]. Dems [23] describes the design shape by a set of prescribed shape functions and a set of shape parameters. The optimization procedure is reduced to the determination of these parameters. Since the use of high-order polynomials usually results in an oscillatory shape, spline representations composed of several low-order elements which are combined to maximize smoothness are used to eliminate this problem [49]. Results obtained in [52] also demonstrate the effectiveness of spline representations in producing better accuracy in the sensitivity information in comparison to an approach in which piecewise liner representation of the boundary is used. The design element concept was first introduced by Imam [16]. The design elements are described by a set of key nodes that control the geometry. Braibant and Fleury [53] used Bezier and B-spine blending functions to represent the boundary of design elements. Many selections of shape design variables are described above. However, the literature does not elaborate on how to 55 choose the order of the polynomial when such a polynomial representation is used. Furthermore, in some problems it may be advantageous to remove material from within the domain. In such situations, it is not always obvious what kind of shape is preferred, what shape such an interior cutout should have, or where such a cutout should be located. To obtain a better understanding of these problems, a guideline for choosing proper shape design variables based on the domain stress contours is studied in this section. This study is limited to a class of design problems where the analysis domain is linear elastostatics . At the very outset of generating a shape design modeling process, the boundary and loading conditions are generally specified. Assuming an initial design domain that satisfies both conditions, the von Mises equivalent stress in this domain can be calculated by any analytical tool. In the present study, the boundary element method was used for this task. An equivalent stress contour plot is a useful representation to gain a better understanding of the load distribution in the structural domain. By connecting the equal equivalent stress ratio points in the structural domain, where the stress ratio is the von Mises stress divided by the allowable stress, the stress distribution on the structural domain can be easily gauged, as shown in figure 4.2. If each equivalent stress ratio line is 56 distribution in the 57 considered to be a stress path, it is apparent that a small amount of stress flow passes through the low stress regions. When reviewing the problem from the standpoint of optimum shape design, the low stress regions may justifiably be treated as regions from which material can be removed. However, there still exist problems, including those related to the existence of low stress regions, and the impracticality of the boundary shapes if all such regions were removed in the shape design problem. As discussed earlier, the eguivalent stress contour plot provides valuable information for optimum shape design. To use the information effectively, several basic shape elements, such as line elements, circle elements, triangular elements and parabolic elements, were defined and made available for optimum shape design. Each element was handled by control parameters, such as locations of the two end points for a line element, the radius or the location of the center point for a circle, or one end point slope and the location of the two end points for a parabolic element. Guided by the equivalent stress distribution, these basic elements were introduced into the low stress regions, and new structural boundaries defined by these elements. The sizes of these elements were determined by the control parameters. In other words, the control parameters were considered to be design variables during the shape optimization procedure. 58 An example of stress contour based representations of shape design boundaries will be shown in Chapter 6, where this approach was combined with grid adaptation and grid refinement techniques. Depending on the structural geometry and the loading conditions, the numerical results of the boundary element method are seriously affected by grid distribution. These problems are most frequently encountered in structures where regions of stress concentration exist. In the following section, several shape design examples are presented using an approach that combined nonlinear programming optimization methods with the boundary element method. These shape design problems include torsional shafts and structural components that are characterized by a state of plane stress, and for which previous finite element based solutions are available in the literature. Design variables chosen in this study are nodal coordinates and coefficients of low-order polynomials. Since the grid adaptation and grid refinement techniques are not applied in these specific examples, the grids used in two-dimensional elasticity problems are evenly and densely located along the structural boundaries. 4 . 2 Optimum Shape Design for Torsional Shaft The first example in this sequence of optimum designs 59 involves the determination of the boundary of a prismatic bar for maximum torsional rigidity, subject to a constraint of a fixed cross sectional area of 50 units. An additional requirement that the moment of inertia about the horizontal axis be greater than 350 units was also investigated in this example. The definitions of design variables, geometry descriptions, and the initial design for this shaft are shown in figure 4.3. Six design variables are used to obtain the optimum design, using both a linear piecewise approximation to the constraints and objective function, and a fully nonlinear treatment of these functions. The latter approach directly converges to the optimum design, while the former approach results in oscillation near the optimum point. The final optimum shapes without and with the constraint on moment of inertia about a horizontal axis, are shown in figures 4.4 and 4.5, respectively. The optimum values of all design variables, the torsional rigidity, cross sectional area, and bending inertia about the horizontal axis for the initial and final designs are also shown in each figure. The second example is an extension of the previous problem, with an interior hole introduced in the torsional shaft. This is an example of a multiply-connected domain. The initial radii of outer and inner circle boundaries are 8 and 3 units, respectively. The optimization problem is simplified to define fourteen design variables in the first 60 Design variables X 1 x 2 x 3 x 4 x 5 x 6 6 6.26 4.97 4.97 6.36 6 Area = 85.26 Bending inertia on horizontal axis = 682.30 Torsional rigidity = 567.14 Figure 4.3 Initial design and design variable definition for torsional shaft 61 Design variables X 1 x 2 x 3 x 4 x 5 x 6 5.45 4.88 4.71 3.36 3.24 3.26 Area = 50.03 Bending inertia on horizontal axis = 350.81 Torsional rigidity = 306.01 Figure 4.4 Final design for torsional shaft with bending inertia constraint. 62 Design variables X 1 x 2 x 3 x 4 x 5 4.43 4.21 4.04 4.04 4.21 Area = 50.10 Bending inertia on horizontal axis Torsional rigidity = 325.54 x 6 4.23 = 206.98 Figure 4.5 Final design for torsional shaft without bending inertia constraint. 63 Y Figure 4.6 Design variable definition for multiply- connected torsional shaft. 64 quadrant, as shown in figure 4.6. The torsional rigidity and bending inertia requirements were similar to the solid shaft with additional requirements to constrain the total and inner hole section areas to 40 and 20 units, respectively. A total of 48 nodes were used in the analysis, with converged results obtained in 10 iterations. The initial and final designs are shown in figure 4.7. 4 • 3 Optimum Shape Design in 2-D Elasticity Problems Two plane stress problems are attempted in this section. A torque arm shown in figure 4.8, was sized for minimum weight and a maximum allowable von Mises stress of 620 Mpa . The maximum stresses in this structure occur at the boundary. A total of 133 linear boundary elements were used in the analysis, with six design variables used to define the geometry. Three additional geometry constraints were imposed to prevent boundary intersections during redesign. The optimization problem was solved using an exterior penalty function approach, with the evolution of the design shown in figure 4.9. The second example in this class of plane stress problems deals with the minimum area sizing of a fillet to sustain a uniaxial load shown in figures 4.10 and 4.11. Two different sets of design variables were used. The first set of results (figure 4.10) were obtained for a single shape 65 TO 03 C Figure 4.7 Initial and final design for multiply-connected torsional shaft. 10,000 N 66 O O o in Figure 4.8 Design variable definition and loading 67 Figure 4.9 Iteration history for torque arm. 68 Figure 4.10 Initial and final design of fillet (parabolic interpolation) Figure 4.11 Initial and final design of fillet (elliptic interpolation) 69 variable which defines the location at which a parabola from point A would meet the horizontal line OB with zero slope. The second set of results shown in figure 4.11 was obtained by determining the location of points A and B on the vertical and horizontal axis that could be connected by an elliptical curve, and satisfy the von Mises stress requirements. This problem is also attempted by defining design variables as radial locations from a fixed reference point figure 4.1. An optimum design could not be obtained for this set because of severe oscillations in the domain boundary. This problem also reinforces the need to use an adaptive node generation capability in the BEM based optimum design strategy. 4 . 4 Discussions In this chapter, we have successfully introduced the boundary element method in the shape optimization problem. A preliminary introduction to the definition of interior cutouts in planar stress problems, which is based on stress distribution in the domain, is also presented. The use of basic boundary shapes is advocated in such problems. The applications of using the boundary element method for shape optimal design has been demonstrated in several numerical examples . The gradient information used in the optimization 70 process for each of the problems was calculated by a finite difference approach. The sensitivity analysis consumes most of the CPU time, especially when the finite difference method is used. Seeking a more efficient approach for sensitivity computations becomes an important task in the shape optimization problem. An alternative approach of utilizing the characteristic of the boundary integral equation to obtain the gradient information, is studied in Chapter 5. The grids used in the two-dimensional elasticity problems in this chapter, were densely and evenly generated along structural boundaries. Therefore, the order of the system of linear equations obtained from the boundary element method becomes very large. Obtaining comparable levels of accuracy by reducing the number of grid points using the boundary element approach is extremely desirable. This would yield substantial savings in the time required the solution of the system equations, a process that is repeated extensively in both analysis and sensitivity computations. Toward this end, a grid refinement and adaptation technique is highly recommended when using the boundary element method, and will be discussed in Chapter 6. CHAPTER 5 SENSITIVITY ANALYSIS FOR BOUNDARY ELEMENT APPLICATIONS 5 . 1 An Overview of Sensitivity Analysis Methods In most structure shape design processes, the design sensitivity analysis comprises an important component. Design sensitivity analysis, which involves the calculation of the derivatives of the objective function and constraints with respect to the design variables, provides essential information required to couple mathematical optimization with structural analysis procedures. Since the calculation of sensitivity information is computationally demanding, the choice of a suitable approach for sensitivity analysis becomes very important. There are several approaches available to perform shape design sensitivity by use of the finite element method, but there is only limited experience in the use of the boundary element method for this purpose. Sensitivity analysis for the optimum shape design problem involving the use of the finite element method for analysis, is generally classified into four groups: (1) the direct difference approach, (2) the implicit differentiation approach, (3) the finite dimensional approach, and (4) 71 72 the continuum derivatives approach. (1) Direct difference approach method The direct difference method is commonly used in member sizing optimization. Botkin [14] used this approach to approximate the shape design sensitivity. By analyzing the structure for both old design variables, x, and new design variables, x + Ax, with given perturbation Ax, the design gradient is then calculated by the finite difference. The procedure of this method is straightforward, but is costly when the number of design variables becomes large, particularly if the calculations of the objective function and the constraints are expensive. The accuracy of the derivatives is dependent upon the step size of perturbation. (2) Implicit differentiation method This approach is usually applied in the finite element method [54]. The finite element equilibrium equation is written as df (x) f ( X+AX) - f(x) (5.1) dx AX [K] {u} = {F} (5.2) where [K] is stiffness matrix, (u) is displacement vector, 73 and {F} is the external load vector. Differentiating eqn. (5.2) with respect to a design variable, x, yields the following expression: 3u 3F <3K [K] { } = { } - [ ] {u} (5.3) ax ax ax The derivative (au/ax) can then be used to obtain the sensitivity of the stress response. The use of such an equation circumvents the necessity of solving for displacements for perturbed and unperturbed values of the design variables. The calculation of [aK/ax] , however, is expensive, especially if changes in the boundary may result in large changes in the stiffness matrix. Furthermore, the changes in shape may result in a distortion of the finite element grid, and the accuracy of response solution and hence also the sensitivity analysis, may deteriorate significantly. (3) Finite dimensional approach If the coordinates of nodes are considered to be the design variables of the shape optimization problem, then the finite element stiffness matrix is clearly dependent on the design variables. One can evaluate the derivatives of the stiffness matrix and the load vector with respect to the shape design variables, and the derivatives of the constraints can then be obtained by the adjoint variable 74 technique [8,9,55]. The number of adjoint equations is the same as the number of constraints and is independent of the number of design variables. Hence, if the number of constraints is less than the number design variables, this approach has an advantage over the finite difference technique . (4) Continuum derivatives Another approach for determining shape design sensitivity for the finite element method stems from the use of the concept of the material derivatives in conjunction with the adjoint method. The approach proceeds by forming material derivatives of the continuum equation [56,57]. Using integration by parts and prescribed boundary conditions, the derivatives can be expressed in terms of boundary integrals, which is a less expensive process than a direct computation of the derivative dK/dx. However, when the finite element method is used for analysis, there are many numerical difficulties associated with the evaluation of boundary integrals, especially for low-order curved boundaries . To overcome these numerical difficulties, a domain method, instead of a boundary method, is developed for shape design sensitivity analysis [53,58]. The domain integrals avoid numerical difficulties associated with boundary integrals but they are as expensive to use as the implicit 75 differentiation method of equation (5.3). The above methods are all well developed for the finite element method. With the exception of the direct difference approach method and the implicit differentiation method, the other three methods are not applicable for the boundary element method. In the BEM, the matrices [H] and [Q] in equation (3.20) are fully populated and nonsymmetric , and efficient solution methods that are applicable to banded and symmetric stiffness matrices, are no longer applicable. More recently, Choi and Kwak [32] have derived a sensitivity expression involvinq the direct boundary integral equation formulation in conjunction with the material derivative concept. In the following section, an alternative semi- analytical approach is discussed for the sensitivity analysis. This method, which utilizes the characteristics of the boundary integral equation and an advantageous use of the LU-decomposition technique, provides an efficient approach for the task at hand. 5 . 2 Semi-analytical approach in BEM An important ingredient in shape optimization is the computation of sensitivity of the design to the shape variables. The structural shape is governed by a group of independent design variables, where such variables affect 76 both the objective function and the structural response from which the design constraints are computed. A careful examination of the eqn.(3.20) indicates that the integral functions H and U, which are functions of r= | £ - x|, are basically determined by the geometry of the structure. When one of the shape design variables is perturbed by a small amount 6x, only a portion of the boundary related to that design variable will be changed. It is therefore logical to expect that not all elements in the [H] and [G] matrices are changed during a perturbation of one shape variable. This concept was used in this study to improve the computational efficiency of the sensitivity analysis . A better understanding of this idea is obtained by considering the boundary of a structural domain r, discretized into several elements as shown in Figure 5.1. If one of design variables d. is disturbed to d.+Sd., the 1 l i' boundary element associated with the design variable d^ will be changed to r ' . From the integral eqn.(3.16), one may note that if the load point x is chosen on the unchanged boundary, the difference in matrix elements of [G] and [H] , and [G'] and [H'], where the latter represent the system matrices of the perturbed domain, depends on the location of field point £ . That is, if the field point £ is also located on the unchanged boundary, these matrix elements 77 original boundary Figure 5.1 Original and perturbed boundaries discretized into N segments. 78 will not be changed. Otherwise, the matrix elements related to the field point £ which is located on the changed boundary, must be recalculated. Of course, if the load point x lies on the disturbed boundary, the matrix elements related to this load in [H'] and [G'J will also change. For simplicity, one may consider a linear interpolation function in eguations (3.18) and (3.19). The boundary of the structural domain is discretized into N linear segments as shown in Figure 5.1. The system eguations of the boundary element method in equation (3.21) can be written as follows . — ' ' A A • . • • A Cl ii A 12 a in A 21 A 22 '• ■ A 2N (X) = • C 2 . 1 > S3 H A N2 " ' a nn_ C N .. j If one of the design variables, dj_, is perturbed to d-L+Sd^, boundary elements related to that design variable will be affected. If for example, the design variable affects boundary elements between those denoted by indices K and L-l (see Figure 5.1), new matrices [A'] and (C'} are obtained as follows. 79 r A 7 A 11 A 7 A 12 A ' A IK • A 7 A 1L • • A 7 A IN C'l A ' 21 A > A 22 A ' A 2K . A 7 2L . . A 7 A 2N C'2 A 7 A K1 A ' .... A K2 A ' A K2 • A ' A KL A ' A KN c' K • (x'l = < : A ' A LI A' .... A L2 A ' A LK . A 7 LL A 7 LN c'l A 7 A N1 A N2 A ' A NK A NL A 7 A NN C'n (5.5) In the equation (5.5), only those elements within the cross- band are changed for the i-th perturbed design variable. By differentiating equation (3.21) with respect to the design variable d^, the following relation is obtained, 3x 3C dA [A] { } = { } - [ ] {x} (5.6) ddi ddi 3di where the matrix [A] on the left hand side is calculated before the design variable is perturbed. A finite difference approximation is used to obtain the matrix [dA/ddjJ and the vector (ac/adi) as follows, 3A [A']-[A] Ad, [ ad ] (5.7) 80 dC {C'}-{C> { } = (5.8) 3d Adi where [A'] and (C'} are as obtained in eqn. (5.5), requiring only those elements of [A] and { C } to be recalculated that are related to the boundary perturbations. The use of these finite difference approximations is the reason for the approach to be considered as a semi-analytical approach. In the present work, the system equation is first subjected to an LU-decomposition as follows. [A] = [L] [U] (5.9) The upper and lower triangular matrices, [U] and [L] , respectively, are only calculated once and are saved for all subsequent design variable perturbations. After the sensitivities [3A/3djJ and {SC/Sd^} are obtained, the solution (3x/3di) of eqn. (5.6) is solved in two stages. First, a forward substitution for the lower triangular matrix of [A] results in a solution for (H'} as, [L] (H'} = (H) (5.10) where , (H) = (3C/3di) - [ 3 A/ 3 d jj • (x) This is followed by a backward substitution for the upper 81 triangular matrix of [A] . [U] {ax/adi} = { H ' } (5.11) In the approach presented above, solution of the system equation for each perturbed design variable to obtain the gradient information, is replaced by the backward and forward substitutions. Since the matrix [A] is fully populated, the efficiency of using the semi-analytical approach for sensitivity analysis is evaluated by an operation count. For a complete set of gradient evaluations, the comparison of operation counts for the finite difference and the semi- analytical approach is shown in Table 5.1. In comparing with a direct application of the finite difference method, the semi-analytical approach produces computational savings of more than 50%. This is even more significant in the presence of a large number of design variables. Additional savings in computational cost can be realized by selectively computing those elements in the system matrix that are related to perturbations in the boundary. 82 Table 5.1 Comparison of operation count for finite difference and semi-analytical approaches to compute design sensitivity. finite difference approach 3 N 2 — ~ (NDV + 1) + N (NDV + 1) +0 (N) semi-analytical approach 3 N 2 — - + N (NDV + 1) + 0 (N) * [A] is an N*N matrix. * NDV is the number of design variable. CHAPTER 6 GRID REFINEMENT AND GRID ADAPTATION IN THE BOUNDARY ELEMENT METHOD 6 . 1 Introduction The boundary element approach offers several advantages over the finite element method in the shape design problem. There is, however, some judgement required in obtaining an appropriate discretization of the boundary. Since the geometry of the structural domain changes continuously during the shape design process, an improper grid distribution would yield inaccuracies in stress evaluation, particularly when the geometry is complex. In the following sections, both automated grid refinement and grid adaptation methods are considered to interface with the optimum shape design problem. The concept of grid refinement technique is first applied to refine the grid points along the boundary and to determine the total number of grid points for the structural analysis. After grid refinement, the grid points are redistributed along the structural boundary by an adaptive scheme. Two 83 84 adaptive schemes are studied in this chapter. The first is a variational approach which provides a simple, explicit means of incorporating desired characteristics into the grid generation scheme. The second approach utilizes existing optimization algorithms to redistribute the grid points on the structural boundary. 6 . 2 Grid Refinement Technique For a two-dimensional structural domain, the grid distribution on the boundary is usually treated as a one- dimensional problem. These grid points not only determine the nodes at which analysis is performed, but also define the structural boundary. Hence, proper caution must be exercised so that the structural geometry is not altered when redistributing the nodes for more accurate response analysis. To achieve this requirement, a set of master nodes are introduced along the structural boundary. The master nodes are not allowed to move during the grid relocation so as to prevent distortion of the original domain. The concept of grid refinement is similar to the h- method of FEM mesh refinement [59]. The grid refinement considered here is basically a local refinement procedure. The structural boundary is assumed to consist of several 85 zones which are defined by the master nodes. On the basis of a solution computed from the uniform grid, the grid distribution can be refined locally in each zone. In the proposed approach for grid refinement, a control function is first defined and computed for each zone, and the gradients of this control function determine the number of points that must be removed or added to that zone. Once the distribution of grid nodes is defined, the total number of analytical nodes is then determined. 6 . 3 Adaptive Scheme Based on Variational Principles In order to obtain more precise analysis results in the boundary element approach, grid points need to be clustered in regions of large response gradients [60]. In this section, a variational approach is used to obtain an adequate grid for response analysis. The desired features of such a grid are represented in terms of a control function defined along the structural boundary, where larger values of the control function gradients require a proportionately larger number of computational grid points to be assigned to the corresponding region. For a problem pertaining to a one-dimensional adaptation of grid points along the boundary, the boundary coordinate s which defines the length along the boundary, and a generalized curvilinear 86 coordinate £ which pertains to the computational grid, are as shown in Figure 6.1. One form of a functional Ip, which can be used as a measure of adaptation of the grid spacing in the f direction, can be written in terms of a control function P as [61], t 2 s S ds (6.1) P where £ s = <3£/ ds, and is referred to as the grid density (number of grid points per unit length) . An example of choice of the control function will be given later. For the one-dimensional problem, the relation between £ s and s^ is expressed as follows. It is observed from eqn. (6.1) that, for the functional I p to be minimum, a region on the boundary with a small value of control function P would have a corresponding low value of £ s . Stated differently, a larger grid spacing is obtained in the region of small P. For a boundary with a fixed number of points, the above implies a denser grid distribution in regions of high P. The variational approach yields the necessary conditions for optimality that a solution for the grid 87 Figure 6.1 Definition of boundary and generalized curvilinear coordinates. 88 density must satisfy to extremize the functional Ip. These are the familiar Euler-Lagrange equations, and for the present problem result in a second order, one-dimensional Poisson equation. ? ss - ^s * P s / P = 0 ( 6 ' 3 > This equation governs the location of the grid coordinates in the physical domain. The solution of this equation is facilitated by a transformation to the curvilinear coordinate system. Using a chain rule, the second derivative terms can be shown to be of the following form. e ss (6.4) Substituting equations (6.2) and (6.4) into eqn.(6.3) yields the transformed one-dimensional equation for adaptation along the £ coordinate. + * P £ 7 P = ° (6 ' 5) This equation can be solved directly for the grid point spacing along the £ -coordinate . Another approach of looking at the adaptation problem is obtained if one integrates eqn. (6.5) once to obtain, 89 where, P • S£ = C ( 6 . 6 ) C I and, represents the total arc length along the boundary. It is obvious from the above expression that the relationship between the control function P, and the grid point spacing s , is reciprocal. A proper selection of the control function P is of essence in an adaptive relocation of grid points along the structural boundary. In the problems considered in this exercise, the control function was defined so as to provide adequate clustering of grid points in regions of high stress and high curvature of the boundary. One such control function definition is of the following form, f = £ * a v ( 1- £) * c k (6.7) where f is a control function; a v is the von Mises stress along the boundary; c^ is the curvature of the boundary; p is a weighting factor that assumes values between zero and unity, and determines the relative importance of the stress or curvature in the adaptation process. The ratio of the maximum to the minimum spacing in the computational grid is 90 a measure that should be controllable by a proper selection of a control function. To avoid numerical problems associated with the control function assuming a very small value, a modified control function is defined as follows, P = 1 + 7. f (6.8) where the function f is normalized to range between 0 and 1, and 7 is a smoothing factor. The minimum and maximum values of P correspond to f =0 and f=l, respectively. It can be shown that the expression for grid spacing from i to i+1 is of the following form. 1 t 1 + 7 • f i As j_ = j 1 + 7 • f j ( 6 . 9 ) where is the total length of the boundary, and P^ is the value of the control function between i and i+1, and can be written as follows. Pi - 1 + 7 . f i ( 6 . 10 ) Hence the ratio of the maximum to the minimum spacing may be 91 expressed as follows. ( ASi ) ( ASi ) max min 1 + 7 ( 6 . 11 ) The value of this ratio is clearly dependent on the smoothing factor 7 which is prescribed in the control function. However, as explained in subsequent sections, the approach used in this work for grid adaptation and refinement is relatively insensitive to the choice of this parameter. 6 . 4 Grid Adaptation Obtained by Nonlinear Programming Methods The variational approach described in the previous section to obtain the computational grid is applicable in several engineering problems. However, if the functional Ip in equation ( 6 . 1 ) is complicated, or if additional constraints must be imposed in the grid adaptation, a solution for grid point spacing is not as easily available. The use of nonlinear programming methods to relocate the grid points is a more generally applicable strategy. This approach can be extended to include additional constraint conditions. Furthermore, in the event that several noncommensurate criteria are to be considered for the grid 92 point locations, a multicriterion optimization approach can be adopted. Instead of solving the partial differential equation for the optimal grid point locations, a nonlinear programming algorithm is utilized to minimize the functional Ip. If the boundary of the structure is discretized into N segments, the functional Ip can be written as follows, N 1 Ip = S (6.12) i=l Pi Asi where, as before, Pi is the control function in the i-th region and Asi represents the length between two adjoining nodes, (i-1) and (i). The mathematical statement of the nonlinear programming based optimization problem for grid adaptation can be written as follows. Find Si, S 2 , ... , sjj to minimize Ip (6.13) Here s^,s 2 , ... , s^ are the length coordinates of all moving grid points. It is important to note that this set does not include the master nodes. This statement is basically that of an unconstrained optimization problem, for which several efficient solution algorithms such as the Fletcher-Reeves method [62] and the Broydon-Fletcher-Goldf arb-Shanno variable metric method [63-66] are available. The gradient of the objective function is needed in such first order 93 methods, and can be easily obtained by forward or backward finite difference method. However, the size of the perturbation step in such a difference approach is very important. Since the master nodes are not allowed to move during the grid redistribution, an arbitrary perturbation of the design variables can cause a grid point to be located in an undesired zone. This is most likely to happen if the perturbation is a fixed fraction of the design variable, as is frequently the case in member sizing problems. In the present work, a fixed design variable perturbation was applied to all design variables. The problem of grid adaptation for more than one criterion function, and subsequence use of the nonlinear programming approach, is explained as follows. Consider a control function that can be separated into two distinct objective criteria, one based on the behavior response such as the von Mises stress, and a second based on the local curvature of the structural boundary. The location of grid points must be optimal from the standpoint of both criteria, resulting in a multicriterion optimization problem. In this study, an implementation of the global criterion approach [67] was used to obtain solutions to the grid adaptation problem. This approach is summarized here for completeness. Consider a situation in which the control function consists of two candidate criteria, fi(s) and f 2 (s), representing the importance of the von Mises stress function 94 (ctv) and the curvature of the boundary (ck) , in the problem of locating grid points. These functions can be written as follows . N 1 f 1 (s) = s i=l aV i ASi (6.14) N 1 f2(s) = S i=l ck-[ Asi (6.15) Each candidate criterion may be optimized individually, ignoring the existence of the second criterion. The resulting optimal values of the objective functions fi lc ^(s) i d and f 2 (s) are referred to as the ideal solutions for each criterion function. The global criterion method is based on the premise that the design vector for which both criteria are simultaneously optimized may be obtained by minimizing a metric measure of the distance between the ideal solutions and the true optimum. As shown in [67], this results in a min-max problem which can be converted into an equivalent scalar minimization problem as follows. Minimize p Subject to f ! ( s ) - fi ld (s) - p < 0 95 f2( s ) ~ f2 ld ( s ) w 2 Id f 2 (s) - p < 0 where , w 1 + w 2 = 1 Here «]_ and w 2 are the weighting coefficients representing the relative importance of each criterion. Figure 6.2 illustrates the feasible and infeasible spaces for a two- criterion function. By choosing different weighting factors and ^ 2 , the ideal solution f ld (s) is located on the boundary between f 1 ld (s) and f 2 ld (s) . 6 . 5 Examples of Grid Refinement and Adaptation The primary motivation for using an adaptive grid is to allocate an appropriate number of grid points to regions of high stress and curvature. In problems of shape design using the BEM, grid points play a dual role of defining the structural shape and providing discrete points for analysis. In computing the sensitivity of the structural response to changes in the boundary, where the latter is linked to the position of grid points, proper care must be exercised to eliminate distortions in the geometry of the boundary. The master nodes uniquely define the shape of the structure and serve as points for response analysis as well as provide the interpolation points for the structural geometry definition. 96 Figure 6.2 Graphical representation of a two-criterion function design space. 97 These nodes, shown in Figure 6.3(a), are not moved when the adaptive scheme is applied. Grid adaptation and refinement is essentially a two- stage process. A uniformly distributed grid is first assumed between the master nodes, and the corresponding approximate solution is determined. This solution allows a computation of the control function and its gradient along the boundary. The approach adopted is one in which more grid points are assigned to the high gradient regions, and less or no grid points are located in regions with a smooth variation of the control function. A relatively simplistic scheme was implemented to realize this goal. The numerical values of the gradients of the control function were classified into several levels, and a proportional number of grid points were assigned to each of these levels. The effect of such a redistribution of points based on control function gradients is illustrated in Figure 6.3(b). Once the total number of grid points for each region is determined, the governing equation for grid adaptation is employed to provide a better grid redistribution. This is shown in Figure 6.3(c). Note that in the solution of grid adaptation, the master nodes are always fixed to maintain the boundary of the structure. Figure 6.4 illustrates the process of grid point distribution on the basis of control function values. After an approximate solution is obtained for a uniform grid » HH» »M »» 98 • G C (3 (/) -HO •P ft) *fH G • (3 TS -P <w tfl -P O 43 g O G 0 -P 0 •H U 4-> 0 C Q) G G -P (0 r TJ g C -p 3 (0 O P 0) 0 O 4-1 ■P P CD H (D C'rl ft a a E <3 -P -P -P rH P W 44 o TS O m a ts <0 CD 3 1C (1) 3 -H ft (0 -p g p n p -p rH u o t?c TS (0 •r- 3 T) ft O -P C <4- P <u -H -p T3 o P O ■P 0 P G CD O -H CD W (0 O -H 0 43 •P Cft (0 CP • (3 — - W CD a-H .-H <3 43 43 W 43 U P ■7-J (0 n H n >,4-i . n > •1 — • f3 • tP rp O W VO <D vo G G D vo G Cr» > -H CD W G O (DO) > > G CD r~> P i — 1 0) O l 3) 43 O -H T3 TS 0) CP-P > 5 CP (0 W •P -P G CD <D p <0 P S -H S3 P CP X! ( c ) grid adapta ti on 99 distribution as shown in Figure 6.4(a), the grid is refined between each pair of neighboring master nodes on the basis of the gradient of the control function. In Figure 6.4(a), we observe that the gradients of the control function around points A, B and C are much higher than elsewhere. Thus, several new grid points are assigned to these regions and evenly distributed within the regions. A smaller number of grid points are assigned to regions of low gradients, as shown in Figure 6.4(b). The points placed in these regions are next relocated by a solution of the adaptive equation, and the results of this adaptation are shown in Figure 6.4(c). It is worthwhile to bear in mind that since master nodes have been introduced along the boundary, the maximum spacing (^ s i) max •*- n the eqn. ( 6 . 11 ) is equal to the space between master nodes. This also renders the adaptation scheme relatively insensitive to the choice of the smoothing factor 7 . Table 6.1 compares the von Mises stress distribution along a structural boundary for a problem solved with 149 evenly distributed analytical points, and with 61 analytical points assigned on the basis of grid refinement and adaptation. The differences in the numerical results are small, even though the number of grid points in the second case are less than half the number in the first. Such a reduction of grid points can result in substantial savings in computational cost in shape design problems. bSS y fUnCtl ° n dist «bution along the ■ 4<b> control f function^ aSe< ^ ° n gradients the Control Function 101 Figure 6.4(c) Grid adaptation in regions of high stress. 102 Table 6.1 Comparison of ct v /(7 ,, between 149 uniform grids and 61 redistributed grids. selected master nodes av/a al] case 1 : uniform grid with 149 analytical nodes case 2 : adaptive grid with 61 analytical nodes 1 0.8500 0.8629 2 0.1207 0.1197 3 0.0999 0.1029 4 0.1144 0.1121 5 0.1234 0.1185 6 0.1144 0.1121 20 0.5154 0.5092 21 0.0596 0.0595 22 0.0493 0.0538 23 0.0650 0.0480 24 0.0411 0.0404 25 0.0259 0.0254 26 0.0019 0.0096 27 0.0470 0.0484 28 0.1166 0.1144 29 0.2640 0.2690 * a v is von Mises stress. * CT a n is allowable stress (62,000 N/cm 2 ) . 103 Another example illustrating grid adaptation, based on a nonlinear programming solution of a multicriterion optimization problem, is presented next. Hypothetical variations for av and ck, defined in preceding sections, were assumed. These can be written as explicit functions of the coordinate variables as follows. av(s)=sin(s)-0.6*sin(s)+0.1*sin(s) ck(s)= |sin(3*s) | Figure 6.5 shows a distribution of nineteen nodes in the interval between 0 and n , with the two end points fixed. Figures 6.5(a) and 6.5(b) show the adapted grid distributions based on av(s) and ck(s) , respectively. The values of the ideal solutions f^ ld (s) and f 2 ld (s) were 113.44 and 109.09, respectively. Figure 6.5(c) shows the grid distribution when crv(s) and ck(s) were considered egually important in the adaptation process. This distribution was obtained as a solution of a multicriterion problem where the weighting factor u>± and o >2 were both set egual to 0.5. Figure 6.5(d) shows the grid distribution for a case where curvature is weighted more than the stress distribution. This distribution corresponds to values of ui=0. 333 and ^2=0. 667. 104 (a) (b) (c) (d) 4 H Milt H-B ■■II 1 1 i 1 1 l"l- i — li 1 <-4- h 1 H 1- I i rri 1 — r-i i i i— i ■ t i ' ■ i i i i i ■ ■ 1 ■ 1 1 1 1 ■ — *- ^ 1 — -i H4-4- — r i j i i i i i i i -H — H-H-H— b- 1 i 1 i -HH— -M- Figure 6.5 Grid distribution obtained by using a control function of the form + u 2 -ck. case case case case (a) (b) (c) (d) w i = w 2 = o, “1 = 0, w 2 = 1, w]_ = 0.5, w 2 = 0.5, — 0.333, w 2 =0.667, 105 6 . 6 Optimization examples using grid refinement and grid adaptation techniques In this section, two optimization examples are used to illustrate the function and the efficiency of applying the grid refinement and grid adaptation schemes in conjunction with BEM for shape optimal design. The framework of the optimization procedure is shown in Figure 6 . 6 . The semi- analytical approach is employed to obtain the gradient information . A flat plate, supported as shown in Figure 6.7, is subjected to a concentrated load of 45,000 N. The allowable stress, Poisson ratio, and Young's modulus for the plate material are 620 Mpa, 0.3, and 226 Gpa, respectively. The dimensions of the plate and a definition of the shape design variables is also shown in this figure. The distance between the support is fixed, as are the points A and B. Other boundaries are defined in terms of half-parabolic curves. Design variables x^, X 2 and X 3 , X 4 control the right outer and the left outer boundaries, respectively. Design variables X 5 and xg define the height of the two parabolic curves of the inner boundary. Symmetry of the loading and support dictate that x^=X 3 and X 2 =X 4 . Figure 6.8 illustrates the design iteration history for this plate. Starting from the initial design of the structure shown in Figure 6 . 8 (a), the design after four 106 System equations are solved by LU-decomposition . Upper and lower triangular matrices are saved for sensitivity evaluation. Compute objective function and constraint information. Use the semi-analytical approach to evaluate the gradient of the objective function and constraints. Forward and backward substitution is used to solve the system equations. Figure 6.6 A fully nonlinear optimization algorithm with grid ridistribution and semi-analytical computation of design sensitivity. 107 45 KN Figure 6.7 Definition of geometry and design variables for the elastic plate problem. 108 cycles is as shown in Figure 6.8(b). At this point, the original node distribution is no longer suitable for analysis, and the application of grid refinement and grid adaptation results in a node distribution shown in Figure 6.8(c). The optimum shape of the structure is obtained in the sixth cycle, and is shown in Figure 6.8(d). The second test problem involved the minimum weight design of a fillet to sustain a uniaxial distributed load of 40,000 N/cm, as indicated in Figure 6.9. The material properties are the same as in the previous example. The structural boundary to be redesigned is included between points A and B, as shown in the Figure 6.9. This boundary was discretized into ten equal segments, and all nodes on this line (excluding points A and B) were allowed to move perpendicular to the line AB during the redesign process. This results in a total of nine design variables for the problem. This is a classical problem in optimal shape design, and is known to present problems in obtaining a converged optimum. Changes in boundary curvature, particularly the introduction of sharp corners on the boundary, results in significant movements of the zones of high stress. The numerical optimization exploits inherent weaknesses in discrete analysis methods related to the fact that analysis is only performed at predetermined points. The design process, therefore, exhibits an oscillatory behavior. There 109 ■p (0 43 -p QJ • c p •P g 0 3 3 <U *H 0 iH rH -p <p ft 43 3 0 13 >p 0 P "P 0 •P ft P ■P •P Cn w <u w c 3 •P -P •P ip (0 T3 c ( 1 ) P • •P ftT3 a) Cn a) •P i — 1 Q) 13 0 P 0 43 ■p •h trN -P 0 -P P [/) iH <v 0 3 (0 13 <w rH -p -p G <U -P P 0 c •H 3 •H tr> P G 0 -P •p 0 -P <W 3 V) <p 43 a) C <D •H T3 c o £ p cn -p .p 1 — I -H T} W 3 w <d <p •H E (D M Ofl ■P •a 3 -P 13 CPTS ft • '-I G • ■H 0 E 3 G -P U • <u •H Cr> G Cr> (U r-H i— i ■P -H -P rH 3 13 •h tn tr> 5 u G 0 G <D < 1 ) 0 ) Jj HQ13! 2 O ft ft <0 43 O T? CO CO CO CO VO VO KO • VO <D P 3 •P ft 110 E u Figure 6.9 Definition of geometry and loading for the Ill are essentially two approaches to counter oscillations in such problem. The first espouses the use of a very large number of points on the boundary to obtain better resolution of regions of high stress. This has the associated drawback of larger computational requirements. The other approach is one which uses grid refinement and adaptation techniques to locate a fewer number of points in the critical regions. The order of the resulting system equations is low, yielding superior computational efficiency. The initial design for the fillet problem is as shown in Figure 6.10(a). An adaptive grid was defined and the design allowed to proceed with prescribed move limits on each design variable during the optimization. The presence of such move limits alleviates the computational burden to some extent, as a new grid refinement is not necessary after each iteration. Figures 6.10(b) and 6.10(c) illustrate the progression of the design to an optimum. 112 B " * — I— * — * 1* * * w Initial Design m » T- | A ,. r * A 4 ► \ ■ 1 1 4 1 * 1 • 1 * Cycle 3 w i —4 4 4 > 1 * 1 • \ ■ \ • 1 ■ • • * * ,, Final Design - Cycle 5 Figure 6.10(a) 6.10(b) 6.10(c) Initial design for the fillet problem. Design at begining of third cycle. Final optimal design for the fillet problem obtained at the fifth cycle. CHAPTER 7 AN INTEGRATED APPROACH FOR SHAPE OPTIMIZATION USING THE BOUNDARY ELEMENT METHOD 7 . 1 Introduction The use of an integrated formulation in structural optimization was first proposed by Schmit and Fox [5,6]. The method proposed collapsing all constraints, including the equality constraints (a set of linear equilibrium equations) and inequality constraints (behavior constraints) , into an unconstrained function of several design variables, and obtaining a solution to the unconstrained problem by the penalty function method. This approach was abandoned at the time for lack of efficient methods to solve the poorly conditioned unconstrained optimization problem. In view of the recent developments in nonlinear- programming-based optimization methodology, the approach of simultaneous analysis and design has again surfaced as a viable option. Haftka [7], uses an element-by-element preconditioned conjugate gradient method in the solution of 113 114 a linear analysis problem, and Newton's method in a design subject to a nonlinear collapse load constraint. When optimization techniques are integrated with structural analysis, the design problem can be viewed as a nested optimization problem. Since the computation of behavior variables such as stresses and displacements is unnecessary for each trial design vector x, the use of this approach reduces the solution of the analysis equations to an iterative process, which proceeds in conjunction with the optimization solution. This approach not only eliminates the evaluation of a large number of intermediate trial designs, but also achieves optimum results of both the response solution and shape design simultaneously. Most integrated approach formulations available in the literature pertain to the structural member sizing problems, and are limited to finite element based analysis. In this study, the integrated approach is extended to shape optimization problems, wherein the BEM analysis equations are represented as a set of equality constraints in the optimization problem. This approach precludes an explicit solution of the set of linear algebraic equations (BEM equations) , at the cost of an increased number of constraints and design variables (response variables are now design variables) in the optimization problem. The generalized reduced gradient (GRG) method is employed to 115 obtain a solution to the combined analysis and design problem. 7 • 2 Integrated Approach Using The Boundary Element Method In the conventional structural optimization approach, the optimization and analysis techniques are independent tools in the design process. The optimization approach provides a mathematical tool to effectively search the design space so as to achieve an optimal objective function value, and simultaneously satisfy a set of constraints. The analysis simply provides the optimizer with necessary information to proceed towards the optimum. In this study, an integrated formulation is adopted, coupling the boundary element analysis equations with the nonlinear programming based optimization problem. The BEM formulation for analysis of a structure typically leads to a system of linear equations, [A] {v} = {C} (7.1) where the response quantities v in equation (7.1), are generally solved by Gaussian elimination. In the integrated approach, these linear equations are treated as additional equality constraints in the optimization problem, combining the components of the unknown response vector v with the 116 vector of shape design variables d. Further, the equations relating the stress response and the response vector v are also included as additional equality constraints. Such equality constraints relate the von Mises stress to the response variables v and are represented as, a i ~ CTV i(v)= 0 (7.2) where cr-^ is the stress design variable which represents the value of von Mises stress at boundary node i, and avi(v) represents the von Mises stress as a function of the response variables v at node i. Additional inequality constraints are obtained due to requirement of allowable limits on the stress response. CT i "iall s 0 (7.3) where CT i al ^ is the i-th allowable stress prescribed at the beginning of the design cycle. A general mathematical statement for the integrated optimization problem involving the boundary element method may be written as follows. 117 Find d, v, and to minimize W Subject to Av - C = 0 a j_ - crvi(v) = 0 ° i ai all (7.4) where d is a vector of shape design variables, and W is the structural weight. For a given set of design variables, the structural weight is only a function of shape design variables d. The matrices A and C in the equality constraints, and the boundary stresses cri(v) in the inequality constraints are obtained from the boundary element formulation. As was stated in Chapter 4, numerical accuracy of the stress response is extremely sensitive to the grid distribution. Improper node point locations would result in misleading information from the BEM based analysis. Although, a uniformly dense grid distribution for the entire boundary may yield better accuracy, the process can become computationally prohibitive. In using an integrated optimization strategy for a two-dimensional elasticity problem, the addition of one node results in the introduction of three equality constraints (two for boundary element equations and one for stress equality equation) and 118 one additional inequality constraint. Furthermore, three more design variables (two v and one ) are added into the design space. From the standpoint of computational efficiency and reduction of memory requirements, the total number of nodes must be minimized. The proper node distribution required in the solution of the problem is determined by the grid refinement and adaptation techniques discussed in the Chapter 6. The shape design variable vector d changes continuously during the optimization process. After just a few iterations, the initial grid distribution may become unsuitable for the new boundary shape. To avoid a new grid definition for every iteration, move limits were imposed on the shape design variables. A prescribed grid distribution can therefore remain effective for an admissible and controlled variation of the boundary. The adaptive grid distribution method was used in conjunction with an integrated analysis/design procedure to determine the minimum weight geometry of a fillet to sustain a uniaxial distributed load of 35,000 N/cm. This minimization was subject to constraints of maximum allowable value of the von Mises equivalent stress on the boundary. The dimensions and material properties of the structure are indicated in Figure 7.1. As shown in this figure, the structural boundary to be redesigned is included between points A and B. This boundary was discretized into eight 119 h* 10 cm i — 5 cm *21 18 Figure 7.1 A planar fillet design problem. 120 equal segments, and all nodes on this line (excluding points A and B) were allowed to move perpendicular to the line AB during the redesign process. This results in a total of seven design variables for the problem. After grid refinement and adaptive procedures, total number of grid points and the grid locations are determined for structural analysis. In this problem a total of 37 nodes were retained for structural analysis. These nodes were distributed as shown in figure 7.2. The number of design variables for the integrated approach is 121, which consists of 7 shape design variables d, 74 response design variables v, and 37 stress design variables a The values of components of the vector v for this problem typically ranged from 10 3 to 10 7 . An improper choice of initial values of v resulted in inordinately high computational expense, and may, in some cases, preclude the determination of an optimum design. A rational choice of initial design variable values is therefore essential in the integrated approach. For a given set of initial shape variables d, exact value estimates of the response design variable vector v and aj can be obtained by an evaluation the BEM based analysis equations. This process can be periodically repeated in the iterative design process. The total number of constraints is dependent on the number of grid points. In the present problem, there were 74 equality constraints corresponding to the linear equations 121 Figure 7.2 37 adaptatively redistributed node points for fillet problem. 122 (7.1), 37 equality constraints for the stress response computations (7.2), and 37 inequality constraints for stress allowable limits (7.3), respectively. Since the boundary element equations are linear with respect to the response variable vector v, the matrix A can be directly used as the gradient of the equality constraints corresponding to the BEM analysis equations. The gradient information for all other constraints may be obtained by the finite difference approach. The integrated optimization problem was solved by the generalized reduced gradient (GRG) method, with prescribed move limits established for each shape design variable. The changes of shape design variables for each piecewise linear cycle were restricted to within 20% of starting values. Within these move limits, first order linear approximations of the constraints were used. The corresponding results were obtained in five cycles and the final shape is shown in Figure 7.3. Such an integrated approach eliminates the necessity for repeated decompositions of the boundary element equations, as would be necessary in the traditional formulation. However, it is worthwhile to note that the integrated approach results in an increase in design variables and constraints, giving rise to an increase in computational costs. This increase is largely due to the fact that tight move limits must be imposed in the optimization problem because of the poor numerical scaling 123 Figure 7.3 Final shape design for fillet problem. 124 of the design space. For instance, the value of equality constraints (7.1) must necessarily be very small to obtain accurate values of the response variables. Multiplying factors of order 10 2 to 10 3 , were used to scale the equality constraints to improve the accuracy of the solution of the response design variables, stress response quantities (a^) , and the shape design variables (d) , in the integrated approach. CHAPTER VIII CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE RESEARCH Some key ideas related to an efficient optimal shape synthesis of elastic structures have been discussed in this study. The boundary element method is used for analysis as it provides better accuracy in response computations at the boundary. Two-dimensional elastic, torsional shaft problems, and small-strain elastostatics problems are successfully solved by the BEM. Both constant and linear elements are used for the boundary element analysis necessary for shape synthesis. A strategy for choosing proper shape design variables is examined. The method is based on the computed load path in the structural domain. Gradient-based, nonlinear programming algorithms are investigated in the shape optimization problems. In general, the direct NLP methods require fewer computations of objective function and gradients than the indirect NLP method. For the NLP problem involving several equality constraints, such as in the integrated approach, the GRG method is effectively employed to obtain optimum solutions. 125 126 The evaluation of gradients in the NLP method is done by semi-analytical methods, where such methods are shown to be significantly more efficient than finite difference approximations. In the present work, this approach was implemented in computation of design sensitivity for a torsional shaft problem and for a class of two dimensional, linear elastostatics problem. The concept of this gradient evaluation can be readily extended to high order elements of BEM for shape optimization problems. The research also describes an effective approach for grid refinement and adaptation in shape synthesis problems. The minimization of the total number of grid points and their proper location is especially significant in the BEM formulation. This is true from the standpoint of both numerical accuracy and computational efficiency. The total number of analytical nodes is determined by grid refinement. Two adaptive schemes, one based on a variational approach, and the other on a nonlinear programming-based unconstrained optimization algorithm, are examined in this work. It is worthwhile to note that the nonlinear programming method allows for addition of constraints as well as an ability to account for multiple criteria in the grid adaptation. In using such grid refinement and adaptive schemes, the order of system equations in the BEM analysis can be reduced significantly. To avoid a new grid definition for every iteration, move limits were imposed on the shape design 127 variables. A prescribed grid distribution can therefore remain effective for an admissible and controlled variation of the boundary. The research also examines the feasibility of using an integrated formulation for a BEM based optimum shape synthesis. The efficiency of the proposed scheme is attributed to eliminating the necessity of decomposing the system matrix in the BEM formulation. This matrix is not banded as in the finite element formulation, and is therefore not amenable to efficient decomposition methods. The inclusion of behavior variables as additional design variables in the integrated optimization formulation often results in a poorly scaled system, and detracts from the efficiency of the process. With regard to future work in this direction, a substructure shape optimization concept which combines both the FEM and the BEM is highly recommended. In many practical design problems encountered in automotive or aircraft structures, the structural domain is very large, and can not be properly discretized by the boundary element method. The FEM can be first applied to analyze the entire structural domain. For a specific portion of the structure where additional benefits may be derived from a shape redesign, a substructure can be created with the pertinent boundary load information obtained from the FEM. The shape optimization problem for redesigning the shape of this 128 substructure may then use BEM based optimization techniques discussed in this study. The current research primarily focused on two- dimensional static structural shape optimization involving the use of BEM. More recently, advances of the BEM in problems related to wave propagation, structural vibrations, and three-dimensional structural analysis have been reported in [3,44]. The methodology developed in this study for efficient BEM-based shape synthesis can be extended to these disciplines with the attendant advantages discussed in earlier sections. 129 APPENDIX In most pratical applications the designer is interested not only in stresses inside the body and tractions on the surface, but also in the complete stress tensor on the boundary. The calculation of this stress tensor proceeds by expressing the tractions in a local coordinate system, and employing the strain-displacement along a direction, tangent to the boundary. Assume a two-dimensional boundary element problem with a local Cartesian coordinate system with axes tangential and normal to the boundary (Figure A.l). The relation between the global coordinate x and local coordinate x is as follows. The above relationship can be expressed in terms of a transformation matrix [T] as follows. x = [T] x For a plane strain problem in isotropic linear elasticity, the stress tensor on the boundary can be 130 Fingure A.l General two-dimensional boundary element and local coordinate system defined at the stress point. 131 expressed in term of tractions and the interior unit extension tangential to the boundary, as follows. -11 - Pi -12 = P2 *22 = — — (" -ii + 2G " 2 2 ) Hence, the stress tensor may be obtained as follows. [-] = [T] T [a] [T] Using the two-dimensional plane strain relations and the expression, [P] = [T] [p] the stress tensor along the boundary can be expressed follows . a 11 2 —2/ = (sin 3 a + COS 2 a 1 - 1 / + + (sin 2 a 2G l-i/ COSa - l-i/ COS 2 a T 22 sina) pi COS 3 a ) P2 (A. 2) (A. 3) (A. 2) (A. 4) as (A. 5) 132 v (7 = ( COSa sin 2 a - COS 3 a ) Pi X Z 1 - 1 / + (sin 3 a - COS 2 a sina) P 2 1 - 1 / 2G + COSa sina ^"22 1 - 1 / a = ( sin 3 a - sina COS 2 a) Pi 1 — 1 / 2 -1/ + (-cos 3 a - COSa sin 2 a ) P2 1 - 1 / 2G 22 (A. 6) + l-i/ sin 2 a e (A. 7) 133 REFERENCES 1. R.T. Haftka and R.V. Grandhi, "Structural Shape Optimization - A Survey", Computer Method in Applied Mechanics and Engineering, 57, pp. 96-106, 1986. 2. N. Kikuchi, "Adaptive Grid Design for Finite Element Analysis in Optimization: Part 3, Shape Optimization," ed. by C.A. Brabbia Martinus Nijhoff publisher, 1983. 3. P.K. Banerjee and R. Butterfield, "Boundary Element Methods in Engineering Science," McGraw-Hill, 1981. 4. G. Kuich, "Starting to Work with Boundary Element Technique in Computer Aided Engineering", ed. by C.A. Brebbia, Martinus Nijhoff publisher, 1983. 5. R.L. Fox and L.A. Schmit, "Advances in the Integrated Approach to Structural Synthesis," AIAA, Vol. 3, No. 6, pp. 858-867, June 1966. 6. L.A. Schmit and R. A. Fox, "An Integrated Approach to Structural Synthesis and Analysis", AIAA, Vol. 3, No. 6, pp. 1104-1112., June 1965. 7. R.T. Haftka, "Simultaneous Analysis and Design," AIAA, Vol. 23, No. 7, pp. 1099-1103. 8. O.C. Zienkiewicz and J.S. Campbell, "Shape Optimization and Sequential Linear Programming," Optimun Structural Design, ed. by R.H. Gallagher and O.C. Zinkiewicz, John Wiley & Sons, New York, pp. 109-126, 1973. 9. C.V. Ramakrishnan and A. Francavilla, "Structure Shape Optimization Using Penality Functions," J. Struct. Mech. , 3(4), pp. 403-422, 1974-1975. 10. A. Francavilla, C.V. Ramakrishan and O.C. Zienkiewicz, "Optimization of Shape to Minimize Stress Concentration," Journal of Strain Analysis, Vol. 10, pp. 63-70, 1975. 11. E.S. Kristensen and N.F. Madsen, "On the Optimum Shape of Fillets in Plates Subjected to Mltiple In-plane Loading Cases," Int. J. for Num. Meth. in Engrg. , Vol. 10, pp.1007- 1019, 1976. 12. Y.W. Chung and E.J. Haug, "Two-Dimensional Shape Optimal Design," Int. J. for Num. Meth. in Engrg., Vol. 13, pp. 311- 336, 1978. 134 13. K. Dems and Z. Mroz, "Multiparameter Structural Shape Optimization by The Finite Element Method," Int. J. for Num. Meth. in Engrg., Vol. 13, pp. 247-263, 1978. 14. M.E. Botkin, "Shape Optimization of Plate and Shell Structures," AIAA Journal, Vol. 20, pp. 268-273, 1982. 15. J.P. Queau and P. Trompette, "Two-dimensional Shape Optimal Design by The Finite Element Method," Int. J. Numer. Eng., No. 15, pp. 1603-1612, 1980. 16. M.H. Imam, "Three-dimensional Shape Optimization," International Journal for Numerical Method in Engineering, Vol. 18, pp. 661-673, 1982. 17. J . A . Bennett and M.E. Botkin, "Structural Shape Optimization with Geometric Description and Adaptive Mesh Refinement," AIAA Journal, Vol. 23, pp. 458-464, 1985. 18. M.E. Botkin and J.A. Bennett, "Shape Optimization of Three-Dimensional Folded Plate Structures," AIAA Journal, Vol. 23, (11), pp. 1804-1810, 1985. 19. N. Kikuchi, K.Y. Chung, T. Torigaki, and J.E. Taylor. "Adaptive Finite Element Methods for Shape Optimization of Linearly Elastic Structures," edited by J.A. Bennett and M.E. Botkin, Plenum Press, New York, 1986, pp. 139-166. 20. H.C. Rodrigues, "Shape Optimal Design of Elastic Body Using A Mixed Variational Formulation," Computer Method in Applied Mechanics and Engineering, 69, 1988, pp. 29-44. 21. K.K. Choi and H.G. Seong, "A Domain Method for Shape Design Sensitivity Analysis of Built-Up Structures," Computer Methods in Applied Mechanics and Engineering, 57, 1986, pp. 1-15. 22. N.V. Banichuk and B.L. Karihaloo, "Minimum Weight Design of Multi-Purpose Cylindrical Bar," International Journal of Solid and Structures, Vol. 12, No. 4, April 1976, pp. 267- 273 . 23. K. Dems, "Multiparameters Shape Optimization of Elastic Bars in Torsion," International Journal for Numerical Methods in Engineering, Vol. 47, No. 10, Oct. 1980, pp. 1517-1539. 24. S. Adali, "Cross-Sectional Shape of An Anisotropic Hollow Bar for Maximum Torsional Stiffness," Engineering Optimization, Vol. 5, 1981, pp. 169-177. 135 25. S. Adali, "Optimization of a Thin-walled Anisotropic Curved Bar for Maximum Torsioal Stiffness," Journal of Structural Mechanics, Vol . 9, No. 4, Oct. 1981, pp. 389- 413 . 26. J.W. Hou and J.L. Chen, "Shape Optimization of Elastic Hollow Bars," ASME, Journal of Transmission and Automation Design, March, 1985, pp. 100-105. 27. J.W. Hou, J.L. Chen and J.S. Sheen, "Computational Method for Optimization of Structural Shapes," AIAA Journal, Vol. 24, No. 6, 1986, pp. 1005-1012. 28. C . A. Mota Soares, H.C. Rodrogues, L.M. Oliveira Faria and E.J. Haug, "Optimization of the Geometry of Shafts using Boundary Elements," J. Mech. Transm. Autom. , Des. 106, pp. 199-203, 1984. 29. C . A. Mota Soares, H.C. Rodroigues, L.M. Oliveira Faria and E.J. Haug, "Optimization of the Shape of Solid and Hollow Shafts Using Boundary Elements," in Boundary Elements (Ed. by C.A. Brebbia) , Springer-Verlag, pp. 883-889, 1983. 30. C.A. Mota Soares, H.C. Rodroigues, L.M. Oliveira Faria and E.J. Haug, "Boundary Elements in Shape Oprimal Design of Shafts," in Optimization in Computer Aid Design (Ed. by J.S. Gero) , North-Holi and, pp. 155-175, 1985. 31. T. Burczynski and T. Adamczyk, "The Boundary Element Formulation for Multiparameter Structure Shape Optimization," Appl. Math. Modeling, 9, pp 195-200, 1985. 32. J.H. Choi and B.M. Kwak, "Boundary Integral Equation Method for Shape Optimization of Elastic Structures," Int. J. for Num. Meth. in Engrg. , Vol. 26, 1988, pp. 1579-1595. 33. E. Sandgren and S.J. Wu, "Shape Optimization Using The Boundary Element Method with Substructuring," Int. J. Num. Meth. Eng. Vol. 26, pp. 1913-1924, 1988. 34. N. Karmarkar, "A New Ploynomial-Time Algorithm for Linear Programming", Vol. 4, No. 4, Combinatorica , 1984. 35. G.N. Vanderplaats, "Numerical Optimization Techniques for Engineering Design : with Applications," McGraw-Hill, Inc. 1984. 36. M.J.D. Powell, "An Efficient Method for Finding the Minimum of a Function of Several Variables Without Calculating Derivatives," Computer J., Vol 7, No. 4, pp. 303-307, 1964. 136 37. R. Fletcher and C.M. Reeves, "Function Minimization by Conjugate Gradients," Br. Computer J. , Vol . 7, No. 2, pp. 149-154, 1964. 38. R. Fletcher and M.J.D. Powell, "A Rapidly Convergent Method for Minimization," Computer J. , Vol. 6, No. 2, pp. 163-168, 1963. 39. L. A. Schmit and H. Miura, "Approximization Concepts for Efficient Structural Synthesis," NASA CR-2552 , 1976 . 40. G.N. Vanderplaats , "Efficient Algorithm for Numerical Airfoil Optimization," AIAA J. Aircraft, Vol. 16, No. 12, pp. 842-847, Dec 1979. 41. A. V. Fiacco and G.P. McCormick, "Nonlinear Programming : Sequential Unconstrained Minimization Techniques," John Wiley and Sons, New York, 1968. 42. K. Imai and L.A. Schmit, "Configuration Optimization of Trusses," ASCE, J. Struct. Div. , Vol. 107, No. ST5 , pp. 745- 756, May, 1981. 43. J.E. Kelley, "The Cutting Plane Method for Solving Convex Programs," J. SIAM, Vol. 8, pp. 702-712, 1960. 44. G.N. Vanderplaats and F. Moses, "Structural Optimization by Methods of Feasible Directions," J. Computers Struct., Vol. 3, pp. 739-755, July 1973. 45. G. A. Gabriele and K.M. Ragsdell, "The Generalized Reduced Gradient Method : A Reliable Tool for Optimal Design," Transactions of the ASME, May 1977, pp. 394-400. 46. G. Hadley, "Nonlinear and Dynamic Programming," Addison- Wesley, 1964. 47. G.N. Vanderplaats, "CONMIN - A Fortran Program for Constrained Function Minimization User"s Manual," NASA, TM- X-6282 , 1973. 48. C . A. Brebbia, "The Boundary Element Method for Engineering," John Wiley & Sons, New York, 1978. 49. V. Sladek and J. Sladek, "Improved Computation of Stresses Using The Boundary Element Method," Applied Mathematic Modelling, Vol. 10, Aug. 1986. 50. P. Pedersen and C.L. Laursen, " Design for Minimum Stress Concentration by Finite Element and Linear Programming", Journal of Strcture Mechanics. Vol. 10, No. 4, 1983, pp. 375-391. 137 51. B. Parasad and J.F. Emerson, "Optimal Structural Remolding of Multi-Objective System," Computer and Structures, Vol . 18, No. 4, 1984, 1984, pp. 619-628. 52. R.J. Yang and K.K. Choi, "Accuracy of Finite Element Based Shape Design Sensitivity Analysis," ASCE Journal of Structural Mechanics, Vol. 13, No. 2, 1985, pp. 223-239. 53. V. Braibant and C. Fleury, "Shape Optimal Design Using B-splines," Comput. Meths. Appl . Mech. Engrg. 44(3), 1984, pp. 247-267. 54. S.Y. Wang, Y. Sun and R.H. Gallangher, "Sensitivity Analysis in shape Optimization of Continuum Structures," Comput. & Struct., Vol 20, No. 5, pp. 855-867, 1985. 55. V. Braibant, "Shape Sensitivity by Finnite Elements," J. Struct. Mech., 14(2), pp. 209-228, 1986. 56. K.K. Choi and E.J. Haug, "Shape Design Sensitivity Analysis of Elastic Structures," J. of Struct. Mechanics., Vol. 11, No. 2, 1983, pp. 231-269. 57. K.K. Choi, "Shape Design Sensitivity Analysis of Displacement and Stress Constraints," J. of Struct. Mechanics., Vol. 13, No. 1, 1985, pp. 27-41. 58. V. Braibant and C. Fleury, "An Approximation Concepts Approach to Shape Optimal Design," Comput. Meths. Appl. Mech. Engrg. 53(2), 1985, pp. 119-148. 59. N. Kikuchi "Adaptive Grid-Design Method for Finite Element Analysis," Computer Method in Applied Mechanics and Engineering, 55, pp. 129-160, 1986. 60. G.F. Carey and S.Kennon, "Adaptive Mesh Redistribution for A Boundary Element (Panel) Method," Int. J. for Num. Meth. in Engrg., Vol. 24, 1987, pp. 2315-2325. 61. R.G. Hindman and J. Spencer, " A New Approach To Truly Adaptive Generation." AIAA paper 83-0450, AIAA 21st Aerospace Science Meeting, Reno, 1983. 62. R. Fletcher and C.M. Reeves, "Function Minimization by Conjugate Gradients," Br. Computer Journal, Vol. 7, No. 2, pp. 149-154, 1964. 63. C.G. Broydon, "The convergence of A Class of Double Rank Minimization Algorithms," Parts I and II, J. Inst. Math. Appl. Vol. 6, pp. 76-90, 1970. 138 64. R. Fletcher, "A New Approach to Vaeiable Metric Algorithms," Comp. J., Vol. 13, pp. 317-322, 1970. 65. D. Goldfard, "A Family of Variable Metric Methods Derived by Variational Means," Math. Compu. , Vol. 24, pp. 23-36, 1970. 66. D.F. Shanno, "Conditioning of Quasi-Newton Method for Function Minimization," Math. Comput. , Vol. 24, pp 647-656. 1970. 67. P. Hajela and C.-J. Shih, "Multiobjective Optimum Design in Mixed Integer and Discrete Design Variables Problems," Proceedings of the 30th AAIA/ASME/ASCE/AHC/ASC SDM conference, Mobile, Alabama, April 3-5, 1989. 139 BIOGRAPHICAL SKETCH Junhaur Jih was born in Taipei, Taiwan, in 1959. He received a B.S. in Civil Engineering from Chung-Yuan University in 1981. After two years of military service in the army, he came to United States in 1984. He earned his M.S. in the Department of Engineering Mechanics and Engineering Sciences in 1986. In 1988 he was married to Mei-Ming, a doctoral candidate in the department of chemistry. He will receive his Ph.D. from the University of Florida Department of Aerospace Engineering, Mechanics and Engineering Science in August 1989. I certify that I have read this study and that in my opinion it conforms to acceptable standards ofk scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of philosophy. P. Hajela, Chairman Associate Professor of Aerospace Engineering, Mechanics and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as dissertation for the degree of Doctor of Philosophy. X £ L.E. Malvern Professor of Aerospace Engineering, Mechanics and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. X C.T. Sun Professor of Aerospace Engineering, Mechanics and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Assistance Professor of Aerospace Engineering, Mechanics and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a,4d-Ssertation for the degree of Doctor of Philosophy. M. Hoit Assistance Professor of Civil Engineering This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. August, 1989 Cl' Dean, College of Engineering Dean, Graduate School