Structural Optimization

4 Structural Optimization
4.1 Introduction
In structural optimization, optimization theory is applied to the design of structures. In Chapter 3, it was emphasized that a design problem should be formulated as in Equation 3.1 in order to apply optimization theory to practical design. It may not be easy to precisely define a design problem as in Equation 3.1. Especially, the equality constraints in Equation 3.1c are difficult to define and manipulate. Nevertheless, the optimization method is frequently utilized in structural design because the formulation in Equation 3.1 can incorporate the objective and the constraint conditions of design.
Structural optimization has been employed to reduce the weight of airplane structures. Aerospace engineering is the field where structural optimization is applied most extensively. Recently, application has expanded to other fields. There are several reasons why optimization is actively applied to structural design.
First, the design process for structures is quite similar to the optimization process. Reduction of weight is important in structural design and it is equivalent to the objective function of optimization in that the reduction is similar to minimization of the objective function. When maximization or minimization of a specific response is pursued, it can be treated as an objective function in optimization. Also, various design conditions for structural design can be easily transformed into constraints in optimization.
Second, the finite element method (FEM) can be exploited in structural optimization. In optimization, the objective function and constraints should be defined in definite forms and they must be evaluated precisely. FEM is an analysis method that facilitates this aspect and it is one of the computational methods that are well developed for analysis of a system. Since many commercial systems are available, it is easy to interface an optimization system with the systems of FEM.
The optimization theories in Chapter 3 are directly used in structural optimization. The process to find optimum sizes of structures is explained by the primal method in Section 3.8. With the current sizes, we evaluate the total weight of a structure (objective function) and the stresses (constraints). A linearized subproblem is defined based on the evaluation. A direction vector is calculated
172 Analytic Methods for Design Practice
from the subproblem and the stepsize is determined. Sensitivity information is needed to make the subproblem. In structural optimization, sensitivity analysis is carried out by using the equilibrium equation (generally formed by FEM). Details will be explained later.
There are two approaches in structural optimization. One is the direct method and the other is the indirect method. In the direct method, an optimization algorithm directly uses the responses from FEM for the functions in Equation 3.1. The functions in Equation 3.1 are approximated to explicit functions of design variables or parameters in the indirect method. An optimization algorithm uses the approximated functions. The functions in Equation 3.1 are replaced by linear or nonlinear functions. The approximation process is conducted based on the characteristics of the responses. Sensitivity information is frequently required in the approximation process.
As mentioned earlier, FEM is utilized to obtain the objective function and constraints. FEM can be extremely expensive for large-scale structures, especially the sensitivity analysis. The process of sensitivity analysis can be directly implemented in FEM software development. Much effort is needed in the implementation, but the resulting numerical cost for the analysis is relatively low. The finite difference method can be used for the sensitivity analysis. In this case, the implementation cost is low, but the numerical cost is relatively high because the finite element analysis is repeatedly utilized. In structural optimization, it is crucial to reduce the numbers of the finite element analyses and sensitivity analyses.
Structural optimization is classified depending upon the characteristics of the design variables. They are size optimization, shape optimization and topology optimization. In size optimization, the domain of FEM is not altered by the design variables. The design variables for shape optimization are defined by the domain of FEM. The distribution or existence of materials serves as design variables in topology optimization. The theories of size and shape optimization are similar, but topology optimization uses quite a different theory.
Structural optimization can be classified into continuous design and discrete design based on the characteristics of the design space. In continuous design, design variables are defined in a continuous space and they can have any real values in the design range. In practical design, the design variables may only have some specific set of values. This is referred to as discrete design. In this chapter, only continuous design is considered and discrete design will be discussed in Chapter 6.
4.2 Finite Element Method
FEM is an analysis method to solve a governing equation that is usually expressed by a differential equation. FEM is generally used for stress analysis, vibrational analysis, etc. in solids. It has been expanded to the areas of heat transfer and fluid mechanics. FEM is not explained in detail and it is introduced from the viewpoint of structural optimization.
Structural Optimization 173
4.2.1 Stress and Strain
Stress is the reaction force per unit area. Each stress component in Figure 4.1 yields the following stress vector:
T zx yz xy z y x
][W W W V V V ı (4.1)
At a point in the space, the stresses in Equation 4.1 should satisfy the equilibrium equation in Equation 4.2
0 w w w w w w x xz
xy x F z y x W W V (4.2a) 0 w w w w w w y yz xy y F z
x
y
W W V (4.2b)
0 w w w w w w z yz
xz z F y
x z W W V (4.2c) where ,,y x F F and z F are components of the body force such as the gravitational force.
Suppose T w v u ][ u is a displacement vector at a point. Then the strains
are
Figure 4.1. Stress components
174 Analytic Methods for Design Practice
x
u
x w w
H (4.3a) y
v
y w w H (4.3b) z
w
z w w
H (4.3c) yx xy x
v
y u J J w w
w w
(4.3d) zy yz y
w
z v J J w w
w w (4.3e)
xz zx z
u x w J J w w w w
(4.3f) In vector form
T zx yz xy z y x
][J J J H H H İ (4.4)
The stress–strain relationship is
E İı (4.5)
where E is a matrix defined by Young’s modulus and Poisson’s ratio.
4.2.2 Formulation of the Finite Element Method
In the finite element analysis for solids, a continuous body is divided into finite elements, and displacements and stresses of elements are evaluated by using solid mechanics. The steps for the finite element method are usually as follows: Step 1. The body for stress analysis is divided into various elements as
illustrated in Figure 4.2.
Step 2. The equilibrium equation is applied to each element and transformed
to a matrix equation by the displacements of nodes in the element. The equation for each element is assembled to form a large equation for the entire structure. This is a simultaneous equation where the unknowns are the displacements of the nodes in elements in Figure 4.3.
Step 3. The solution of the simultaneous equation is obtained. The solution
includes the displacements at the nodes.
Step 4. Other responses such as stresses are obtained from the displacements. Other than the above process, there are many other methods or processes in the finite element method.
Structural Optimization 175
The stiffness matrix is made in Step 2. There are a few methods to make a stiffness matrix as follows:
(1) The matrix method: A simultaneous equation is defined from the
equilibrium equation.
(2) The principle of minimum potential energy: The principle that the
condition of minimum potential energy satisfies the equilibrium equation is utilized.
(3) The principle of virtual work: The principle of virtual work is utilized. It
indicates that the sum of the work by infinitesimal virtual displacements satisfying the boundary conditions and the internal strain energy is zero. Table 4.1 shows the processes for making the element stiffness matrix by the above methods. The matrix method can only make stiffness matrices for simple elements such as truss, beam and triangular elements. The principle of minimum potential energy or the principle of virtual work should be utilized to attain the stiffness matrices for modern complicated elements. It is noted that the stiffness matrices for truss, beam and triangular elements have explicit expressions. However, the stiffness matrix is generally expressed by functionals as follows:
V V
T L d ³ EB B k (4.6)
where V is the domain of each element. As shown in Table 4.1, the equilibrium equation of each element is
L L L f z k (4.7)
Figure 4.2.
A structure defined by finite elements
176 Analytic Methods for Design Practice
Figure 4.3. Elements and nodes
(a) Truss or beam element
(b) Rectangular element
(d) Hexahedral solid element
(c) Triangular element
(e) Tetrahedral solid element
T a b l e 4.1. F o r m u l a t i o n p r o c e s s o f t h e f i n i t e e l e m e n t m e t h o d (C h o i 2002)
Structural Optimization 177
M a t r i x m e t h o d P r i n c i p l e o f m i n i m u m p o t e n t i a l e n e r g y P r i n c i p l e o f v i r t u a l w o r k 1. B y u s i n g t h e s h a p e f u n c t i o n N , t h e d i s p l a c e m e n t f i e l d i n a n e l e m e n t i s a p p r o x i m a t e d a s f o l l o w s : L
N z u w h e r e e L R c z i s a d i s p l a c e m e n t v e c t o r a t t h e n o d e s o f a n e l e m e n t
2. T h e a b o v e e q u a t i o n i s d i f f e r e n t i a t e d a c c o r d i n g t o c o o r d i n a t e v a r i a b l e s a n d m a t r i x B i s o b t a i n e d . T h e s t r a i n –d i s p l a c e m e n t r e l a t i o n s h i p i s d e f i n e d a s L
B z İ 3. B y u s i n g m a t r i x E , w h i c h c h a r a c t e r i z e s t h e m a t e r i a l p r o p e r t i e s , t h e s t r e s s f i e l d i s o b t a i n e d a s L
L S z E B z E İı 4. T h e r e l a t i o n s h i p e q u a t i o n b e t w e e n ı a n d L f i s d e f i n e d a s A ıf L w h e r e ı i s t h e s t r e s s f i e l d o n t h e b o u n d a r y o f a n e l e m e n t a n d L f i s t h e e q u i v a l e n t n o d a l f o r c e v e c t o r 4. T h e p o t e n t i a l o f e x t e r n a l f o r c e s i s o b t a i n e d a s P z T L V 5. T h e s t r a i n e n e r g y i s ³ V T d V U İı³ V T d V
E İİL V T T L d V z E B B z ³ =L L T L z k z 21w h e r e ³ V
T L
d V
E B B k 4. U s i n g t h e s h a p e f u n c t i o n N , v i r t u a l d i s p l a c e m e n t s a n d v i r t u a l s t r a i n s a r e L L z B İz N u ,5. T h e p o t e n t i a l o f t h e f o r c e s b y t h e v i r t u a l d i s p l a c e m e n t i s L
T L V f z
T a b l e 4.1. F o r m u l a t i o n p r o c e s s o f t h e f i n i t e e l e m e n t m e t h o d (C h o i 2002) (c o n t i n u e d )
5. F r o m t h e r e l a t i o n s h i p b e t w e e n ı a n d t h e n o d a l d i s p l a c e m e n t v e c t o r ,L z t h e r e l a t i o n s h i p b e t w e e n L f a n d t h e n o d a l d i s p l a c e -m e n t s i s o b t a i n e d a s L L L L z k A E B z f w h e r e A E B k L
6. T h e f u n c t i o n a l o f t h e t o t a l p o t e n t i a l e n e r g y i s o b t a i n e d f r o m t h e s u m o f t h e d e f o r m a t i o n e n e r g y a n d t h e p o t e n t i a l o f t h e e x t e r n a l f o r c e s
7. T h e f u n c t i o n a l i s m i n i m i z e d . I t i s d i f f e r e n t i a t e d w i t h r e s p e c t t o L z a s f o l l o w s : 0z w 3w L
P 0f z z k z L T L L L T L d d
I f T L z d i s e l i m i n a t e d
L L L f z k
6. T h e s t r a i n e n e r g y b y t h e v i r t u a l s t r a i n i s e v a l u a t e d a s ³³ V V T T V V U d d E İİİı L V T T L V z E B B z ³ d L L T L z k z w h e r e ³ V T L V d E B B k
7. U s e t h e p r i n c i p l e o f v i r t u a l w o r k 0 V U 0 L T L L L T L f z z k z I f T L z i s e l i m i n a t e d L L L f z k
178 Analytic Methods for Design Practice
Structural Optimization 179
where 'l L R z is the nodal displacement vector at each element and 'l is the degrees of freedom for each element. The vector L f is a vector of the external loads imposed on the nodes.
The equilibrium equation of the entire structure is obtained by assembling the element equilibrium equation in Equation 4.7. Equation 4.7 is defined in the global coordinate. If the boundary condition is applied to the assembled equation, the equation yields
f Kz (4.8)
where l R z is a displacement vector for all the degrees of freedom and l is the total degrees of freedom. K is an )(l l u matrix and f is an )1(u l column vector and they are the stiffness matrix and the external load vector for the entire structure, respectively.
The governing equation for vibration of structures can be formulated in the finite element method as follows:
0My Ky [ (4.9)
where K is the stiffness matrix element, M is the )(l l u mass matrix, y is the )1(u l eigenvector and [ is the natural frequency or eigenvalue. The method for obtaining the mass matrix M is not explained.
When the finite element method is numerically implemented, the Galerkin method, which is a weighted residual method, is frequently used. The Galerkin method only uses the equilibrium equation (usually a differential equation, not functionals). Therefore, the Galerkin method is convenient because we do not have to use the functionals. The results of the Galerkin method are the same as those from the principle of minimum potential energy or the principle of virtual work. Details of the finite element method are explained in many books (Bathe 1996, Cook et al . 1989).
Example 4.1 [Symmetric Three Bar Truss Structure]
In Figure 3.3, m,1 l q 60T and the load .N 000,40 P Young’s modulus is GPa.207 E Perform the stress analysis for the three bar truss by the finite
element method. The areas of the cross sections are 2431m 100.4 u b b and
242m 100.2 u b .
Solution
The stiffness matrix iL k of the i th element in the local coordinate system is
»¼
º
«¬ª
1111i
i iL L E
b k (4.10)
180 Analytic Methods for Design Practice
where i L is the length of the i th member. Equation 4.10 should be transformed to the stiffness matrix in the global coordinate system. For the transformation, the following transformation matrix T is utilized:
»¼
º
«
» D D
D D sin cos 0
00sin cos T (4.11) where D
is the angle between a truss element and the x -axis of the global coordinate system.
Using the transformation matrix T, the matrix in Equation 4.10 is transformed as follows:
T k T k iL T i (4.12)
where i k is the stiffness matrix in the global coordinate system for the i th element. Using Equation 4.12, the product of the stiffness matrix and the displacement vector in the global coordinate system for each member is as follows:
»¼º«¬ª »»
»»»»»»»»
¼
º
«««««
«««««¬ª¸¹·¨©§ ¸¹·¨©§ ¸¹·¨©§ ¸¹·¨©§ 111143sin 043cos 0043sin 043cos 1111S S S S L E b z k
»»»»¼
º
««««¬ª»»»»¼º««««¬
ª¸¹·¨©§ ¸¹·¨©§ ¸
¹·
¨©§ ¸¹·¨©
§ 221143sin 43cos 000043sin 43cos v u v u S S S S
»»
»
»¼
º««««¬ª»»»»¼º««««¬ª
22111111111111111111112v u v u L E b (4.13)
»¼º«¬ª »»»»»»
»»»»
¼º
««««««««««¬ª¸¹·¨©§ ¸
¹·¨©§ ¸¹
·¨©§ ¸¹·¨©§ 11112sin 02cos 002sin 02cos 2222S S S S L E b z k »»»»¼º««««¬ª»»»»¼º«««
«¬
ª¸¹·¨©§ ¸¹·¨©§ ¸
¹·
¨©§ ¸¹·¨©§ 33112sin 2cos 00002sin 2cos v u v u S S S S »»»
»
¼
º
««««¬ª»»»»»¼º«««««¬ª
33
11221010
0000
1010000
v u v u L E b (4.14)
»¼º«¬ª »»»»»»
»»»»¼º«««««««««
«¬ª¸¹·¨©§ ¸
¹·¨©§ ¸¹
·¨©§ ¸
¹·
¨©§ 11114sin 04cos 004sin 04cos 3333S S S S L E b z k »»»»¼º««««¬ª»»»»¼º«««
«¬
ª¸¹·¨©§ ¸¹·¨©§ ¸
¹·
¨©§ ¸¹·¨©§ 44114sin 4cos 00004sin 4cos v u v u S S S S »»
»
»¼
º««««¬ª»»»»¼º«
«««
»
44113311111111111111112v u v u L E b (4.15)
where i z is the displacement vector of the i th element in the global coordinate system. The product of the global stiffness matrix K and the displacement z is obtained by assembling the element stiffness matrices as follows:
»»»»
»
»»
»»
»»
»»»
»»
»
»»¼º
«
««««««««««««««««««¬ª»»»»»»»»»»»»»»»»»»»¼º«««
««
«
««
«««««««««««¬
ª 44332211333
33
333333333332222
1111111
1
11111111
33332
211113
32211331
1333311113
311
331
1
220
222200002200000000000
0000000000022220000222222000222220222000
022202202v u v u v u v u L b
L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b L b
L b L b L b L b L b L b L b L b L b L
b L b E Kz (4.16)
If we apply the boundary conditions, components are eliminated by lines as shown in Equation 4.16. Then Equation 4.16 yields
»¼
º«¬ª»¼º«¬ª u u u 119422002
10207102v u Kz (4.17) The external load f is
»¼
º«¬ªu u 2/3400002/140000f (4.18) Substituting Equations 4.17 and 4.18 into Equation 4.8, the following equilibrium
equation is obtained:
»¼
º«¬ªu u »¼º«¬ª»¼º«¬ª u u u 2/3400002/14000022002
102071021194v u (4.19) The simultaneous equation in Equation 4.19 yields
»»¼
º««¬ªu u »¼º«¬ª 4411104659.3104160.3v u z (4.20) The stress ı in a member is defined as
L EBz ı (4.21)
where B is the matrix with a differential operator and L z
is the displacement vector of a member in the local coordinate system. The stress of a truss member can also be obtained by the following equation:
b p / V (4.22)
where b is the area. Equation 4.22 gives the same solution as Equation 4.21. Equation 4.22 is utilized here for stress calculation.
Equation 4.23 is the equilibrium equation in the local coordinate for the truss element in Figure 4.4
»¼º«¬ª»¼º«¬ª »¼
º«¬ª21211111L L z z L bE p p (4.23)
Equation 4.23 is solved as
L L
bE
p p z ]11[12
(4.24) If we substitute Equation 4.24 into Equation 4.22, the stress is
L L
E
z ]11[
V (4.25) Equation 4.25 is a detailed expression of Equation 4.21. Equation 4.20 shows the displacement vector z in the global coordinate system. We need the
displacement vector L z
in the local coordinate system in Equation 4.25. The relationship between z and L z is
Tz z L (4.26)
where the transformation matrix T represents the relationship between the two coordinate systems as shown in Equation 4.11. By substituting Equation 4.26 into Equation 4.25, the stress can be obtained by the displacement vector in the global
1
2
p Figure 4.4. Degrees of freedom for a truss element
coordinate system as follows:
i i L
E
Tz ]11[
V (4.27) Using Equation 4.27, the stresses of the members are
»»»»»
¼
º
««««
«¬ªu u »»»»¼º«««
«¬
ª¸¹·¨©§ ¸¹·¨©§ ¸¹·
¨©§ ¸¹·¨©§ 00
104659.3104160.343sin 43cos 000043sin 43cos ]11[44
11S S S S V L E MPa 227.71 (4.28)
»»»»»¼
º«««««¬ªu u »»»»¼º«««
«¬ª¸¹·¨©§ ¸¹·¨©§ ¸¹·
¨©§ ¸¹·¨©§ 00104659.3104160.32sin 2cos 00002sin 2cos ]11[4422S S S S V L E MPa 744.71 (4.29)
»»»»»¼
º«««««¬ªu u »»»»¼º««««¬ª¸¹·¨©§ ¸¹·¨©§ ¸¹·
¨©§ ¸¹·¨©
§ 00104659.3104160.34sin 4cos 00004sin 4cos ]11[44
33S S S S V L E kPa 5.516 (4.30)
The aforementioned analysis details the process for obtaining the stresses for a
truss structure. Generally, the first process in the finite element method is the calculation of the displacements. Other responses such as strains and stresses are obtained from the displacement information.
Example 4.2 [A Structure with Elements for Plane Stress]
Figure 4.5 shows a four node element for plane stress. Loads are imposed as illustrated in Figure 4.5. Calculate the stress of the element. The thickness is mm,1 h Young’s modulus is GPa 200 E and Poisson’s ratio is .25.0 Q
Solution
The stiffness matrix of the four node element with thickness h is
y x h T L d d ³ EB B k (4.31)
The matrix E in Equation 4.31 is
»»»¼º«««¬ª »»»»¼º«««
«¬
ª
4.000
006667.126667.0026667.006667.12100010112E v v v v E E (4.32) Using the s –t natural coordinate system within an element, Equation 4.31 is
t s y x t s t s t s h T L d d |),,,(|),(),(111
1
J B E B k ³³ (4.33)
where the natural coordinate system is an artificial coordinate system to define the
shape function in an isoparametric formulation. The matrix J is the Jacobian matrix between the s–t coordinate and the x –y coordinate systems. Details are explained in the references for the finite element method. Using four point Gaussian quadrature, the integration of Equation 4.33 is approximated as
¦# 4
1
),,,(),(),(i i i i i i i T L y x t s t s t s h J EB B k (4.34)
From Equation 4.34, it is noted that thickness h is the only property that can be a design variable in size optimization. It is separated from the other terms. Therefore, partial differentiation with respect to h is easy and the sensitivity analysis is relatively simple.
Figure 4.5. Four node element
x
y
For this problem, the element stiffness matrix is
»»»
»»»
»
»»»»¼º««
«««
««
«
««
«¬ª 488667.0166667.0044333.0033.0244333.0166667.0288667.0033
.0166667.0488667.0033.0288667.0166667.0244333.0033.0044333.0044333.0033.0488667.0166667.0288667.0033.0244333
.0166667.0033.0288667.0166667.0488667.0033.0044333.0166667.0244333.0244333.0166667.0288667.0033.0488667.0166667.0044333.0033.0166667.0244333.0033.0044333.0166667.0488667.0033.0288667.0288667.0033.0244333.0166667.0044333.0033.0488667.0166667.0033.0044333.0166667.0244333.0033.0288667.0166667.0488667
.0Eh L
k
(4.35)
The stiffness matrix in the global coordinate system can be obtained by using
Equation 4.12. Since the transformation matrix T is the unit matrix and we have one element in this problem, the element stiffness matrix in Equation 4.35 is identical to the assembled matrix in the global coordinate system. If we apply the boundary condition, the final stiffness matrix is
»»
»»¼º«
«««¬
ª 488667.0166667.0288667.0033.0166667.0488667.0033.0044333.0288667.0033.0488667.0166667.0033.0044333
.0166667.0488667.0Eh
K (4.36) Now we solve the following equilibrium equation to obtain displacements:
f Kz (4.37)
where K is the matrix from Equation 4.36 and T v u v u ][332
2
z (4.38)
T ]10000010000000[ f (4.39)
From Equation 4.37, the displacements are obtained as
410]0012.150759.49649.90419.4[ u T z m (4.40)
Plane stresses are obtained. Stresses of different places in an element are
different for a plate or a shell element. Therefore, the stresses at all the nodes are evaluated from Equation 4.21 and the average value is considered as the final stress. The stress is frequently adopted from the stress at the center of an element (s = t = 0). The latter one is used here. In this case, Equation 4.21 becomes
EBz ı (4.41)
If we substitute B (0,0) at the center, Equation 4.32 and Equation 4.40 into Equation 4.41, the following stresses are obtained:
»»»¼º
««
«¬ª 4.000
006667.126667.0026667.006667.1E ı410000012.150759.49649.90419.400412
14
12
14
12
14
12
121021021021
004
1041041041
u »
»»
»»»
»»»
»»
¼º
««
«
««
«««
««
«¬ª»»»»»¼º««««««¬ª
(4.42) From Equation 4.42, the stresses are
T ]068.50545.64726.56[ ıMPa (4.43)
4.3 Formulation of Structural Optimization
To use the optimization theory in a design problem, the design problem should be formulated as in Equation 3.1. In structural optimization, Equations 4.8 and 4.9 are treated as equality constraints. A structural design problem is formulated as follows:
Find 1,, R R R l n [z b (4.44a)
to minimize ),,([z b f (4.44b)
subject to f z b K )( (4.44c)
0y b M y b K )()([ (4.44d)
m j g j ,,1 0, ),,( d [z b (4.44e)
U L b b b d d (4.44f)
where Equation 4.44c is the equilibrium equation in the finite element method and it is assumed to be linear. Equation 4.44d is the eigenvalue equation. In Equation。