Maximizing the fundamental eigenfrequency of geometrically nonlinear structures by topology optimization based on element connectivity parameterization

Computers and Structures 88 (2010) 120–133

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Maximizing the fundamental eigenfrequency of geometrically nonlinear structures by topology optimization based on element connectivity parameterization Gil Ho Yoon School of Mechanical Engineering, Kyungpook National University, Republic of Korea

a r t i c l e

i n f o

Article history: Received 10 July 2009 Accepted 20 July 2009 Available online 13 August 2009 Keywords: Topology optimization Internal element connectivity parameterization method Modal analysis Nonlinear structure

a b s t r a c t This paper pertains to the use of topology optimization based on the internal element connectivity parameterization (I-ECP) method for nonlinear dynamic problems. When standard density-based topology optimization methods are used for nonlinear dynamic problems, they typically suffer from two main numerical difficulties, element instability and localized vibration modes. As an alterative approach, the I-ECP method is employed to avoid element instability and a new patch mass model in the I-ECP formulation is developed to control the problem of localized vibration modes. After the I-ECP based formulation is developed, the advantages of the proposed method are checked with several numerical examples. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Topology optimization methods are used for design purposes in civil and mechanical engineering [1–4]. In particular, layouts of linear dynamic structures, whose stiffness and mass matrices are independent on load and displacements, have been topologically optimized to reduce vibration or noise at target frequencies within pre-assigned design limits (see [5–16] and references therein for more details). In other words, by controlling the first eigenfrequency or some of the lowest eigenfrequencies of a linear structure, its dynamic characteristics could be improved indirectly. However, a review of current topology optimization methods for nonlinear dynamic problems highlights the difficulties related to unstable elements in geometrical nonlinear static analysis [2,17–24] and difficulties associated with the highly localized vibration modes in the modal analysis [12–16,25–29]. This study examines how these difficulties can be resolved and investigates the use of topology optimization for maximizing the fundamental eigenfrequency of geometrically nonlinear structures. In order to calculate eigenfrequencies of a geometrically nonlinear structure, it is necessary to perform a two-step numerical analysis; a modal analysis should be performed after a nonlinear static analysis. First, the deformation of a structure subject to a load at which nonlinear responses are detected in the structure should be computed using a standard nonlinear static solver commonly E-mail addresses: [email protected], [email protected] 0045-7949/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2009.07.006

with Newton–Raphson iteration. Second, we perform a modal analysis with the mass matrix and the tangent stiffness matrix computed at the last Newton–Raphson step of the first nonlinear static analysis [30,31]. In density-based topology optimization, flipped elements having negative Jacobian values – referred to as the unstable elements – inevitably appear around void areas simulated by elements with weak Young’s moduli during the Newton– Raphson iteration of the first nonlinear static analysis. Because these unstable elements make the Newton–Raphson iteration diverge, it is difficult to carry out topology optimization with the nonlinear static analysis. These difficulties have been overcome by developing numerical methods such as the modified Newton– Raphson iteration, the element removal and reintroduction method, and the displacement loading method [17,20,21]. In addition to the first difficulty of the unstable elements, the modal analysis for topology optimization problem is also suffering from highly localized vibrating modes in low-density regions having relatively low eigenvalues. Considering one or some of the lowest modes in topology optimization, these highly localized vibrating modes are troublesome. To resolve them, an optimization formulation maximizing localized buckling modes with a mesh independent filter have been proposed [12,26,27,29]. Alternatively several techniques involving mass density interpolation have also been proposed without changing optimization formulation [14,25]. In this research, to overcome the above mentioned difficulties, an alternative topology optimization approach called the internal element connectivity parameterization (I-ECP) method was

121

G.H. Yoon / Computers and Structures 88 (2010) 120–133

employed. Unlike the standard density-based approach, this new method defines a structural layout by defining connectivity among solid elements using zero-length links, as shown in Fig. 1 [19,32–

34]; it may be possible to interpret the element connectivity parameterization (ECP) method as a numerical method imposing weak continuity constraints among elements using the penalty

(a)

Weak links

Design variable Weak element Strong links

Solid element

SIMP approach (b)

Patch Design variable

I-ECP approach (c)

Fig. 1. Modeling by the density-based SIMP and the I-ECP (internal-element connectivity parameterization).

Fig. 2. Solution procedure for the nonlinear modal problem.

122

G.H. Yoon / Computers and Structures 88 (2010) 120–133

or the Lagrange multiplier (see related discussions in [20]). Because the material properties of solid elements do not change during optimization iterations, the unstable element of nonlinear structures can be eliminated. Depending on how solid elements

are connected with zero-length links, two ECP methods exist [32]. One is the external ECP (E-ECP) method in which solid elements are directly connected, and the other is the internal ECP (I-ECP) method, where the solid elements are indirectly connected

Fig. 3. Effects of the applied load on eigenfrequencies for a slender beam. (a) Problem definition, (b) eigenmodes and eigenvalues for zero force, (c) eigenmodes and eigenvalues for nonzero force, and (d) curves of eigenfrequency versus applied load for the beam.

G.H. Yoon / Computers and Structures 88 (2010) 120–133

123

Fig. 4. e-th planar rectangular patch consisting of a plane finite element and four zero-length links.

