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