Skip to main content

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

See other formats







To Mei-Ming and my parents 


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. 








1.1 Motivation and Objectives 1 

1.2 Scope of Present Work 5 


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 


2.3 Overview of Mathematical Programming 




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.1 Selection of Design Variables 52 

4.2 Optimum Shape Design for Torsional 


4.3 Optimum Shape Design in 2-D Elasticity 


4.4 Discussions 69 




5.1 An Overview of Sensitivity Analysis 

Methods 71 

5.2 Semi-analytical Approach 75 



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 




7.1 Introduction 113 

7.2 Integrated Approach Using the 

Boundary Element Method 115 







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 



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 


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 

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. 




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. 



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 

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 


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 

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 


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. 


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. 


Finally, conclusions drawn from this study and 
recommendations for future research are presented in Chapter 
8 . 


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 



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 


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 

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 


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] 


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 


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. 


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 


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: 


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 


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. 



Figure 2 . 1 

The relation among optimizer, grid generator, 
and analyzer. 


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 


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. 


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 


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 


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 . 


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 


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 


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 ) 


Here, the subscripts z and y identify the gradients with 
respect to the independent and dependent variables, 

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> 



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. 


V r f(x)i = 

> 0 if zi =zj/ 1 ^ 
< 0 if zi =zj/ u ^ 
= 0 otherwise 

i 1,2,. ..,Q 


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. 



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 



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 



• u) dr 


Since the governing equation of the torsional problem is an 


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 


Figure 3.1 Definition of the boundary element method for 
potential type problem. 


Jn u(x) 8 (x,x Q ) dx 

u(x Q ) , x D e o 
1/2 u(x Q ) , x Q e r 


The fundamental solution of the above problem is of the 

w = i 


4 r 


for 3-D 
for 2-D 


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 


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 



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 


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 


f aw 

H ij = Jr dr 


G^j — J* r w dr 


Figure 3.2 Boundary surface assumed to be hemispherical for 
integration purpose. 


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 


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 


Figure 3 . 3 

A multiply-connected domain. 


K = 2 


n u dft + 2 u Q • D 0 


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 



un it force applied at x in the i direction. These functions 
can be expressed as follows, 



t (3-4./) 

ln(— £) -fiij + 

dr dr 

ax. dXjJ 


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 



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) 


where C. . 


can be expressed as follows, 


C. . = S . . - I . . 
13 ID ID 


Here, 1^^ is defined in terms of 9 ]_ and 9 2 (Figure 3.4), and 

is written as, 


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 = 




Substitution of expressions for u and p into equations 


Figure 3.4 

Definition of domain boundary orientation as 
defined by angles 8 ± and 82 - 


(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) 


[A] (x) = (C) 



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. 



" 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 


2Gv aF 3F. aF. 

{ 6., + G( — + — 1 — ) } U k dr 

c i- 2 i^ J a xi axj axj 


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- 


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 






( computational ) 








4 . 05 


3 .3632 

3 . 3749 


2 . 5 

2 . 6 

4 . 0435 




2 . 98 

4 .3964 

4 .4116 



Figure 3.5 (a) Solution of stress function for elliptic 
torsional shaft. 


Figure 3.5 (b) Shear stress distribution along the 


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 



Figure 3.6 Stress distributation along the multiply- 

connected elliptic domain of the torsioal shaft. 


Figure 3 . 6 Continued 


Figure 3. 



Description of cylindrical shell subjected to 
internal pressure P a and external pressure Pfc,. 


Figure 3.8 (a) Stress distribution r „ versus 

cylindrical shell. xx 


Figure 3.8 (b) Stress distribution r versus 

cylindrical shell. ^ 


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 ] . 



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 



Design generated when node locations are used as 
design variables 

Figure 4 . 1 


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 


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 

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 


distribution in the 


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. 


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 


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 


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 


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. 


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 


= 206.98 

Figure 4.5 Final design for torsional shaft without bending 
inertia constraint. 



Figure 4.6 Design variable definition for multiply- 
connected torsional shaft. 


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 





Figure 4.7 Initial and final design for multiply-connected 
torsional shaft. 

10,000 N 






Figure 4.8 Design variable definition and loading 


Figure 4.9 Iteration history for torque arm. 


Figure 4.10 

Initial and final design of fillet (parabolic 

Figure 4.11 Initial and final design of fillet (elliptic 


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 


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. 



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) 



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) 




[K] {u} = {F} 


where [K] is stiffness matrix, (u) is displacement vector, 


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 

