Structural Design for Engineers
Structural Design for Engineers
INTRODUCTION TO
STRUCTURAL DESIGN
design are now at hand. Properly integrated, they can resolve uncertainties and significantly impact
mechanical system design.
An example of an integrated concurrent engineering environment for development of large-scale
wheeled and tracked vehicle systems is illustrated in Fig. 1.1.1 [3]. It comprises simulation and mod-
eling tools and an integration infrastructure to: (1) support design analysis, supportability analysis, op-
eration analysis, and development process control; (2) establish connectivity between all application
tools, with tool interactions transparent to the user; (3) refine product requirements; and (4) conduct
trade-off analyses and make informed decisions to yield a robust optimal design. The integrated test-
bed shown supports concurrent design, operator-in-the-loop driving simulation, dynamic performance
analysis, durability prediction, structural design sensitivity analysis and optimization, maintainability
analysis, and design process management. The test-bed permits all members of the development team
to simulate the performance and effectiveness of product and process designs, at a level of fidelity com-
parable to that, which would be achieved in physical prototyping.
Using the test-bed, an integrated simulation-based design process for the fatigue life of vehicle
components can be developed [3], as shown in Fig. 1.1.2. The process includes CAD-modeling, dynamic
analysis, fatigue analysis, design sensitivity analysis, and design optimization. The CAD-based design
model is critically important for multidisciplinary analysis and design optimization. The process allows
engineers to create a CAD model of a vehicle system and automatically translate the CAD model into
a dynamics model. Dynamic simulations of the vehicle model are then carried out over typical road
profiles to obtain load histories at selected components. In the meantime, CAD models of the selected
vehicle components are created for design parameterization and translated into finite element analysis
(FEA) models. The computation of the fatigue life of a component consists of two parts: dynamic stress
computation and fatigue life prediction. The dynamic stress can be obtained either from experiments
(mounting sensors or transducers on a physical component) or from simulation. Fatigue analysis is
performed using the low cycle fatigue approach. Once the fatigue life of the vehicle components is
obtained, design sensitivity analysis with respect to shape design variables defined in the CAD model is
performed. With the design sensitivity information, design optimization can be carried out to obtain an
optimum design.
As shown in Figs. 1.1.1 and 1.1.2, modern developments of structural design are closely related to
concurrent engineering environments by which multidisciplinary simulation, design, and manufacturing
are possible. Even though concurrent engineering is not the focus of this text, we want to emphasize
structural design as a component of concurrent engineering. Figure 1.1.1 shows an example of concur-
rent engineering environments used in structural design. An important feature of Fig. 1.1.1 is database
management using the CAD tool. Structural modeling and most interfaces are achieved using the CAD
tool. Thus, design parameterization and structural model updates have to be carried out in the CAD
model. Through the design parameterization, CAD, CAE, and CAM procedures are interrelated to form
a concurrent engineering environment.
The engineering design of the structural system in the simulation-based design process consists
of: structural modeling, design parameterization, structural analysis, design problem definition, design
sensitivity analysis, and design optimization. Figure 1.1.3 is a flow chart of the structural design process
in which computational analysis and mathematical programming play essential roles in the design. The
success of the system-level, simulation-based design process shown in Fig. 1.1.2 strongly depends on a
consistent design parameterization, an accurate structural and design sensitivity analysis, and an efficient
mathematical programming algorithm.
A design engineer simplifies the physical engineering problem into a mathematical model that can
represent the physical problem up to the desired level of accuracy. A mathematical model has param-
eters that are related to the system parameters of the physical problem. A design engineer identifies
those design variables to be used during the design process. Design parameterization, which allows the
design engineer to define the geometric properties for each design component of the structural system
being designed, is one of the most important steps in the structural design process. The principal role
of design parameterization is to define the geometric parameters that characterize the structural model
1.1 Elements of Structural Design 5
Physical engineering
problem
Structural modeling
Design parameterization
Performance definition
(cost, constraints)
Structural analysis
Structural model
(FEM, BEM, CFD…)
update
Design sensitivity
analysis
Optimized?
No
Yes
Stop
and to collect a subset of the geometric parameters as design variables. Design parameterization forces
engineering teams in design, analysis, and manufacturing to interact at an early design stage, and sup-
ports a unified design variable set to be used as the common ground to carry out all analysis, design, and
manufacturing processes. Only proper design parameterization will yield a good optimum design, since
the optimization algorithm will search within a design space that is defined for the design problem. The
design space is defined by the type, number, and range of design variables. Depending on whether it is a
concept or detailed design, selected design variables could be non-CAD based parameters. An example
of such a design variable is a tire stiffness characteristic in vehicle dynamics during an early vehicle
design stage.
Structural analysis can be carried out using experiments in actual or reduced scale, which is a
straightforward and still prevalent method for industrial applications. However, the expense and the inef-
ficiency involved in fabricating prototypes make this approach difficult to apply. The analytical method
may resolve these difficulties, since it approximates the structural problem as a mathematical model and
solves it in a simplified form. In this text, a mathematical model is used to evaluate the performance
measures of a structural problem. However, the analytical method has limitations even for very simple
structural problems.
With the emergence of various computational capabilities, most analytical approaches to mathemat-
ical problems have been converted to numerical approaches, which are able to solve very complicated,
real engineering applications. Finite element analysis (FEA), boundary element analysis (BEA), and
meshfree analysis are a short list of mathematical tools used in structural analysis. The development of
FEA is one of the most remarkable successes in structural analysis. The governing differential equation
of the structural problem is converted to its integral form and then solved using FEA. Vast amounts
of literature are published regarding FEA; for example, refer to [4] and references therein. The com-
plex structural domain is discretized by a set of non-overlapping, simple-shaped finite elements, and an
equilibrium condition is imposed on each element. By solving a linear system of matrix equations, the
performance measures of a structure are computed in the approximated domain. The accuracy of the
1.1 Elements of Structural Design 7
approximated solution can be improved by reducing the size of finite elements and/or increasing the
order of approximation within an element.
Selection of a design space and an analysis method must be carefully determined since the analysis,
both in terms of accuracy and efficiency, must be able to handle all possible designs in the chosen design
space. That is, the larger the design space, the more sophisticated the analysis capability must be. For
example, if larger shape design changes are expected during design optimization, mesh distortion in
FEA could be a serious problem and a finite element model that can handle large shape design changes
must be used.
A performance measure in a simulation-based design is the result of structural analysis. Based on
the evaluation of analysis results, such engineering concerns as high stress, clearance, natural frequency,
or mass can be identified as performance measures for design improvement. Typical examples of perfor-
mance measures are mass, volume, displacement, stress, compliance, buckling, natural frequency, noise,
fatigue life, and crashworthiness. A definition of performance measures permits the design engineer to
specify the structural performance from which the sensitivity information can be computed.
Cost and constraints can be defined by combining certain performance measures with appropri-
ate constraint bounds for interactive design optimization. Cost function, sometimes called the objective
function, is minimized (or maximized) during optimization. Selection of a proper cost function is an im-
portant decision in the design process. A valid cost function has to be influenced by the design variables
of the problem; otherwise, it is not possible to reduce the cost by changing the design. In many situa-
tions, an obvious cost function can be identified. In other situations, the cost function is a combination
of different structural performance measures. This is called a multi-objective cost function.
Constraint functions are the criteria that the system has to satisfy for each feasible design. Among
all design ranges, those that satisfy the constraint functions are candidates for the optimum design.
For example, a design engineer may want to design a bridge whose weight is minimized and whose
maximum stress is less than the yield stress. In this case, the cost function, or weight, is the most
important criterion to be minimized. However, as long as stress, or constraint, is less than the yield
stress, the stress level is not important.
Design sensitivity analysis, which is a main focus of this text, is used to compute the sensitiv-
ity of performance measures with respect to design variables. This is one of the most expensive and
complicated procedures in the structural optimization process. Structural design sensitivity analysis is
concerned with the relationship between design variables available to the engineer and the structural
response determined by the laws of mechanics. Design sensitivity information provides a quantitative
estimate of desirable design change, even if a systematic design optimization method is not used. Based
on the design sensitivity results, a design engineer can decide on the direction and amount of design
change needed to improve the performance measures. In addition, design sensitivity information can
provide answers to ‘what if’ questions by predicting performance measure perturbations when the per-
turbations of design variables are provided.
Substantial literature has emerged in the field of structural design sensitivity analysis [5]. Design
sensitivity analysis of structural systems and machine components has emerged as a much needed de-
sign tool, not only because of its role in optimization algorithms, but also because design sensitivity
information can be used in a computer–aided engineering environment for early product trade–off in a
concurrent design process.
Recently, the advent of powerful graphics–based engineering workstations with increasing compu-
tational power has created an ideal environment for making interactive design optimization a viable
alternative to more monolithic batch–based design optimization. This environment integrates design
processes by letting the design engineer create a geometrical model, build a finite element model, param-
eterize the geometric model, perform FEA, visualize FEA results, characterize performance measures,
and carry out design sensitivity analysis and optimization.
Design sensitivity information can be used during a post–processing of the interactive design pro-
cess. The principal objective of the post–processing design stage is to utilize the design sensitivity in-
formation to improve the design. Figure 1.1.4 shows the four–step interactive design process: (1) to
8 1 INTRODUCTION TO STRUCTURAL DESIGN
visually display design sensitivity information, (2) to carry out what–if studies, (3) to make trade–off
determinations, and (4) to execute interactive design optimization. The first three design steps, which
are interactive modes of the design process, help the design engineer improve the design by providing
structural behavior information at the current design stage. The last design step, which could be either
interactive or a batch mode of the design process, launches a mathematical programming algorithm to
perform design optimization. Depending on the design problem, the design engineer could use some or
all of the four design steps to improve the design at each iterative step. As a result, new designs could
be obtained from what–if, trade–off, or interactive optimization design steps.
For the purposes of design optimization, a mathematical programming technique is often used to
find an optimum design that can best improve the cost function within a feasible region. Mathemat-
ical programming generates a set of design variables that require performance values from structural
analysis and sensitivity information from design sensitivity analysis to find an optimum design. Thus,
the structural model has to be updated for a different set of design variables supplied by mathematical
programming. If the cost function reaches a minimum with all constraint requirements satisfied, then an
optimum design is obtained.
In the following sections, each element of the design process is discussed in detail.
Structural Analysis
Design Sensitivity
Analysis
Yes
Stop
When engineers analyze a structural problem, they need to convert the physical problem into a mathe-
matical representation. Many analysis tools can be used to solve this ideal mathematical problem. After
arriving at a solution of the mathematical problem, the meaning of the solution has to be correctly inter-
preted in its physical sense. Thus, if there is an error in the mathematical representation of the physical
1.2 Structural Modeling and Design Parameterization 9
problem, then it is impossible to properly analyze the physical problem no matter what analysis tools
are used. This mathematical representation of the structural problem is called structural modeling.
The reliability of the analysis results strongly depends on the assumptions and idealization used
in structural modeling. However, a too–complex representation of the physical problem may make it
difficult to solve the mathematical problem. For example, when an engineer wants to determine the
height and width of a bridge, it may not be important to model every bolt, because the desired results
will consist of global flexibility and the bridge’s maximum degree of deflection. If each bolt is modeled
for the entire bridge structure, then the analysis cost dramatically increases. It may be presumed that
sections are constructed continuously without any breaks. However, after determining the size of the
bridge, the engineer may want to design the number and size of each section of the bridge. The maximum
magnitude of load carried by each bolt would then be of major concern. In this case, the size and
distance between bolts would be important and structural modeling would need to include bolt strength.
Consequently, different concerns require different structural models. It is the engineer’s responsibility
to find an appropriate trade–off between accuracy and computational costs of analysis.
In structural modeling, the physical problem is represented by mathematical expressions, which contain
parameters for defining that problem. For example, the cantilever beam in Fig. 1.2.1 has parameters
including the length l, the radius of cross-section r, and Young’s modulus E. These parameters define the
system and are called design variables. If design variables are determined, then the structural problem
can be analyzed. Obviously, different design variable values usually yield different analysis results. The
aim of the structural design process is to find the values of design variables that satisfy all requirements.
All design variables must satisfy the physical requirements of the problem. For example, length l of
the cantilever beam in Fig. 1.2.1 cannot have a negative value. Physical requirements define the design
variable bounds. Valid design variables may have to take into account various manufacturing require-
ments. For example, the radius of a cantilever beam satisfies its physical requirement if r is a positive
number. However, in real applications, the circular cross-sectional beam may not be manufactured if its
radius is bigger than r0 . Thus, the range of feasible design can be stated as 0 < r ≤ r0 . 1 In addition, the
design engineer may want to impose certain design constraints on the problem. For example, the max-
imum stress of the beam may not exceed σ 0 and the maximum tip displacement of the beam must not
be greater than z 0 . A set of design variables that satisfy the constraints is called a feasible design, while
a set that does not satisfy constraints is called an infeasible design. It is difficult to determine whether
a current design is feasible, unless the structural problem is analyzed. For complicated structural prob-
lems, it may not be simple to choose the appropriate design constraints so that the feasible region is not
empty.
F
E
r
There are two types of design variables: continuous and discrete. Many design optimization algo-
rithms consider design variables to be continuous. In this text, we presume that all design variables are
1
In general the bounds of design variables are denoted as rL ≤ r ≤ rU where rL is called the lower bound and
rU is called the upper bound, respectively.
10 1 INTRODUCTION TO STRUCTURAL DESIGN
continuous within their lower and upper bound limits. However, discrete design problems often appear
in real engineering problems. For example, due to manufacturing limitations, the structural components
of many engineering systems are only available in fixed shapes and sizes. Discrete design variables can
be thought of as continuous design variables with constraints. As a result, it is more expensive to ob-
tain an optimum design for a problem with discrete design variables. It is possible, however, to solve
the problem assuming continuous design variables. After obtaining an optimum solution for the design
problem, the nearest discrete values of the optimum design variables can be tested for feasibility. If the
nearest discrete design variables are not feasible, then several iterations can be carried out to find the
nearest feasible design.
It is convenient to classify design variables according to their characteristics. In the design of struc-
tural systems made of truss, beam, membrane, shell, and elastic solid members, there are five kinds of
design variables: material property design variables such as Young’s modulus; sizing design variables
such as thickness and cross-sectional area; shape design variables such as length and geometric shape;
configuration design variables such as orientation and location of structural components; and topological
design variables.
In structural modeling, the material property is used as a parameter of the structural problem. Young’s
modulus and Poisson’s ratio, for example, are required in the linear elastic problem. If these material
properties are subject to change, then they are called material propertydesign variables. These kinds of
design variables do not appear in regular design problems, since in most cases material properties are
presumed to be constant. Analysis using such a constant material property is called the deterministic
approach. Another approach uses probability and assumes that material properties are not constant but
distributed within certain ranges. This is called the probabilistic approach and is more practical, since
a number of experiments will usually yield a number of different test results. In this case, material
properties are no longer considered to be constant and can therefore be used as design variables.
Sizing design variables are related to the geometric parameter of the structure. For example, most au-
tomotive and airplane parts are made from plate/shell components. It is natural that a design engineer
wants to change the thickness (or gauge) of the plate/shell structure in order to reduce the weight of the
vehicle. For a structural model, plate thickness is considered a parameter. However, the global geometry
of the structure does not change. Plate thickness can be considered a sizingdesign variable. The sizing
design variable is similar to the material property design variable in the sense that both variables change
the parameters of the structural problem.
Another important type of sizing design variable is the cross-sectional geometry of the beam and
truss. Figure 1.2.2 provides some examples of the shapes and parameters that define these cross-sections.
In the structural analysis of truss, for example, the cross-sectional area is required as a parameter of the
problem. If a rectangular cross-section is used, then the area would be defined as A = b × h. Thus,
without any loss of generality, b and h can be considered design variables of the design problem. De-
tailed discussions of sizing design problems are discussed in Chapter 5 using the distributed parameter
approach.
While material property and the sizing design variables are related to the parameters of the structural
problem, the shape design variable is related to the structure’s geometry. The shape of the structure
does not explicitly appear as a parameter in the structural formulation. Although the design variables
in Fig. 1.2.2 determine the cross-sectional shape, they are not shape design variables, since these cross-
sectional shapes are considered parameters in the structural problem. However, the length of the truss or
1.2 Structural Modeling and Design Parameterization 11
t b b
r b r t
h h
h
t w
w
beam should be treated as a shape design variable. Usually, the shape design variable defines the domain
of integration in structural analysis. Thus, it is not possible to extract shape design variables from a
structural model and to use them as sizing design variables.
Consider a rectangular block with a slot, as presented in Fig. 1.2.3. The location and size of the slot
is determined by the geometric values of Cx , Cy , Dy , r1 , and r2 , which are shape design variables.
Different values of shape design variables yield different structural shapes. However, these shape design
variables do not explicitly appear in the structural problem. If the finite element method is used to
perform structural analysis, then integration is carried out over the structural domain (the gray area),
which is the shape design variable. Since shape design variables do not explicitly appear in the structural
problem, the shape design problem is more difficult to solve than the sizing design problem. Detailed
discussions of the shape design problem are presented in Chapter 6 using the material derivative concept
of continuum mechanics.
Cy Cy
Cx r1 Cx r1
r2 r2
Dy Dy
For those built-up structures made of truss, beam, and shell components, there is another type of de-
sign variable in addition to shape design called the configuration design variable, which is related to the
structural component’s orientation. These components have local coordinate systems fixed on the body
of the structure, and state variables of the problem are described in local coordinate systems. If several
different components are connected together for the built-up structure, the state variables described in
the local coordinate system are transformed to the global coordinate system. If the structural compo-
nents change their orientation in space, the transformation between the local and global coordinates also
changes. Thus, this transformation can be considered the configuration design variable. Since configu-
ration design variables are defined for built-up structures, they are inherently coupled with shape design
variables. That is, in order to allow one member of the built-up structure to rotate, another member’s
12 1 INTRODUCTION TO STRUCTURAL DESIGN
shape needs to be changed. The configuration design variable is not applicable to solid components in
which all rotations can be expressed in terms of shape changes. A simple configuration design variable
will be explained using the example of a three-bar truss in Section 1.2.3. More detailed discussions of
the configuration design problem are presented in Chapter 7 using the material derivative concept in
continuum mechanics.
If shape and configuration design variables represent changes in structural geometry and orientation,
then topology design determines the structure’s layout. For example, in Fig. 1.2.3, shape design can
change the size and location of the slot within the block. However, shape design cannot completely
remove the slot from the block, or introduce a new slot. Topology design determines whether the slot
can be removed or an additional slot is required.
The choice of the topology design variable is nontrivial compared to other design variables. Which
parameter is capable of representing the birth or death of the structural layout? Early developments in
topology design focused on truss structures. For a given set of points in space, design engineers tried to
connect these points using truss structures, in order to find the best layout to support the largest load.
Thus, the on–off types of topology design variables are used. These kinds of designs, however, could
turn out to be discontinuous and unstable.
Recent developments in topology design are strongly related to finite element analysis. The candidate
design domain is modeled using finite elements, and then the material property of each element is
controlled. If it is necessary to remove a certain region, then the material property value (e.g., Young’s
modulus) will approach zero, such that there will be no structural contribution from the removed region.
Thus, material property design variables could be used for the purpose of topology design variables. The
on–off type of design variable can be approximated by using continuous polynomials in order to remove
the difficulties associated with discrete design variables.
In many applications, topology design is used at the concept design stage such that the layout of the
structure is determined. After the layout is determined, sizing and shape designs are used to determine
the detailed geometry of the structure.
A final comment on design parameterization: it is desirable to have a linearly independent set of de-
sign variables. If one does not, then relations between design variables must be imposed as constraints,
which may make the design optimization process expensive, as the number of design variables and con-
straints increase. Furthermore, if design variable constraints are not properly established, meaningless
design results will be obtained after an extensive amount of computational effort. As mentioned before,
this problem is strongly related to structural modeling, since a well–defined structural model should
have an independent set of parameters to define the entire system. Even if defining a good model is not
an easy task for a complicated design problem, the design engineer nevertheless has to define a proper
and independent set of parameters as much as possible in the structural modeling stage.
In this section, a simple example is introduced to discuss design parameterization, which includes ma-
terial property, sizing, shape, and configuration design variables. This example will be used repeatedly
in subsequent sections to explain the structural analysis and design process. The three–bar truss consists
of three truss components, as shown in Fig. 1.2.3.
For truss components, only one material parameter is involved, which is Young’s modulus E. Thus,
the material design parameter is u = [E]. On the other hand, the cross-sectional area of each component
can be chosen to represent the sizing design variables, stated as u = [b1 , b2 , b3 ]T . As explained in
Fig. 1.2.2, the dimensions that determine the cross-sectional shape can be represented as a sizing design.
However, for general truss structures, it is purely a matter of convenience whether the cross-sectional
dimensions or the cross-sectional area is chosen as the sizing design. As far as analysis is concerned,
1.2 Structural Modeling and Design Parameterization 13
only cross-sectional area information is required. For example, if the cross-section of component 1 is
a solid circular shape with radius r, then the relation between the cross-sectional area and the radius
would be b1 = πr2 .
z4
z3
α
b2 θ z6
z5
l b1
b3
z1, f1
z2, f2
δx4 δx6
δx3 δx5
θ+δθ
δl
δx1
δx2
Fig. 1.2.5. Shape and Configuration Design Variable of Truss
As previously pointed out, shape and configuration design variables are closely coupled in the three–
bar truss structure. The shape and configuration design variable is
(Fig. 1.2.5). The change of length involves the shape design, while the rotation of a member affects the
configuration design.
x, z
E, A(x)
F1 F2
f(x)
l
For given x, we assume that stress σ is constant over the cross-section A(x). Thus, stress is a function
of x alone. The same assumption is given to strain ε such that ε(x) = dz/dx. The stress–strain relation
is linear elastic, stated as
When the structure is in equilibrium, the forces generated by the structural deformation are the same
as the externally applied loads. It is equally true that the total potential energy in Eq. (1.3.4) becomes
stationary, which means that the first variation vanishes, so that
1.4 Finite Element Analysis 15
Z l µ ¶µ ¶ Z l
dz dz̄
δΠ = EA dx − f z̄ dx − F1 z̄(0) − F2 z̄(l) = 0 (1.3.5)
0 dx dx 0
where z̄ is the first–order variation of displacement z. A detailed explanation of the variation is provided
in Chapter 2. For the moment, z̄ can be thought of as a small, arbitrary perturbation of z. If z is fixed at a
point x, then z̄ vanishes at the same point. The structural problem is to solve for z in a way that satisfies
Eq. (1.3.5) for all arbitrary z̄.
If the kinematic boundary conditions are given at some point for displacement z, then the possible
candidates for the solution are limited to those that satisfy the displacement boundary conditions. For
example, if the truss structure in Fig. 1.3.1 is fixed at x = 0, then the candidates for the solution belong
to the following solution space:
© ª
Z = z ∈ H 1 (0, l) | z(0) = 0 (1.3.6)
1
where H is a Sobolev space [6] of order one whose elements are continuous functions, and the first
derivative of the function is square integrable in the domain. For a more detailed discussion of basic
function spaces, refer to Appendix A.2. For the moment, readers can think of H 1 as a space of smooth
functions. In addition, the displacement variation z̄ should satisfy Eq. (1.3.6). Thus, the term F1 z̄(0) in
Eq. (1.3.5) would vanish. Physically, if the displacement is fixed, then there would be no work done to
the structure by the applied load.
The analytical solution to Eq. (1.3.5) is non-trivial, even for a simple built-up structure. For a general
shaped structure, an approximation of Eq. (1.3.5) is required by using finite elements. The finite element
method approximates the domain of the structure as a simple geometry set, and then establishes the
equilibrium conditions for each finite element. By combining all finite elements, a global system of
matrix equations is obtained.
Consider a truss finite element with a constant cross-sectional area as shown in Fig. 1.4.1. For sim-
plicity, the distributed load is removed. Displacement of the truss element is represented by two end–
displacements, namely, z1 and z2 .
x, z1 E, A z2
F1 F2
dz 1
= (z2 − z1 )
dx l
dz̄ 1
= (z̄2 − z̄1 ) (1.4.2)
dx l
in which the second equation is the derivative of the displacement variation. By using Eqs. (1.4.1) and
(1.4.2), the variation of the total potential energy in Eq. (1.3.5) can be written as
EA
δΠ = (z2 − z1 )(z̄2 − z̄1 ) − F1 z̄1 − F2 z̄2 = 0 (1.4.3)
l
for all z̄1 and z̄2 in Z. To express Eq. (1.4.3) systematically, it is necessary to define the following vectors
· ¸ · ¸ · ¸ · ¸
z1 z̄1 F1 EA 1 −1
z= , z̄ = , f= , k= (1.4.4)
z2 z̄2 F2 l −1 1
where z is the nodal displacement vector and z̄ its variation; f is the nodal force vector; and k is the
element stiffness matrix. The variational equation (1.4.3) for the truss element can be written as
kz = f (1.4.6)
which is called the local finite element equation. Equation (1.4.6) is applicable to one finite element.
However, for a built-up structure, many truss elements are connected together to make a complete struc-
ture. In this case, the local finite element equation in Eq. (1.4.6) has to be combined to construct the
global finite element equation, and this process is called assembly.
To see the assembly and solution processes of the truss element, consider the three–bar truss example
with a multipoint boundary condition, as shown in Fig. 1.2.3. The displacement and load vectors can be
written in the global coordinate system as
zg = [ z1 z2 z3 z4 z5 z6 ]T (1.4.7)
Fg = [ f1 f2 0 0 0 0 ]T (1.4.8)
The element stiffness matrix in the body-fixed local coordinate system must transformed to the
global coordinate system. For this purpose, let di denote the element local coordinate, and qi represent
the globally oriented element coordinate, as illustrated in Fig. 1.4.2. In each element, the transformation
between the body-fixed and the globally oriented coordinate is
d1 = q1 ,
2 2
d = q ,
c s 0 0
−s c 0 0 3
d =
3
0 0 c s q
0 0 −s c
where c= cosθ and s= sinθ. In addition, the globally oriented coordinates are transformed into the global
coordinate zg by using Boolean matrices as
1.4 Finite Element Analysis 17
010000
1 0 0 0 0 0
q1 =
0 0 0 1 0 0 zg
0 0 1 0 0 0
001000
0 0 0 1 0 0
2
q = z
0 0 0 0 1 0 g
0 0 0 0 0 1
100000
0 1 0 0 0 0
q3 =
0 0 0 0 1 0 zg
000001
After the assembly process, the generalized global stiffness matrix is expressed in terms of the global
displacement coordinate zg as
b3 c2 s b3 cs2 0 0 −b3 c2 s −b3 cs2
b3 cs2 b1 + b3 s3 0 −b1 −b3 cs2 −b3 s3
E 0 0 b2 s/c 0 −b2 s/c 0
Kg (u) = (1.4.9)
l 0 −b1 0 b1 0 0
−b3 c s −b3 cs −b2 s/c 0 b2 s/c + b3 c s b3 cs
2 2 2 2
q14 , d 41 q33
θ
q22 , d 22 q42 , d 42
q12 , d12 q32 , d32
1
3
2
q23 d 3
q11 , d11 d 23 1
q12 , d 21 q13
Fig. 1.4.2. Body-Fixed and Globally Oriented Element Coordinate System
18 1 INTRODUCTION TO STRUCTURAL DESIGN
In actual computation, the global variational equation is modified to explicitly eliminate the bound-
ary condition. However, this is not the general case. In many FEA codes, the size of Kg is retained.
Instead of explicitly removing rows and columns corresponding to the boundary conditions, equivalent
relations are substituted to make Kg positive definite. Since z3 and z4 are prescribed, they can be elim-
inated from the variational equation. In addition, since z5 and z6 has a relation, z6 can be expressed in
terms of z5 , as in Eq. (1.4.10).
Consider the case in which θ = 45o and α = 30o, and define the reduced global displacement vector
as z = [z1 , z2 , z5 ]T . Accordingly, the reduced global load vector is defined as F = [f1 , f2 , 0]√
T
. By
removing the third and the fourth rows and columns of Kg , and by substituting the relation z6 = − 3z5 ,
the reduced stiffness matrix in this example would be
√
b3 √ b3 (√3 − 1)b3
E
K(u) = √ √ b3 2 2b1 + b3 ( 3 − 1)b3 (1.4.12)
2 2l ( 3 − 1)b (√3 − 1)b 2√2b + (4 − 2√3)b
3 3 2 3
Kz = F (1.4.13)
If f1 = f2 = 1 and l = 1, then the solution of the reduced global matrix Eq. (1.4.13) is obtained as
h √ √ √ iT
z= 4−2 3 2 2 1− 3
Eb2 + Eb3 0 Eb2
(1.4.14)
The solution method in Eq. (1.4.13) is easier than that of Eq. (1.4.11). However, for general, complex
kinematic constraints, it may not be easy to explicitly construct the reduced matrix K. In addition, many
FEA codes do not generate a reduced matrix K during the solution procedure. Thus, the use of Eq.
(1.4.11) is clear. In the next section, we will discuss how solution z in Eq. (1.4.14) can be changed as a
function of design vector b.
W (r) = πr2 l
where u = r is the radius and l is the length of the beam. If the radius is a design variable, then the
design sensitivity of W with respect to r would be
dW
= 2πrl
dr
This type of function is explicitly dependent on the design, since the function can be explicitly written
in terms of that design. Consequently, only algebraic manipulation is involved and no finite element
analysis is required to obtain the design sensitivity of an explicitly dependent performance measure.
1.5 Structural Design Sensitivity Analysis 19
However, in most cases, a structural performance measure does not explicitly depend on the design.
For example, when the stress of a beam is considered as a performance measure, there is no simple way
to express the design sensitivity of stress explicitly in terms of the design variable r. In the linear elastic
problem, the stress of the structure is determined from the displacement, which is a solution to the finite
element analysis. Thus, the sensitivity of stress σ(z) can be written as
dσ dσ T dz
= (1.5.1)
dr dz dr
where z is the displacement of the beam. Since the expression of stress as a function of displacement is
known, dσ/dz can be easily obtained. The only difficulty is the computation of dz/dr, which is the state
variable (displacement) sensitivity with respect to the design variable r.
When a design engineer wants to compute the design sensitivity of performance measures such as
stress σ(z) in Eq. (1.5.1), structural analysis (finite element analysis, for example) has presumably al-
ready been carried out. Assume that the structural problem is governed by the following linear algebraic
equation
ψ = ψ (z(u), u) (1.5.4)
The sensitivity of ψ can thus be expressed as
¯ ¯T
dψ(z(u), u) ∂ψ ¯¯ ∂ψ ¯¯ dz
= + (1.5.5)
du ∂u ¯z=const ∂z ¯u=const du
The only unknown term in Eq. (1.5.5) is dz/du. Various computational methods to obtain dz/du are
introduced in the following sub–sections.
Various methods employed in design sensitivity analysis are listed in Fig. 1.5.1. Three approaches are
used to obtain the design sensitivity: the approximation, discrete, and continuum approach. In the ap-
proximation approach, design sensitivity is obtained by either the forward finite difference or the central
finite difference method. In the discrete method, design sensitivity is obtained by taking design deriva-
tives of the discrete governing equation. For this process, it is necessary to take the design derivative of
20 1 INTRODUCTION TO STRUCTURAL DESIGN
the stiffness matrix. If this derivative is obtained analytically using the explicit expression of the stiffness
matrix with respect to the design variable, it is an analytical method, since the analytical expressions
of K(u) and f (u) are used. However, if the derivative is obtained using a finite difference method, the
method is called a semi-analytical method. In the continuum approach, the design derivative of the vari-
ational equation is taken before it is discretized. If the structural problem and sensitivity equations are
solved as a continuum problem, then it is called the continuum–continuum method. However, only very
simple, classical problems can be solved analytically. Thus, the continuum sensitivity equation is solved
by discretization in the same way that structural problems are solved. Since differentiation is taken at the
continuum domain and is then followed by discretization, this method is called the continuum–discrete
method. These methods will be explained in detail in the following sections.
Semi–Analytical Method
Discrete Approach
Analytical Method
Continuum–Discrete Method
Continuum Approach
Continuum–Continuum Method
The easiest way to compute sensitivity information of the performance measure is by using the finite
difference method. Different designs yield different analysis results and, thus, different performance
values. The finite difference method actually computes design sensitivity of performance by evaluating
performance measures at different stages in the design process. If u is the current design, then the
analysis results provide the value of performance measure ψ(u). In addition, if the design is perturbed
to u + ∆u, where ∆u represents a small change in the design, then the sensitivity of ψ(u) can be
approximated as
(1.5.6) and (1.5.7) are virtually independent of the problem types considered. Consequently, this method
is still popular in engineering design.
However, sensitivity computation costs become the dominant concern in the design process. If n
represents the number of designs, then n+1 analyses have to be carried out for either the forward or
backward difference method, and 2n+1 analyses are required for the central difference method. For
modern, practical engineering applications, the cost of structural analysis is rather expensive. Thus, this
method is infeasible for large-scale problems containing many design variables.
dψ4/du
dψ3/du
dψ2/du
dψ1/du
u0 u1 u2 u3 u4 u
Fig. 1.5.2. Influence of Step–size in Forward Finite Difference Method
Another major disadvantage of the finite difference method is the accuracy of its sensitivity results.
In Eq. (1.5.6), accurate results can be expected when ∆u approaches zero. Figure 1.5.2 shows some
sensitivity results using the finite difference method. The tangential slope of the curve at u0 is the exact
sensitivity value. Depending on perturbation size, we can see that sensitivity results are quite different.
For a mildly nonlinear performance measure, relatively large perturbation provides a reasonable esti-
mation of sensitivity results. However, for highly nonlinear performances, a large perturbation yields
completely inaccurate results. Thus, the determination of perturbation size greatly affects the sensitiv-
ity result. And even though it may be necessary to choose a very small perturbation, numerical noise
becomes dominant for a too small perturbation size. That is, with a too small perturbation, no reliable
difference can be found in the analysis results. For example, if up to five digits of significant numbers are
valid in a structural analysis, then any design perturbation in the finite difference that is smaller than the
first five significant digits cannot provide meaningful results. As a result, it is very difficult to determine
design perturbation sizes that work for all problems.
Consider the three–bar truss example shown in Fig. 1.2.3. In this example, the finite element matrix
equation, Eq. (1.4.13), can be solved analytically with the solution given in Eq. (1.4.14) as a function of
the design variable vector u = [b1 , b2 , b3 ]T , which is the cross-sectional area of the truss elements. For
simplicity, if the current value of the design is u = [1, 1, 1]T , and E = 1, then the solution becomes
h √ √ √ iT
z = 4 − 2 3 + 2 2, 0, 1− 3 (1.5.8)
Let us compute the design sensitivity of z1 by using the finite difference method. Since the depen-
dence of z1 on the design is explicitly given, z1 can be straightforwardly computed at different design
stages. Table 1.5.1 shows the sensitivities of z1 with different perturbation sizes. As the perturbation
size decreases, the sensitivity value using finite difference method approaches an exact sensitivity value.
In many cases, the central finite difference method is more accurate than the forward/backward finite
22 1 INTRODUCTION TO STRUCTURAL DESIGN
difference method, as shown in Table 1.5.1, although for the central finite difference method, two per-
formance measure evaluations are involved. Note that, in the latter case, since Eq. (1.4.14) is the exact
solution of the matrix equation in Eq. (1.4.13), there is no concern with numerical noise. That is, as
the design perturbation size decreases, the finite difference results will converge to the exact design
sensitivities.
A structural problem is often discretized in finite dimensional space in order to solve complex problems,
as shown with the finite element method in Section 1.4. The discrete method computes the performance
design sensitivity of the discretized problem, where the governing equation is a system of linear equa-
tions, as in Eq. (1.5.2). If the explicit form of the stiffness matrix K(u) and the load vector f (u) are
known, and if solution z of matrix equation Kz = f is obtained, then the design sensitivity of the dis-
placement vector can also be obtained, as
dz df dK
K(u) = − z (1.5.9)
du du du
This is a discrete approach to the analytical method, since the explicit expressions of K(u) and f (u) are
used to obtain design derivatives of the stiffness matrix and load vector. Even if the expression of Eq.
(1.5.9) is in the global system matrix, actual computation of these derivatives can still be carried out on
the element level in order to avoid a massive amount of calculation related to global stiffness matrix K.
An in–depth discussion of this method is presented in Chapter 4.
It is not difficult to compute df /du, since the applied force is usually either independent of the de-
sign, or it has a simple expression. However, the computation of dK/du in Eq. (1.5.9) depends on the
type of problem. In addition, modern advances in the finite element method use numerical integration
in the computation of K. In this case, the explicit expression of K in terms of u may not be available.
Moreover, in the case of the shape design variable, computation of the analytical derivative of the stiff-
ness matrix is quite costly. Because of this, the semi–analytical method is a popular choice for discrete
shape design sensitivity analysis approaches. However, Barthelemy and Haftka [7] show that the semi–
analytical method can have serious accuracy problems for shape design variables in structures modeled
by beam, plate, truss, frame, and solid elements. They found that accuracy problems occur even for a
simple cantilever beam. Moreover, errors in the early stage of approximation multiply during the ma-
trix equation solution phase. As a remedy, Olhoff et al. [8] proposed an exact numerical differentiation
method when the analytical form of the element stiffness matrix is available.
To obtain discrete design sensitivity with respect to sizing design variables, consider the three bar truss
problem. More complicated design variables, namely, shape and configuration, will be considered in
Chapters 6 and 7. Let us begin with the reduced global stiffness matrix given in Eq. (1.4.12). Since the
explicit form of K(u) is given as a function of design variable u, its derivative can be obtained as
1.5 Structural Design Sensitivity Analysis 23
√
000 000 1 1 √3 − 1
dK dK dK 1
= 0 1 0 , = 0 0 0 , = √ √ 1 √ 1 3 −√1 (1.5.10)
db1 db2 db3 2 2
000 001 3−1 3−1 4−2 3
and load vector f is independent of the design, i.e., df /du = 0. The right side of Eq. (1.5.9) can be
obtained by multiplying Eq. (1.5.10) with Eq. (1.5.8) for each design variable. If Fu denotes the right
side of Eq. (1.5.9), then its explicit expression would be
0 0 −1
Fb1 = 0 , Fb2 = √ 0 , Fb3 = −1√ (1.5.11)
0 3−1 1− 3
Since the right side, corresponding to the first design variable, vanishes, the sensitivity of displacement
with respect to b1 also vanishes. By solving Eq. (1.5.9) with respect to dz/du, we obtain
√ √
0 2 3−4 −2 2
dz dz dz
= 0 , = √ 0 , = 0 (1.5.12)
db1 db2 db3
0 3−1 0
which is the same result as the direct computation of sensitivity in Eq. (1.4.14).
Consider the cantilever beam in Fig. 1.5.3 with point load p at x = l. If one finite element is used to
discretize the structure, as shown in Fig. 1.5.3(b), then displacement z(x) can be approximated as
z1
θ1
z(x) = N zg = [ N1 N2 N3 N4 ]
T
z2 (1.5.13)
θ2
where z1 and z2 are nodal displacement, θ1 and θ2 are nodal rotations, and Ni ’s are corresponding shape
functions of the approximation, defined as
where E is Young’s modulus, I = bu3 /12 is the second moment of inertia, δ(x − l) is the Dirac delta
measure which has a value of infinity at x = l and zero otherwise, and Z is the space of kinematically
admissible displacements. For the moment, Z can be thought of as the space of smooth functions that
satisfy the boundary condition z(0) = θ(0) = 0.
To obtain a finite element equation, the approximation in Eq. (1.5.13) is substituted into variational
Eq. (1.5.15) to obtain
p
E
x
u
b
l
z
(a) Continuum Model
1 2
θ1 θ2
z2
z1
(b) Finite Element Model
Fig. 1.5.3. Cantilever Beam
12 6` −12 6`
2
2`2
Kg = EI 6` 4` −6`
`3 −12 −6` 12 −6` (1.5.17)
6` 2`2 −6` 4`2
£ ¤T
Fg = 0 0 p 0
For the cantilever beam shown in Fig. 1.5.3, the boundary condition is given such that z1 (0) = θ1 (0) =
0. Thus, matrix Kg and vector Fg can be reduced only for z2 and θ2, as
· ¸· ¸ · ¸
EI 12 −6` z2 p
Kz ≡ 3 = ≡F (1.5.18)
` −6` 4`2 θ2 0
The solution z, and thus zg , can be obtained by solving Eq. (1.5.18) as
h iT
zg = 0 0 p`3 p`2 (1.5.19)
3EI 2EI
and by using approximation in Eq. (1.5.13), displacement function z(x) can be obtained as
p
z(x) = NT zg = (−x3 + 3`x2 ) (1.5.20)
6EI
Note that z(x) in Eq. (1.5.20) is the exact solution in this special example.
The discrete method of design sensitivity can be obtained by differentiating the finite element matrix
Eq. (1.5.18) with respect to the design. Consider height u of the cross-sectional dimension as a design
variable. The design sensitivity equation can be obtained from Eq. (1.5.18) using a procedure similar to
that in Eq. (1.5.9), as
dz dF dK
K(u)= − z ≡ Fu (1.5.21)
du du du
where dF/du = 0, since F is independent of the design, and dK/du is calculated from Eq. (1.5.18) as
· ¸
dK 3EI 12 −6`
= 3 (1.5.22)
du ` u −6` 4`2
1.5 Structural Design Sensitivity Analysis 25
and the design sensitivity of the displacement vector can be solved from Eq. (1.5.21) as
· 2 ¸· ¸ " p`3 #
dz ` 4` 6` − 3p − EIu
= u = 3p`2
(1.5.24)
du 12EI 6` 12 0 − 2EIu
Since the shape function of the finite element approximation is independent of the design, the interpola-
tion in Eq. (1.5.13) is valid for displacement sensitivity. Thus, the displacement function sensitivity can
be obtained as
dz(x) dz p
= NT = x2 (x − 3`) (1.5.25)
du du 2EIu
Equation (1.5.25) can be verified by directly differentiating the exact solution in Eq. (1.5.20).
In the continuum method, the design derivative of the variational equation (the continuum model of the
structure) is taken before discretization. Since differentiation is taken before any discretization takes
place, this method provides more accurate results than the discrete approach. In addition, profound
mathematical proofs are available regarding the existence and uniqueness of the design sensitivity. Most
discussions in this text focus on the continuum method, in which analytical expressions of design sensi-
tivity are obtained in the continuum setting.
Sizing design variables are distributed parameters of the continuum equation. For shape design vari-
ables, the material derivative concept of continuum mechanics is used to relate variations in structural
shape to the structural performance measures [5]. Using the continuum design sensitivity analysis ap-
proach, design sensitivity expressions are obtained in the form of integrals, with integrands written in
terms of such physical quantities as displacement, stress, strain, and domain shape change. If exact so-
lutions to the continuum equations are used to evaluate these design sensitivity expressions, then this
procedure is referred to as the continuum–continuum method. On the other hand, if approximation meth-
ods such as the finite element, boundary element, or meshfree method are used to evaluate these terms,
then this procedure is called the continuum–discrete method. The continuum–continuum method pro-
vides the exact design sensitivity of the exact model, whereas the continuum–discrete method provides
an approximate design sensitivity of the exact model. When FEA is used to evaluate the structural re-
sponse, then the same discretization method as structural analysis has to be used to compute the design
sensitivity of performance measures in the continuum–discrete method.
Continuum–based design sensitivity analysis is used to differentiate the variational Eq. (1.5.15) for the
cantilever beam discussed in Example 1.5.3. As with the discrete method, the right side of Eq. (1.5.15)
is independent of the design. Let us define differentiation or variation as
dz
z0 = δu (1.5.26)
du
where δu is the amount of perturbation. The left side of Eq. (1.5.15) can be differentiated with respect
to design u as
26 1 INTRODUCTION TO STRUCTURAL DESIGN
"Z # Z Z
l l l
d 0 3EI
EIz,xx z̄,xx dx δu = EIz,xx z̄,xx dx + z,xx z̄,xx δu dx (1.5.27)
du 0 0 0 u
Thus, the continuum–based design sensitivity equation is obtained as, with δu= 1,
Z l Z l
0 3EI
EIz,xx z̄,xx dx = − z,xx z̄,xx dx, ∀z̄ ∈ Z (1.5.28)
0 0 u
which yields the solution z 0 = dz/du. The continuum–continuum method solves Eq. (1.5.28) to obtain
z 0 directly, whereas the continuum–discrete method first discretizes Eq. (1.5.28), following the same
procedure as the finite element method. If the same approximation is used for z 0 as displacement function
z in Eq. (1.5.13), then the left side of Eq. (1.5.28) becomes equivalent to Eq. (1.5.15) by considering zg
as z0g . The discretized design sensitivity equation therefore becomes
As explained in previous sections, the design sensitivity analysis method has been developed along two
fundamentally different paths, as shown in Fig. 1.5.1. In the discrete method, design derivatives of a
discretized structural FEA equation are taken to obtain design sensitivity information. In the contin-
uum method, design derivatives of the variational governing equation are taken to obtain explicit design
sensitivity expressions in integral form with integrands written in terms of the following variations: ma-
terial property, sizing, shape, and configuration design variables, and such natural physical quantities
as displacement, stress, and strain [5]. The explicit design sensitivity expressions are then numerically
evaluated using the analysis results of FEA codes. Unlike the finite difference method, the continuum
method provides accurate design sensitivity information without recourse to the uncertainties of pertur-
bation size. In addition, the continuum method does not require the derivatives of stiffness, mass, and
damping matrices, as with the discrete method shown in Fig. 1.5.4. Another advantage of the continuum
method is that it provides unified, structural design sensitivity analysis capability, so that it is possible to
develop one design sensitivity analysis software system that works with a number of well–established
analysis methods, such as FEA, the boundary element method, the p-method of FEA, and the meshfree
method.
One important advantage of the continuum–based design sensitivity analysis method is that it is
possible to compute the results of design sensitivity analysis of the established FEA/BEA/Meshfree
1.6 Second–order Design Sensitivity Analysis 27
K(b)z=F T
dΩ = z T f B dΩ + z T f S dΓ
Ω Ω ΓS
Continuum Adjoint
Algebraic Adjoint or
State Sensitivity
ψ′ ≈ h( z , )δ u d Ω
Ω
∂
ψ ′ ≈ [ K (b) z ]δ b
∂b
K =G or Other
Analysis Methods
Evaluate Model
Parameter Sensitivity Evaluate Design Sensitivity
codes with respect to the geometric design variables employed by CAD tools. For example, connections
to CAD tools can be made by providing the design sensitivity of performance measures with respect
to those design variables defined on the CAD tool. Using the same CAD-based design parameters in
manufacturing tools lays the foundation for concurrent engineering. And once models are based on
the same CAD tool, an integrated CAD–FEA–DSA system can be used to develop a design tool for a
concurrent engineering method, such that design and manufacturing engineers can perform trade–off
analysis in the early stages of the design process. A connection can also be made to multi-body dynamic
simulation, computational fluid dynamics, and other computer–aided engineering (CAE) tools, if these
tools use the same geometric CAD modeler as explained in Figs. 1.1.1 and 1.1.2.
First–order design sensitivity analysis, which was introduced in the previous section, is a linear approx-
imation of the performance measure in terms of the design. However, if a higher–order approximation
is used, then the accuracy of the approximation obviously increases. Second–order design sensitivity
analysis uses a quadratic formula in order to approximate the performance change. Let us consider a
Taylor series expansion of ψ(u) up to the quadratic terms, as
µ ¶T
dψ 1 d2 ψ
ψ(u + ∆u) ≈ ψ(u) + ∆u + ∆uT 2 ∆u (1.6.1)
du 2 du
This approximation is exact if ψ is a quadratic function of u. For general nonlinear performance, the
quadratic approximation in Eq. (1.6.1) is much more accurate than linear approximation. In Eq. (1.6.1)
the term d2 ψ/du2 is called the second–order design sensitivity of ψ. For a general n dimensional design
vector u, the second–order design sensitivity becomes a n×n symmetric Hessian matrix, which involves
n(n+1)/2 calculations.
28 1 INTRODUCTION TO STRUCTURAL DESIGN
Consider the three–bar truss example given in the previous section. For second–order design sensitivity,
the design derivative is taken from the first–order sensitivity equation in Eq. (1.5.3), to obtain
d2 z d2 f d2 K dK dz
K 2
= 2− z−2 (1.6.2)
du du du2 du du
The derivatives of stiffness matrix with respect to design in Eq. (1.5.10) are constants. Thus, the second-
order derivative of stiffness matrix d2 K/du2 vanishes, along with d2 f /du2 , and Eq. (1.6.2) is simplified
to
d2 z dK dz
K = −2 (1.6.3)
du2 du du
By using Eqs. (1.5.10) and (1.5.12), Eq. (1.6.3) can be solved for d2 z/du2 , yielding
√ √
4(2 − 3)
2
d z d z 4 2
2
2 = 0 √ , = 0 (1.6.4)
db2 db23
2(1 − 3) 0
which can be verified by differentiating the explicit expression of displacement in Eq. (1.4.14) twice
with respect to the design. All other terms are zero.
Second–order sensitivity information is very useful for an optimization algorithm, since quadratic
convergence can be achieved near the solution points if the Hessian information is available. However,
computing second–order sensitivity results in quite large computational costs. Thus, the design engineer
has to decide between computational cost and the optimization algorithm’s efficiency.
The linear programming method can be used when cost and constraints are linear functions of the design
variables [13]. Most structural design problems, however, are nonlinear with respect to their design
1.7 Design Optimization 29
variables. Thus, the linear programming method is not of much use for structural problems. However,
a nonlinear problem can be solved by approximating a sequence of linear problems, which will be
discussed in Section 1.7.3. The standard form of a linear programming problem is
minimize f = cT u
subject to Au = b (1.7.1)
u≥0
where c = [c1 , c2 , . . . , cn ]T is the coefficient of the cost function, A is the m × n matrix, and b is the
m × 1 vector. Inequality constraints can be treated as equality constraints by introducing slack variables.
Since all functions are linear, the feasible regions defined by linear equalities are convex, along with the
cost function.. Thus, if any optimum solution of Eq. (1.7.1) exists, then it is a global minimum solution
of the problem. The reason for introducing the linear problem here is that a very efficient method exists
for solving linear programming problems, namely the simplex method. A positive feature of a linear
programming problem is that the solution always lies on the boundary of the feasible region. Thus, the
simplex method finds a solution by moving each corner point of the convex boundary.
When cost and/or constraints are nonlinear functions of the design, the design problem is called a non-
linear programming method, as contrasted to the linear programming method discussed in the previous
section. Most engineering problems fall into the former category. Because the properties of nonlinear
programming are nonlinear, this method is frequently solved using the numerical, rather than the ana-
lytical, method.
When there are no constraints on the design problem, it is referred to as an unconstrainedoptimization
problem. Even if most engineering problems have constraints, these problems can be transformed into
unconstrained ones by using the penalty method, or the Lagrange multiplier method. The unconstrained
optimization problem sometimes contains the lower and upper limits of a design variable, since this type
of constraint can be treated in simple way. The standard form of an unconstrained optimization problem
can be written as
minimize f (u)
subject to uL U
k ≤ uk ≤ uk , k = 1, · · ·, n (1.7.2)
In the following sub–sections, numerical methods for solving Eq. (1.7.2) are discussed.
The numerical procedure for solving Eq. (1.7.2) is an iterative update of design u. If uk is the value of
the design at the k-th iteration, then the new design at the (k+1)-th iteration can be obtained by
∂f (uk )
dk+1 = − = −∇f (uk ) (1.7.4)
∂uk
30 1 INTRODUCTION TO STRUCTURAL DESIGN
which is the design sensitivity of the cost function. This method suffers from a slow convergence near
the optimum design, since it does not use any information from the previous design, and only first order
information of the function is used. Note that dk and dk+1 are always orthogonal, such that a zigzagging
pattern appears in the optimization process.
The conjugate gradient method developed by Fletcher and Reeves [14] improves the rate of slow con-
vergence in the steepest descent method by using gradient information from the previous iteration. The
difference in this method is the computation of dk+1 in Eq. (1.7.3). The new descent direction is com-
puted by
Newton Method
The previous methods we have examined use first-order information (first-order design sensitivity) of
the cost function to find the optimum design, which is called linear approximation. The Newton method
uses second–order information (second-order design sensitivity) to approximate the cost function as a
quadratic function of the design. The major concern is how to compute the second-order design sensi-
tivity (or Hessian) matrix. Let us define the Hessian matrix as second-order design sensitivity, defined
in Eq. (1.6.1) as
" #
k ∂f (uk )
H(u ) ≡ , i, j = 1, · · ·, n (1.7.7)
∂uki ∂ukj
The new design can then be determined, as
Quasi–Newton Method
Although Newton’s method has a quadratic convergence, the cost of computing the Hessian matrix, and
the lack of a guaranteed convergence, are drawbacks to this method. The quasi–Newton method has an
advantage over the steepest descent method and Newton’s method: it only requires first–order sensitivity
information, and it approximates the Hessian matrix to speed up the convergence.
The DFP (Davidon–Fletcher–Powell [14]) method approximates the inverse of the Hessian matrix
using first–order sensitivity information. By initially assuming that the inverse of the Hessian is the
identity matrix, this method updates the inverse of the Hessian matrix during design iteration. A nice
feature of this method is that the positive definiteness of the Hessian matrix is preserved.
The BFGS (Broydon–Fletcher–Goldfarb–Shanno [15]) method updates the Hessian matrix directly,
rather than updating its inverse as with the DFP method. Starting from the identity matrix, the Hessian
matrix remains positive definite if an exact line search is used.
Most engineering problems have constraints that must be satisfied during the design optimization pro-
cess. These two types of constraints are handled separately: equality and inequality constraints. The
standard form of the design optimization problem in constrained optimization can be written as
minimize f (u)
subject to hi (u) = 0, i = 1, · · ·, p
gj (u) ≤ 0, j = 1, · · ·, m (1.7.10)
L U
ul ≤ ul ≤ ul , l = 1, · · ·, n
The computational method to find a solution to Eq. (1.7.10) has two phases: first, to find a direction d
that can reduce the cost f (u), while correcting for any constraint violations that are violated; and second,
to determine the step size of movement αgn the direction of d.
The SLP method approximates the nonlinear problem as a sequence of linear programming problems
such that the simplex method in Section 1.7.1 may be used to find the solution to each iteration. By
using function values and sensitivity information, the nonlinear problem in Eq. (1.7.10) is linearized in
a similar way as Taylor’s expansion method in the first order, as
Since all functions and their sensitivities at uk are known, the linear programming problem in Eq.
(1.7.11) can be solved using the simplex method for ∆uk . Even if the sensitivity information is not used
to solve a linear programming problem, design sensitivity information is required in order to approxi-
mate the nonlinear problem as a linear one with SLP. In solving Eq. (1.7.11) for ∆uk , the move limit
∆ukL ≤ ∆uk ≤ ∆ukU is critically important for convergence.
32 1 INTRODUCTION TO STRUCTURAL DESIGN
Compared to previous methods, which use first–order sensitivity information to determine the search
direction d, SQP solves a quadratic sub–problem to find that search direction, which has both quadratic
cost and linear constraints.
1
minimize f (uk ) + ∇f T d + dT Hd
2
subject to hi (uk ) + ∇hTi d = 0, i = 1, · · ·, p
k T
gj (u ) + ∇gj d ≤ 0, j = 1, · · ·, m (1.7.12)
This special form of the quadratic problem can be effectively solved, for example, by using the
Kuhn–Tucker condition and the simplex method. Starting from the identity matrix, the Hessian matrix
H is updated at each iteration by using the aforementioned methods in unconstrained optimization
algorithms. The advantage of solving Eq. (1.7.12) in this way is that for positive definite H the problem
is convex and the solution is unique. Moreover, this method does not require the move limit as in SLP.
In the unconstrained optimization process detailed in Section 1.7.2, the descent direction d is obtained
from the cost function sensitivity. When constraints exist, this descent direction has to be modified in
order to include their effect. If constraints are violated, then these constraints are added to the cost
function using a penalty method. Design sensitivity of the penalized cost function combines the effects
of the original cost function and the violated constraint functions.
If the linear approximation of constraints in SQP is substituted for a quadratic approximation, then
the convergence rate of Eq. (1.7.12) will be improved. However, solving the optimization problem for
quadratic cost and constraints is not an easy process. The constrained quasi–Newton method combines
the Hessian information of constraints with the cost function by using the Lagrange multiplier method.
Nevertheless, it is still necessary to compute the constraint function Hessian. The main purpose of the
constrained quasi–Newton method is to approximate the Hessian matrix by using first–order sensitivity
information. The extended cost function is
p
X m
X
L(u, v, w) = f (u) + vi hi (u) + wi gi (u) (1.7.13)
i=1 j=1
where both v = [v1 , v2 , . . . , vp ]T and w = [w1 , w2 , . . . , wm ]T are the Lagrange multipliers for equality
and inequality constraints, respectively. Note that w > 0. Let the second–order design sensitivity of L
be ∇2 L. The extended quadratic programming problem of Eq. (1.7.12) thus becomes
1
minimize f (uk ) + ∇f T d + dT ∇2 Ld
2
subject to hi (uk ) + ∇hTi d = 0, i = 1, · · ·, p
gj (uk ) + ∇gjT d ≤ 0, j = 1, · · ·, m
wl ≥ 0, l = 1, · · ·, m (1.7.14)
1.7 Design Optimization 33
The feasible direction method is designed to allow design movement within the feasible region in each
iteration. Based on the previous design, the updated design reduces the cost function and remains in the
feasible region. Since all designs are feasible, a design at any iteration can be used, even if it is not an
optimum design. Since this method uses the linearization of functions as in SLP, it is difficult to maintain
nonlinear equality constraints. Thus, this approach is used exclusively for inequality constraints. Search
direction d can be found by solving the following linear sub–problem:
minimize β
subject to ∇f T d ≤ β
∇giT d ≤ β, i = 1, · · ·, mactive
−1 ≤ dj ≤ 1, j = 1, · · ·, n (1.7.15)
where mactive is the number of active inequality constraints. After finding a direction d that can reduce
cost function and maintain feasibility, a line search is used to determine step size α.
The feasible direction method solves the linear programming problem to find the direction of the design
change. The gradient projection method, however, uses a simpler method for computing this direction.
The direction obtained by the steepest descent method is projected on the constraint boundary, such that
the new design can move along the constraint boundary. Thus, the direction of the design change reduces
the cost function while maintaining the constraint along its boundary. For a general nonlinear constraint,
however, a small movement along the tangent line of the boundary will violate this constraint. Thus,
in actual implementation, a correction algorithm has to be followed. The gradient projection method
behaves well when the constraint boundary is moderately nonlinear.