through outer nodes. Because of the computational advantage in the latter, this study employs the I-ECP method [32,35]. For successful analysis and optimization, this study also develops a novel way of mass modeling formulation for the I-ECP method. So far, because of the ambiguity in the definition of the formulation of a mass matrix, the I-ECP method has not been applied to dynamic problems. It might have been obvious to define a mass matrix for a solid element inside a patch, but it was not apparent that a mass matrix could be defined for zero-length links; here ‘‘define” means calculating and assembling the mass matrix in the framework of the finite element (FE) procedure. Moreover, from several empirical test results discussed in Section 2.3, it appears that assembling the mass matrix of a solid element to the degrees of freedom of the solid element complicates the optimization process due to the presence of highly localized modes in patches

Fig. 5. Mass modelings using the I-ECP method. (a) A straightforward model (assigning the mass matrix at the solid element to the degrees of freedom of the inner nodes), and (b) the proposed patch mass matrix model (assigning the mass matrix to the degrees of freedom of the outer nodes).

[12–14,26–28]. To avoid the complication, an alternative approach in which the mass matrix of a solid element is assembled to the degrees of freedom of the outer nodes is presented in this study. Adopting the new mass modeling method proposed here and the mass interpolation functions proposed in [14,25], the difficulty arising from highly localized modes inside patches could be overcome and topology optimization could be carried out for nonlinear dynamic system. This paper is organized as follows. In Section 2, we provide an overview of the basic notations and governing equations for nonlinear modal analysis. We use the I-ECP method and study the mass modeling of the method. The parameterization of design variables and the sensitivity analysis of the I-ECP method are dealt with in Section 3. In Section 4, two numerical examples in which the fundamental eigenfrequency of two-dimensional structures is maximized are presented to show the potential of the proposed method in Section 4. In Section 5, some observations are made regarding designs of nonlinear dynamics system and future work is discussed.

Fig. 6. Solid element: (a) a finite element model, and (b) four eigenvalues and eigenmodes.

124

G.H. Yoon / Computers and Structures 88 (2010) 120–133

displacements at time t + Dt of a generic point of a body by DU and U, the following update rules can be obtained:

2. Modal analysis of nonlinear structures using the I-ECP method

tþDt

tþDt

2.1. Formulation of nonlinear vibration problem

t

The vibration motion of a structure with practical dynamic loads is represented by sinusoidally varying eigenmodes as well as static equilibrium displacements [30,31,36,37]. In the case of a linear structure, the static equilibrium displacements are zero and only the sinusoidal eigenmodes remain. Thus, the linear eigenvalue problems can be solved by using load-independent stiffness and mass matrices can be solved. In contrast, to calculate eigenvalues and their associated eigenmodes for a nonlinear structure, the modal analysis depicted in Fig. 2 should be used along with the tangent stiffness matrix computed for the deformed domain for given static loads and the mass matrix [30,31]. Consequently, eigenfrequencies of a geometrically nonlinear structure are dependent on current static displacements even for a structure made of linear elastic material. For nonlinear static analysis, the notations given in [30,31] were followed. By respectively denoting the updated displacements and

structure

UðkÞ ¼ tþDt Uðk1Þ þ DUðkÞ ;

ðk1Þ KT DUðkÞ

¼ Rð

tþDt

ðk1Þ

U

Uð0Þ ¼ t U

ð1Þ

Þ ðsee ½19; 20Þ

ð2Þ

tþDt

where the superscript (k) denotes the kth iteration step in the Newton–Raphson method. The incremental residual and tangent stiffðk1Þ ness matrix are denoted by RðtþDt Uðk1Þ Þ and t KT , respectively. After solving the nonlinear static equation (Eq. (2)), the following modal analysis should be carried out, where the displacementdependent tangent stiffness matrix tþDt KT and the displacementindependent mass matrices, M are used

ðtþDt KT  x2 MÞU ¼ 0; UT MU ¼ I

ð3Þ

Here, x and U are eigenfrequencies and associated eigenmodes, respectively. The identity matrix is I is of the same size as the as the mass or stiffness matrix. In order to investigate how eigenvalues of nonlinear structures are affected significantly by applied forces in practice, we considered the simple beam shown in Fig. 3a. The left side of the beam

Fig. 7. Comparison of mass models obtained using the ECP method ðlmax =kdiagonal ¼ 104 Þ. (a) A straight mass model, (b) a present patch mass model, and (c) modes of two methods.

125

G.H. Yoon / Computers and Structures 88 (2010) 120–133

is clamped and a concentrated force is applied to the right side. The deformed shape of the beam determined by nonlinear static analysis in Fig. 3c and the associated eigenmodes obtained from the modal analysis of Eq. (3) indicate that the eigenfrequencies vary nonlinearly with the magnitude of the applied force; the first and second angular speeds at under the applied force of 0.003 N are increased approximately by 130% and 14%. However, the third angular speed is decreased by 8%. As observed in the figure, the eigenfrequencies are influenced significantly by the applied loads and associated deformations. Fig. 3c plots the eigenmodes of the undeformed structure; if the modes were plotted at the deformed structure, it would be difficult to compare them with those of linear analysis as shown in Fig. 2. Thus, later in the document, we present eigenmode plots of examples at undeformed structures.