(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 


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 


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 


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 




Figure 5.1 Original and perturbed boundaries discretized 
into N segments. 


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 12 

a in 

A 21 

A 22 '• 

■ A 2N 

(X) = • 

C 2 






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. 



A 7 

A 11 

A 7 

A 12 

A ' 


• A 7 
A 1L 

• • A 7 


A ' 


A > 

A 22 

A ' 

A 2K 

. A 7 


. . A 7 
A 2N 


A 7 

A K1 

A ' .... 

A K2 

A ' 

A K2 

• A ' 


A ' 


c' K 

• (x'l = < 


A ' 


A' .... 

A L2 

A ' 


. A 7 


A 7 



A 7 

A N1 

A N2 

A ' 



A 7 




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, 









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 


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. 


Table 5.1 Comparison of operation count for finite difference 
and semi-analytical approaches to compute design 

finite difference approach 


N 2 

— ~ (NDV + 1) + N (NDV + 1) +0 (N) 

semi-analytical approach 


N 2 

— - + N (NDV + 1) + 0 (N) 

* [A] is an N*N matrix. 

* NDV is the number of design variable. 



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 



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 

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 


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 


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) 


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 


Figure 6.1 Definition of boundary and generalized 
curvilinear coordinates. 


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. 




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, 



P • S£ = C 

( 6 . 6 ) 



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 


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. 


t 1 + 7 • f i 

As j_ = 

j 1 + 7 • f 


( 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 


expressed as follows. 

( ASi ) 
( ASi ) 



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 

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 


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 


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 


(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 


N 1 

f2(s) = S 

i=l ck-[ Asi 


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 


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. 


Figure 6.2 Graphical representation of a two-criterion 
function design space. 


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 »» 


• G C 

(3 (/) -HO 

•P ft) *fH 

G • 

(3 TS -P <w tfl -P 

O 43 

g O G 0 -P 0 

•H U 


0 C Q) G G 

-P (0 


TJ g C -p 3 

(0 O 

P 0) 0 O 4-1 

■P P 


H (D C'rl ft 

a a 


<3 -P -P -P rH 

P W 44 o TS O 

m a 
ts <0 


3 1C (1) 3 -H ft 


-p g p n p -p 


u o t?c 

TS (0 


3 T) ft O 

-P C 


P <u -H -p T3 o 

P O 

■P 0 P G CD 

O -H 


W (0 O -H 0 43 


Cft (0 CP 

• (3 

— - W CD a-H 


<3 43 43 W 43 

U P 



n H n >,4-i . 

n > 

•1 — 

• f3 • tP rp O W 

VO <D vo G G D 

vo G 


> -H CD W G 


(DO) > > G CD 


P i — 1 0) O l 

3) 43 O -H T3 



CP-P > 5 CP (0 


•P -P G CD <D p 


P S -H S3 P CP 


( c ) grid adapta ti on 


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 


Figure 6.4(c) Grid adaptation in regions of high stress. 


Table 6.1 Comparison of ct v /(7 ,, between 149 uniform grids and 61 
redistributed grids. 

master nodes 

av/a al] 

case 1 : 

uniform grid with 
149 analytical nodes 

case 2 : 

adaptive grid with 
61 analytical nodes 

















































* a v is von Mises stress. 

* CT a n is allowable stress (62,000 N/cm 2 ) . 


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. 


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. 







H Milt H-B 

■■II 1 1 i 1 1 l"l- 

i — li 1 <-4- 



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— 


Figure 6.5 Grid distribution obtained by using a control 
function of the form + u 2 -ck. 









w i = w 2 = o, 

“1 = 0, w 2 = 1, 

w]_ = 0.5, w 2 = 0.5, 

— 0.333, w 2 =0.667, 


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 


System equations are solved by 
LU-decomposition . 

Upper and lower triangular matrices 
are saved for sensitivity evaluation. 
Compute objective function and constraint 

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. 


45 KN 

Figure 6.7 Definition of geometry and design variables for 
the elastic plate problem. 


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 

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 







• c 



g 0 



<U *H 



rH -p 



43 3 

0 13 



P "P 



ft P 





<u w 



•P -P 



(0 T3 


( 1 ) 

P • 


ftT3 a) 



•P i — 1 



0 P 0 43 


•h trN 

-P 0 -P 


[/) iH 



3 (0 13 


rH -p -p 


<U -P P 



•H 3 



P G 0 -P 


0 -P <W 






C <D 



c o £ 


cn -p .p 

1 — I 

-H T} 



w <d <p 



(D M Ofl 


•a 3 


13 CPTS 

ft • 

'-I G • 


0 E 

3 G -P 

U • 


•H Cr> G 

Cr> (U r-H i— i 

■P -H -P 


3 13 

•h tn tr> 5 u 

G 0 

G <D < 1 ) 

0 ) Jj 


2 O 

ft ft 

<0 43 
















Figure 6.9 Definition of geometry and loading for the 


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. 



" * — I— 

* — * 1* 

* * w 

Initial Design 

m » T- | 

A ,. r * A 

4 ► 

\ ■ 
1 1 
4 1 

* 1 • 1 * 

Cycle 3 

w i —4 




1 * 1 • 

\ ■ 
\ • 

1 ■ 

• • * * ,, 

Final Design - Cycle 5 

Figure 6.10(a) 

Initial design for the fillet problem. 
Design at begining of third cycle. 

Final optimal design for the fillet problem 
obtained at the fifth cycle. 



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 



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 


obtain a solution to the combined analysis and design 

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 


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 


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. 


Find d, v, and to minimize W 

Subject to Av - C = 0 

a j_ - crvi(v) = 0 

° i 

ai all 


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 


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 


h* 10 cm i — 5 cm 

*21 18 

Figure 7.1 A planar fillet design problem. 


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 


Figure 7.2 37 adaptatively redistributed node points 

for fillet problem. 


(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 

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 


Figure 7.3 Final shape design for fillet problem. 


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 



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. 



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 


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 


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. 



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 


Fingure A.l General two-dimensional boundary element 

and local coordinate system defined at the 
stress point. 


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 . 



2 —2/ 

= (sin 3 a + COS 2 a 

1 - 1 / 



(sin 2 a 


COSa - 


COS 2 a T 22 

sina) pi 

COS 3 a ) P2 

(A. 2) 

(A. 3) 
(A. 2) 

(A. 4) 

(A. 5) 



(7 = ( COSa sin 2 a - COS 3 a ) Pi 

X Z 

1 - 1 / 

+ (sin 3 a - COS 2 a sina) P 2 

1 - 1 / 


+ 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 / 



(A. 6) 



sin 2 a e 

(A. 7) 



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. 


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. 

24. S. Adali, "Cross-Sectional Shape of An Anisotropic 
Hollow Bar for Maximum Torsional Stiffness," Engineering 
Optimization, Vol. 5, 1981, pp. 169-177. 


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. 


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. 


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. 


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. 

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. 



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. 


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 


Dean, College of Engineering 

Dean, Graduate School