2.2. Nonlinear static analysis using the ECP method The I-ECP method enables the realization of a layout different from that of the element density-based method. To underline the basic concepts of these two methods and some fundamental differences between them, consider the layout in Fig. 1a, which may be observed during topology optimization iterations. When the standard element density method is employed, design variables defined for each element and their corresponding Young’s moduli are varied to realize the layout of Fig. 1a in terms of the potential energy. Although this method is robust and simple, it suffers from numerical instability called the unstable element with negative

Jacobian values in topology optimization for geometrically nonlinear structures [2,17–21,23,24]. As opposed to the density-based methods, the I-ECP method does not change the material properties of plane or cubic finite elements discretizing a domain during topology optimization as shown in Fig. 1c. To define a layout, all elements are disconnected from other elements and nodes at the same locations are connected using one-dimensional links. To reduce the computation time, the static condensation scheme, which condenses out for the degrees of freedom of inner nodes of the I-ECP patch from the global stiffness matrices, was proposed in [32]. A detail condensation implementation of the I-ECP method can be found in [35]. To formulate the I-ECP method for nonlinear static analysis, the eth patch shown in Fig. 4 is considered along with the assumption of geometrical nonlinearity. The nodes connecting elements are named the outer nodes and those defining plane elements are named the inner nodes. The displacements of the outer nodes ðkÞ ðkÞ and inner nodes are denoted by tþDt ue;out and tþDt ue;in in the kth iteration, respectively. The displacements are then updated as follows:

" tþDt

ðkÞ

ue;out

tþDt

#

ðkÞ

ue;in ðkÞ

" tþDt ¼

tþDt

ðk1Þ

ue;out

ðk1Þ

ue;in

#

" þ

DuðkÞ e;out DuðkÞ e;in

# ð4Þ

ðkÞ

where Due;out and Due;in denote the updated displacements for the outer nodes and the inner nodes, respectively, and are calculated by the following equation:

Fig. 8. Ratios of the condensed stiffness matrix to mass matrix with various mass interpolation functions for one element in Fig. 7 (SIMP penalty: 3). (a) The stiffness matrix behavior, (b) a linear mass interpolation function of Eq. (25), (c) a power mass interpolation function of Eq. (26), and (d) a C1 continuous mass interpolation function of Eq. (27) proposed in [14].

126

G.H. Yoon / Computers and Structures 88 (2010) 120–133

kI;e ¼ le ðce ÞI88 (

ðwhere I88 is a 8  8identity matrixÞ

kI;e

kI;e

kI;e

kI;e

"

 þ

0 0

#)"

0 t

DuðkÞ e;out

structure;ðk1Þ

kT;e

#

DuðkÞ e;in

" ¼

ð5Þ ðk1Þ

Re;out

#

ðk1Þ

Re;in

ð6Þ The link stiffness of the eth patch, le, is a function of the design variable ce. The stiffness matrix and the residual force terms of the outer and the inner nodes of the eth patch are denoted by ðk1Þ ðk1Þ t structure;ðk1Þ kT;e , Re;out and Re;in , respectively. ðk1Þ ðk1Þ For the I-ECP, Re;out and Re;in can be formulated as

"

# ðk1Þ

Re;out

ðk1Þ

Re;in 2 4

tþDt

 tþDt ¼

Re

0

link;ðk1Þ 0 f e;out

tþDt link;ðk1Þ f e;in 0

3 5¼









0 tþDt

kI;e kI;e



structure;ðk1Þ 0fe

kI;e kI;e

" tþDt tþDt

2 4

tþDt link;ðk1Þ 0 f e;out tþDt

link;ðk1Þ

0 f e;in

3 ð7Þ

# ðk1Þ

ue;out

ð8Þ

ðk1Þ

ue;in

In Eq. (7), the externally applied force on the outer nodes and the internal force acting on the inner nodes are denoted by tþDt Re and tþDt structure;ðk1Þ , respectively. Because the degrees of freedom of 0fe the inner nodes are independent of those of the other nodes, the static condensation scheme can be applied. t

ð9Þ ð10Þ

The global tangent matrix is assembled as t

ðk1Þ

KCon ¼

Np X

t

ðk1Þ

kCon;e

h

i structure kT;e  x2 me ½/ ¼ |{z} 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |{z} 81 tþDt

ð11Þ

e¼1

ð14Þ

81

88

For the I-ECP patch of Fig. 7a:

("

kI;e

kI;e

#  x2



0

0

) "

ue;out

#

structure ue;in kI;e kI;e þ tþDt kT;e 0 me |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl} 1616

5

 1 ðk1Þ structure;ðk1Þ kCon;e ¼ ðkI;e  kI;e kI;e þ tþDt kT;e kI;e Þ  1 structure;ðk1Þ ðkÞ ðk1Þ ðk1Þ t ðk1Þ kCon;e Due;out ¼ Re;out þ kI;e kI;e þ t kT;e Re;in

patch assembling the mass matrix of the solid element to the inner nodes is constructed in Fig. 7a. The modal analysis equations of the solid element and the patch are formulated as follows For the solid element of Fig. 6:

¼ |{z} 0 161

161

ð15Þ Here, u, ue,out, and ue,in are the eigenmodes of the solid element and the sub-vectors of the eigenmodes for the outer and the inner nodes of the patch in Fig. 7a, respectively. The stiffness matrix and the structure and me, mass matrix of the solid element are denoted by t kT;e respectively. In Eq. (15), note that the mass matrix of the solid element is assembled for the degrees of freedom of the inner nodes. structure is assembled for the Further, the tangent stiffness matrix t kT;e degrees of freedom of the inner nodes. Because four new nodes are added when constructing a patch in the I-ECP method, there are eight additional modes between the inner and the outer nodes in addition to the eight eigenmodes of the solid element. As seen in Fig. 7c, the first four eigenvalues and eigenmodes of the outer nodes become almost identical with the high stiffness value of the links. In contrast to the case of links with a high stiffness value, highly localized vibrating modes inside the patch of Fig. 7a have been found as side effects when the stiffness value of links becomes small. In other words, with a small stiffness value, the fundamental eigenmode of the patch will be the motion observed between the solid element and links inside the patch. As this phenomenon

where Np is the total number of patches. Then the following system of equations is solved iteratively for DUðkÞ out , which is the global displacement vector for the outer nodes: t

ðk1Þ

ðk1Þ

KCon DUðkÞ out ¼ RCon ðk1Þ

where RCon ¼

Np  X

ðk1Þ



structure;ðk1Þ

Re;out þ kI;e kI;e þ t kT;e

1

ðk1Þ

ð12Þ 

Re;in

e¼1

ð13Þ 2.3. Modeling of mass matrix in I-ECP method for modal analysis Eigenfrequencies of the nonlinear system that can be determined using the I-ECP method are also calculated by modal analysis (Eq. (3)). The tangent stiffness matrix can be easily computed by Eq. (2). But it is not clear how the mass matrix in Eq. (3) is constructed for the I-ECP method, which is one of the limitations of this method. The construction of the mass matrix is examined in the following sections. 2.3.1. Formulation 1: straightforward method for mass matrix modeling For dynamic analysis, assigning mass matrices to each solid elements is straightforward as shown in Fig. 5a. If the stiffness values of the links of a patch are sufficiently large, some eigenvalues of an I-ECP patch become equivalent to the eigenvalues of a solid element in principle. For example, consider the two-dimensional solid element in Fig. 6a under the assumption of small displacements. Because the degrees of freedom of the bottom nodes are clamped, there are 4 real and 4 infinite eigenvalues in Fig. 6b. Here an I-ECP

Fig. 9. Cantilever problem in the case where the ends of the cantilever are clamped. (a) Problem definition, (b) a symmetric linear result with an evenly distributed initial guess, and (c) an unsymmetrical result with an unsymmetrical initial guess.

G.H. Yoon / Computers and Structures 88 (2010) 120–133

can be observed for any patch with a small stiffness value, many highly localized modes will appear. Furthermore, because the topology optimization formulation for dynamic structures usually involves the consideration of the fundamental eigenfrequency or some of the lowest eigenfrequencies, these localized modes of patches are not desirable. Consequently, although the proposed method, which assigns the mass matrices to solid elements is simple, it is not applicable to topology optimization. 2.3.2. Formulation 2: present mass matrix approach To resolve the above-mentioned problems related to the localized modes in patches, this paper presents a new concept assembling the mass matrices of solid elements in patches into the degrees of freedom of the outer nodes, as shown in Fig. 5b. Unlike the straightforward method, the solid element of the patch is modeled as massless, which makes the eigenmodes between the inner and the outer nodes numerically infinite. To formulate this patch mass matrix approach, Eq. (15) is modified as below.

For the patch mass matrix : # (" kI;e kI;e kI;e

kI;e þ

tþDt

structure kT;e

 x2



me

0

0

0

)"

ue;out ue;in

# ¼0

ð16Þ

127

This approach has the following features. (1) Compared to Eq. (15), the mass matrix of the solid element is simply assembled into the degrees of freedom of the outer nodes, which makes it simple to implement. (2) This simple modification of the mass matrix resolves the local modes observed between the inner and the outer nodes completely. For example, Fig. 7b and c shows the result of the reanalysis of the eigenmodes using this patch matrix approach, where four real eigenvalues and eight infinite eigenvalues are existing. (3) It is important to note that the local modes observed among patches, which correspond to local modes among void elements in element-density-based approaches, still exist. Therefore, links with very small stiffness values can produce highly localized vibrating modes which cause non-convergence in topology optimization. For the efficient construction of the present patch mass model, the following condensation scheme is used for the global stiffness matrix and mass matrix:

ðkI;e ue;out  kI;e ue;in Þ  x2 me ue;out ¼ 0

Fig. 10. Eigenfrequency histories for (a) a design with an initial uniform density, and (b) an optimized design with an unsymmetrical density distribution.

ð17Þ

128

G.H. Yoon / Computers and Structures 88 (2010) 120–133

  structure  kI;e ue;out þ kI;e þ tþDt kT;e ue;in ¼ 0   1  structure;ðk1Þ t ðk1Þ kCon;e ¼ kI;e  kI;e kI;e þ tþDt kT;e kI;e t

KCon ¼

Np X

t

ðk1Þ

kCon;e ; MCon ¼

e¼1

Np X

me

ð19Þ

Material model 2 : ( m0 ce ; ce > 0:1 me ¼ m0 ðc1 c6e þ c2 c7e Þ; ce  0:1; c1 ¼ 6  105 ; c2 ¼ 5  106

ð20Þ

e¼1

ð21Þ

3. Topology optimization formulation

For the sake of simplicity, we consider topology optimization to maximize only the fundamental eigenfrequency frequency of a geometrically nonlinear structure, which is defined as c

For the topology optimization of dynamic problems, the values of two mechanical material properties, i.e., link stiffness and material density, are interpolated with respect to the design variable (c). In this paper, the following interpolation function is used for the link stiffness in Eq. (5) (see [32,35] for more detail):

as structure kdiagonal

a ¼ lmax  lmin ; b ¼ lmin cmin  ce  1; cmin ¼ 0:001

k

! ð22Þ ð23Þ ð24Þ

where k is the number of degrees of freedom per node and s and n are penalties. A diagonal term of the linear stiffness matrix is destructure noted by kdiagonal . The upper and lower bounds of the stiffness of links are denoted by lmax and lmin, respectively. Numerical examples indicate that with sufficiently large and small values for lmax and lmin and appropriate penalties, similar results can be obtained (see [35] for the effects of these links in linear static structure cases and Section 4 for numerical examples in dynamic structure cases). structure structure Here, lmax and lmin are set to 104  kdiagonal and 104  kdiagonal , respectively. In all numerical examples, we assumed n = 3, k = 2, and s = 10 assumed that because optimal layouts similar to those obtained in the solid isotropic material with penalization (SIMP) method could be obtained. Fig. 8a shows the ratio of the condensed stiffness matrix to the normal stiffness. The mass matrix can be interpolated as follows:

me ¼ m0 ce

3.2. Topology optimization formulation

Max

3.1. Material Interpolation

s¼

ð26Þ

ð27Þ

½t KCon  x2 MCon Uout ¼ 0

ce n þb 1 þ ð1  cne Þs

m0 ce ; ce > 0:1 m0 c6e ; ce  0:1

Material model 1 : me ¼

Finally, the following modal equation is solved for obtaining the outer nodes, Uout

le ¼ a



ð18Þ

s:t:



Min fxj g



j¼1;...;j Np X

q c v  V

eð eÞ e e¼1  tþDt

Rð

ð28Þ

UÞ ¼ 0

ðtþDt KT ðtþDt U Þ  x2 MÞU ¼ 0 where qe, ve and V* are the element density, element volume, and prescribed volume limit, respectively. The converged displacements obtained by solving the nonlinear static equation are denoted by t+Dt * U , and the number of candidate eigenfrequencies is denoted by J. To determine optimal topologies using a gradient-based optimizer, the sensitivity analysis of the jth eigenfrequency, xj, should be performed with respect to the design variable. In this study, the following sensitivity equation can be derived for the patch mass method (see [26–28,38]):

  dxj 1 dle dme ¼ ðue;out  ue;in ÞT ðue;out  ue;in Þ  x2j uTe;out ue;out dce 2xj dce dce ð29Þ

ð25Þ

where m0 is the nominal mass stiffness matrix. Numerical tests show that this interpolation function causes highly localized vibration modes inside I-ECP patches having links with low stiffness values [12,26–28].When the present patch mass method is used, it is observed that highly localized vibrating modes do not appear between the inner and outer nodes because the solid plane element inside an I-ECP patch only has stiffness and is massless, as shown in [16]. However, the localized modes continue to appear at the outer nodes, which are locally vibrating when the design variables have low values. That is because the ratio of the condensed stiffness matrix of the Eq. (19) to the mass matrix for the I-ECP patch becomes almost zero, similar to the case of the SIMP method; this results in the presence of localized vibration modes in patches with low link stiffness values (Fig. 8). (See [12,14,26–29] for localized vibration modes and sensitivity analysis of the SIMP method.) To suppress these local modes, some mass interpolation functions that have been verified in [12,14,25] are also tested here. By employing the following material interpolations, the locally vibration modes can be suppressed as expected

Fig. 11. Optimized results considering the geometrical nonlinearity. (a) A model with an arbitrary force, (b) with F = 10 N, and (c) with F = 20 N.

G.H. Yoon / Computers and Structures 88 (2010) 120–133

129

Fig. 12. Eigenfrequency histories for (a) a design with F = 10 N, and (b) a design with F = 20 N.

4. Topology optimization of nonlinear dynamic problems In order to illustrate the potential of the developed method, two-dimensional optimization problems with different loads are considered. Unless stated otherwise, uniform initial guesses of c satisfying given mass constraints are used. The dimension, material properties and magnitude of the loads are arbitrarily chosen to show the effectiveness of the present I-ECP method. For successful optimization, it is important to use a proper optimization algorithm such as sequential linear programming, sequential quadratic programming, and the method of moving asymptotes. In this study, the method of moving asymptotes is used as an optimization algorithm; however other optimization algorithms can also be used [39]. To get rid of checkerboard patterns inside the design domain, there are many regulation methods such as mesh-independent filtering, slope-constraint, density-filtering, and morphology filtering. In this paper, the mesh-independent sensitivity filtering is used (see Ref. in [40] for more details). 4.1. Example 1: cantilever beam A cantilever beam with clamped ends (Fig. 9a) is considered first. A finite element model of the design domain is constructed by using 5000 (200  25) patches containing bilinear Q4 elements. The volume considered is constrained to be less than 50% of the domain volume. Fig. 9b and c, and Fig. 10 show topological layouts

obtained by the proposed ECP method, and iteration histories. With a zero force (F = 0 N), the topological layout in Fig. 9b, which is comparable to the layouts presented in [14,15] using the SIMP method, could be obtained using a constant initial guess. Results of empirical tests revealed that another local optimum design could be obtained using an unsymmetrical initial guess in Fig. 9c.1 To consider the effect of large displacements on the optimal layouts, a concentrated point load is now applied at the center of the bottom surface of the beam in the downward direction (Fig. 11a). With sufficiently large loads (F = 10 N and 20 N) that elicit nonlinear responses in the beam, the unsymmetrical designs in Fig. 11b and c, whose eigenfrequency histories are shown in Fig. 12, could be obtained. Table 1 and Fig. 13 compare the eigenfrequencies of the four solutions with respect to the applied loads. From the analysis of the numerical results presented in Fig. 13, it is observed that the optimal layouts considering the geometrical nonlinearity in Fig. 11 are optimized for the given loads. At F = 10 N, the fundamental eigenfrequency of the nonlinear design 1 Up the best of our knowledge, the benchmark problem has been solved without considering a fixed mass by an evenly distributed initial design only [10,14]. Consequently, it appears that the first eigenfrequency, whose mode is the bending mode, is maximized. Therefore, from our numerical examples, we could obtain the unsymmetrical design shown in Fig. 9c with a randomly distributed initial guess. These empirical numerical tests imply that the optimization problem referred to has many local optima. Furthermore, eigenfrequencies and associated eigenmodes of a few initial iterations influence the solutions of the optimization problem.

130

G.H. Yoon / Computers and Structures 88 (2010) 120–133

Table 1 Eigenfrequencies of optimized designs for the clamped beam. Design

Fig. 9b

Fig. 9c

Fig. 11b

Fig. 11c

Linear

223.06 (rad/s)

176.26 (rad/s)

256.99 (rad/s)

241.72 (rad/s)

Loading 1 (F = 10 N)

225.86 (rad/s)

189.92 (rad/s)

304.33 (rad/s)

297.24 (rad/s)

Loading 2 (F = 20 N)

222.63 (rad/s)

196.57 (rad/s)

322.41 (rad/s)

343.16 (rad/s)

Fig. 13. Load and eigenfrequency curves for the optimized results for the clamped beam.

Fig. 14. Numerical tests for examining the effects of link stiffness values. (The values in parentheses angular speed.)

of Fig. 11b, optimized for F = 10 N, are better than that of the design of Fig. 11c, optimized for F = 20 N, whereas the eigenfrequency of Fig. 11c is better than that of Fig. 11b at F = 20 N. Furthermore, a larger force (F = 20 N) produces a straight beam whose length is longer than the length of a straight beam of Fig. 11b. This example shows that a physically sound structure can be obtained using the present ECP method. Using Fig. 14 that shows numerical test results with different upper and lower bounds for the link stiffness values, similar local optimized results can be obtained in the case of both sufficiently large and sufficiently small link stiffness values (see [35]).

equivalent to those layouts obtained by solving the static compliance minimization problem (see Ref. [19] for more details). This might be explained by the fundamental structural eigenfrequency being proportional to the square root of the global stiffness. Fig. 16d plots the first frequencies against the applied loads for each design. After exceeding near F = 0.5 N, the lower beam in

4.2. Example 2: beam structure with a point load or bending load In the second example, a beam with clamped boundary conditions is considered (Fig. 15a). The design domain is discretized by 200  30 patches. To avoid obtaining a trivial void structure, the densities of two elements at the loading point are fixed as 20 kg/m3. The volume is also constrained to be less than 50% of the design domain. Fig. 15b shows an obtained linear design with a uniform initial density distribution for F = 0 N. It is found that this result is similar to the results obtained using the SIMP method too [14]. Subsequently, optimization problems with different loads were considered (Fig. 16). The layouts of the results maximizing the first eigenfrequency in Fig. 16 are, in a certain sense,

Fig. 15. Beam structure with a point load. (a) Problem definition, and (b) result 1 (a linear optimized result).

131

G.H. Yoon / Computers and Structures 88 (2010) 120–133

Fig. 16. Optimized results for point loading. (a) Problem definition, (b) result 2 (a result with a point load of 1.5 N), (c) result 3 (a result with a point load of 5 N), and (d) the curves of eigenfrequencies with respect to the applied point load.

Fig. 15b exhibits the buckling phenomenon. Consequently, the fundamental structural eigenfrequency of Fig. 15b decreases. In contrast, because the nonlinear designs of Fig. 16b and c can support larger forces without bucklings, their frequencies do not decrease at the given loads. It was observed that small density perturbations in the linear design of Fig. 15b can cause bucklings of the internal

Fig. 17. Optimized results for the bending loading case. (a) Problem definition, (b) result 2 (result for a bending moment of 1.25 Nm), (c) result 3 (a result for a bending moment of 2.5 Nm), and (d) the curves of eigenfrequencies with respect to the applied bending moment.

bar structure. Therefore, investigations on methods to control structural behaviors after post-buckling and eigenfrequency should be carried out in future. Table 2 compares the eigenfrequencies of the present designs.

Table 2 Eigenfrequencies of optimized designs for a cantilever beam subject to a point load. Design

Fig. 15b

Fig. 16b

Linear

8.43 (rad/s)

8.25 (rad/s)

Fig. 16c 7.74 (rad/s)

Shear (F = 1.5 N)

N/A (buckling of solid elements)

8.12 (rad/s)

7.77 (rad/s)

Shear (F = 5 N)

N/A (buckling of solid elements)

3.83 (rad/s)

7.87 (rad/s)

132

G.H. Yoon / Computers and Structures 88 (2010) 120–133

Table 3 Eigenfrequencies of optimized designs for a cantilever beam subject to a bending moment. Design

Fig. 15b

Fig. 17b

Fig. 17c 6.81 (rad/s)

Linear

8.43 (rad/s)

6.61 (rad/s)

Bending (BM = 1.25 Nm)

N/A (buckling of solid elements)

8.51 (rad/s)

8.28 (rad/s)

Bending (BM = 2.5 Nm)

N/A (buckling of solid elements)

8.46 (rad/s)

8.99 (rad/s)

In addition to the point load, changes in the layout designs with the bending moments were investigated (Fig. 17). Similar to the case of the point loads, a large bending moment makes the leftend-lower beam prone to buckling as shown in Fig. 17d. Consequently, the first structural eigenfrequency of the linear design of Fig. 15b is decreasing with respect to the magnitude of bending moment. Considering the geometrical nonlinearity, different topological layouts, which make the left-end parts stiffer and can support given larger moments, can be obtained in Fig. 17. The eigenfrequencies for these designs are given in Fig. 17d and Table 3. An investigation of the deformed shapes indicates that the curled left-end-upper beam of Fig. 17 becomes straight for the given bending moment, thereby increasing the overall stiffness. This example also shows that the present I-ECP method can be used to design for nonlinear dynamic structures.

5. Conclusions This paper pertains to topology optimization using the internal element connectivity parameterization (I-ECP) method in order to maximize the first eigenfrequency of a geometrically nonlinear structure. To calculate the eigenfrequencies of a nonlinear structure, modal analysis should be conducted after nonlinear static analysis using the Newton–Raphson method. However, in the framework of the standard density-based method, the modal and static analysis of nonlinear structures are hindered due to the instability of flipped elements among weak elements and the presence of locally vibrating modes. To the best of our knowledge, owing to these difficulties, there has not been any research on topology optimization in the context of nonlinear dynamic structures. To overcome the shortcomings of density-based methods, this paper proposes the use of the ECP method. This method was developed for compliance minimization problems involving large displacements, yet it was not applied to dynamic problems because of the ambiguity in the mass model. To examine whether distributed mass models can be constructed using the ECP method, this paper investigates the construction of two computational mass models. One mass model relates the mass matrix, formulated by the finite element method, to the degrees of freedom of a solid element inside a patch, while the other relates the same mass matrix to the degrees of the freedom of outer nodes forming a patch. Showing good agreement with each other for a solid element with a high stiffness value for links, the first mass model with a small stiffness value for links inherently has highly localized vibrating modes inside patches which can be troublesome on topology optimization. Therefore, this paper employs the second mass model (the patch mass modeling), which is relatively free from local modes, in preference to the direct method shown in an example with one element (Fig. 7). To illustrate the potential of the proposed approach, two numerical examples with sufficiently larger loads that elicit nonlinear responses from the structure considered are presented. For the sake of simplicity in optimization and sensitivity analysis, only the first eigenfrequency is considered. The optimized results show that it is possible to solve topology optimization problems by con-

sidering a geometrically nonlinear structure and using the I-ECP method. Moreover, we verified that results similar to those obtained in the case of the compliance minimization problem can be obtained by maximizing the fundamental structural eigenfrequency. It is suggested that in future studies, the proposed approach can be extended to layout designs by considering the post buckling phenomena, the crashworthiness, and other dynamic aspects along with material nonlinearities. A study should also be performed to compare the proposed ECP method and the discontinuous Galerkin method may be required. Acknowledgement This research was supported by the Grant of the Korean Ministry of Education, Science and Technology – The Regional Core Research Program and by the Kyungpook National University Research Fund, 2008. References [1] Bendsøe MP, Sigmund O. Topology optimization theory, methods and applications. Springer; 2003. [2] Yoon GH, Kim YY. The element connectivity parameterization formulation for the topology design optimization of multiphysics systems. Int J Numer Methods Eng 2005;64:1649–77. [3] Eschenauer HA, Lund E, Olhoff N. Topology optimization of continuum structures: a review. Appl Mech Rev 2001;54(4):331–91. [4] Yoon GH, Jensen JS, Sigmund O. Topology optimization of acoustic-structure interaction problems using a mixed finite element formulation. Int J Numer Methods Eng 2007;70:1049–76. [5] Lee JW, Wan SM, Altay D. Topology optimization for the radiation and scattering of sound from thin-body using genetic algorithms. J Sound Vib 2004;276:899–918. [6] Jog CS. Topology design of structures subjected to periodic loading. J Sound Vib 2002;253(3):687–709. [7] Diaz AR, Kikuchi N. Solutions to shape and topology eigenvalue optimization using a homogenization method. Int J Numer Methods Eng 1992;35:1487–502. [8] Yoon GH, Kim YY. Optimal design of the optical pickup suspension plates using topology optimization. Am Inst Aeronautics Astronaut 2003;41(9):1841–3. [9] Keong L, Hejun D. Topology optimization of head suspension assemblies using modal participation factor for mode tracking. Microsyst Technol 2005;11:1243–51. [10] Ma ZD, Kikuchi N, Hagiwara I. Structural topology and shape optimization for a frequency response problem. Comput Mech 1993;13(3):157–74. [11] Ma ZD, Kikuchi N, Cheng HC, Hagiwara I. Topology optimization technique for free vibration problems. J Appl Mech 1995;62:200–7. [12] Pedersen NL. Maximization of eigenvalues using topology optimization. Struct Multidisciplinary Optim 2000;20(1):2–11. [13] Kim TS, Kim YY. Mac-based mode-tracking in structural topology optimization. Comput Struct 2000;74:375–83. [14] Du J, Olhoff N. Topological design of freely vibrating continuum structures for maximum values of simple and multiple eigenfrequencies and frequency gaps. Struct Multidisciplinary Optim 2007;34:91–110. [15] Achtziger W, Kocvara M. Structural topology optimization with eigenvalues. Soc Ind Appl Math 2007;18(4):1129–64. [16] Maeda Y, Nishiwaki S, Izui K, Yoshimura M, Matsui K, Terada K. Structural topology optimization of vibrating structures with specified eigenfrequencies and eigenmode shapes. Int J Numer Methods Eng 2006;67:597–628. [17] Cho SH, Jung HS. Design sensitivity analysis and topology optimization of displacement-loaded non-linear structures. Comput Methods Appl Mech Eng 2003;192:2539–53. [18] Bruns TE, Tortorelli DA. An element removal and reintroduction strategy for the topology optimization of structures and compliant mechanisms. Int J Numer Methods Eng 2003;57:1413–30. [19] Yoon GH, Kim YY. Element connectivity parameterization for topology optimization of geometrically nonlinear structures. Int J Solids Struct 2005;42(7):1983–2009.

G.H. Yoon / Computers and Structures 88 (2010) 120–133 [20] Bruns TE. Topology optimization by penalty (TOP) method. Comput Methods Appl Mech Eng 2007;196:4430–43. [21] Buhl T, Pedersen CBW, Sigmund O. Stiffness design of geometrically nonlinear structures using topology optimization. Struct Multidisciplinary Optim 2000;19(2):93–104. [22] Pedersen CBW, Buhl T, Sigmund O. Topology synthesis of large displacement compliant mechanism. Int J Numer Methods Eng 2001;50:2683–705. [23] Schwarz S, Maute K, Ramm E. Topology and shape optimization for elastoplastic structural response. Comp Methods Appl Mech Eng 2001;190:2135–55. [24] Kemmler R, Schwarz S, Ramm E. Topology optimization including geometrically nonlinear response. In: Proceedings of the third world congress of structural and multidisciplinary optimization, Buffalo, USA; 1999. [25] Bruyneel M, Duysinx P. Note on topology optimization of continuum structures including self-weight. Struct Multidisciplinary Optim 2005;29(4):245–56. [26] Neves MM, Sigmund O, Bendsøe MP. Topology optimization of periodic microstructures with a penalization of highly localized buckling modes. Int J Numer Methods Eng 2002;54(6):809–34. [27] Neves MM, Rodrigues H, Guedes JM. Generalized topology design of structures with a buckling load criterion. Struct Optim 1995;10:71–8. [28] Rodrigues HC, Guedes JM, Bendsøe MP. Necessary conditions for optimal design of structures with a non-smooth eigenvalue based criterion. Struct Optim 1995;9:52–6. [29] Mateus HC, Rodrigues HC, Mota Soares CM, Mota Soares CA. Sensitivity analysis and optimization of thin laminated structures with a nonsmooth eigenvalue based criterion. Struct Multidisciplinary Optim 1997;14(4): 219–24.

133

[30] Bathe KJ. Finite element procedures. New Jersey: Prentice hall; 1996. [31] Cook RD, Malkus DS, Plesha ME, Witt RJ. Concepts and applications of finite element analysis. 4th ed. USA: John wiley & sons; 2001. [32] Yoon GH, Joung YS, Kim YY. Optimal layout design for three dimensional geometrical nonlinear structures using the element connectivity parameterization. Int J Numer Methods Eng 2007;69:1278–304. [33] Yoon GH, Kim YY. Topology optimization of material-nonlinear continuum structures by the element connectivity parameterization. Int J Numer Methods Eng 2007;69:2196–218. [34] Langelaar M, Yoon GH, Kim YY, Keulen FV. Topology optimization of shape memory alloy actuators using element connectivity parameterization. In: Proceedings of the sixth world congress of WCSMO, Seoul; 2005. [35] Yoon GH, Kim YY, Langelaar M, Keulen FV. Theoretical aspects of the internal element connectivity parameterization approach for topology optimization. Int J Numer Methods Eng 2002;76:77–797. [36] Chiu JT, Li YY. Modal analysis of multi-layer structure for chemical mechanical polishing process. Int J Adv Manuf Technol 2008;37:83–91. [37] Xue Y, Jairazbhoy VA, Niu X, Qu J. Large deflection of thin-plates under certain mixed boundary conditions-cylindrical bending. J Electron Packaging 2003;125:53–8. [38] Seyranian AP, Lund E, Olhoff N. Multiple eigenvalues in structural optimization problems. Struct Multidisciplinary Optim 1994;8:207–27. [39] Svanberg K. The method of moving asymptotes – a new method for structural optimization. Int J Numer Methods Eng 1987;24:359–73. [40] Sigmund O. Morphology-based black and white filters for topology optimization. Struct Multidisciplinary Optim 2007;33(4–5):401–24.