Stiffness design of heterogeneous periodic beam by topology optimization with integration of commercial software

Computers and Structures 172 (2016) 71–80

Contents lists available at ScienceDirect

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

Stiffness design of heterogeneous periodic beam by topology optimization with integration of commercial software Sinan Yi a,b, Gengdong Cheng a,⇑, Liang Xu a a

Department of Engineering Mechanics, Dalian University of Technology, Dalian 116024, PR China Beijing Institute of Mechanical and Electrical Engineering, Beijing 100074, PR China b

a r t i c l e

i n f o

Article history: Received 27 May 2015 Accepted 9 May 2016 Available online 31 May 2016 Keywords: Periodic beam Topology optimization Sensitivity analysis Homogenization method Commercial software

a b s t r a c t A topology optimization method is developed for microstructure design of heterogeneous periodic beam structure aiming at its extreme or specified effective stiffness. The effective stiffness is calculated using a FEM formulation of asymptotic homogenization method for heterogeneous periodic beam. Sensitivity of stiffness to the density design variable is derived analytically based on the solution of unit cell problems under corresponding generalized strain fields. Implementation of optimization procedure is generalized to take full advantage of commercial FEM software capabilities, with several examples presented to demonstrate its effectiveness. It is shown here the proposed method is extendible to periodic truss beam design. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Recent development of manufacturing technology make it possible to realize many different variety of engineering structure configurations to meet the requirement of multifunctional, lightweight, and efficient performance. One important category is heterogeneous periodic continuum: structure consisting of unit cells arranged periodically in one, two or three dimensions (Fig. 1). Among them, slender periodic heterogeneous beams such as sandwich beams, ribbed boxes and stranded ropes are widely applied in architectural, aeronautical and automotive industries. Structural analysis and optimum design of slender periodic heterogeneous beams become important research topics. For example, several recent research papers address structural analysis and optimization of composite corrugated skins, metallic corrugated core sandwich panels, beams with periodically variable cross-sections and graded corrugated truss core composite sandwich beams [1–4]. Periodic continuum with unit cell arranged in two or three dimensions behaves generally as a solid elastic continuum with effective material properties. The configuration and composition of their unit cell microstructure determine the overall performance of the macro structure, and the microstructure design of the unit cell is also known as material design. Sigmund [5] introduced the inverse homogenization method into material design. Combining 2D or 3D AH (Asymptotic Homogenization) theory, extreme or ⇑ Corresponding author. Tel.: +86 041184706599. E-mail address: [email protected] (G. Cheng). http://dx.doi.org/10.1016/j.compstruc.2016.05.012 0045-7949/Ó 2016 Elsevier Ltd. All rights reserved.

unique material properties can be obtained through topology optimization of the unit cell, e.g. material with negative Poisson’s ratio [5]. Many researches followed subsequently with the objective function of mechanical or multifunctional performance. For example, Guest et al. [6,7] and Xu and Cheng [8] obtained material microstructure considering both stiffness and isotropic permeability coefficient by the inverse homogenization method. Torquato et al. [9] combined the thermal and electrical conductivities with equal weight for the material design problem. However, all the above work assumes that the continuum has periodicity in three dimensions or within certain plane, where well developed AH (Asymptotic homogenization) method applies. There were several earlier investigations on the topology optimization of beam structures. Kim and Kim [10] formulated a section topology optimization technique to find the optimal cross section configuration with the objective function of a weighted sum of bending and torsional rigidities. Liu and An [11] established a topology optimization formulation based on an anisotropic beam theory considering warping of sections and coupling among deformations. Blasques and Stolpe [12] presented a novel framework for simultaneous optimization of topology and laminate properties in structural design of laminated composite beam cross sections. However, the above beam optimizations are limited to the cross section design, which assumes that the cross section along the beam axis keeps constant. If cross section of the heterogeneous beam structure varies periodically in its axial direction, and its out-axial effective properties is to be optimized, several difficulties need to be addressed in a general method.

72

S. Yi et al. / Computers and Structures 172 (2016) 71–80

Solid

Plate

Beam

Fig. 1. Various periodic structures and their homogeneous models.

The first difficulty is model reduction and effective stiffness prediction of periodic heterogeneous, beam-like structure. The cross sectional dimensions of periodic heterogeneous beam are significantly smaller than their length along the axial direction, the conventional numerical methods to analyze and optimize the overall behavior of these structures lead to heavy computations. To reduce the computational effort, this kind of periodic beam structures needs to be simplified as Euler–Bernoulii beam or Timoshenko beam by dimension reduction in addition to homogenization. Representative volume element (RVE) method and asymptotic homogenization (AH) method are two popular methods to predict effective properties of periodic structures. The RVE method has clear physical or mechanical meaning but no rigorous mathematical foundation, so it provides approximate solutions and is easy to implement. The AH method has rigorous mathematical foundation and has been successfully used in predicting effective modulus of 3D and in-plane effective modulus of 2D periodic materials both analytically and numerically [13–15]. Moreover, the sensitivity information can be obtained from the AH method, which is the foundation of the gradient-based optimization algorithm. However, FEM implementation of AH for periodic beam and plate, which does not have periodicity in its cross section or its thickness direction, has been a challenging subject. Cheng et al. [16] developed a New Implementation of the Asymptotic Homogenization (NIAH) method to predict effective properties of periodic materials with periodicity in 3D or 2D in plane, and this new FEM implementation has been extended to the homogenization method for periodic heterogeneous plate and shell structures by Cai et al. [17]. The new implementation has a rigorous mathematical foundation of the AH method, and can be simply implemented in commercial software as a black box. Based on the AH theory of heterogeneous periodic beams developed by Kolpakov and Kalamkarov [18–20], Yi et al. [21] established a FEM formulation for the AH theory of periodic heterogeneous beam structures in the same framework of the NIAH method. With this FEM implementation of AH method of periodic beam structures, one can obtain the effective stiffness of beam with complicated microstructure in its 3D unit cell with reduced computational cost. The topology optimization of the unit cell microstructure is another challenge because structural analysis of three dimensional elasticity problem is computationally intensive. Structural optimization algorithms are generally classified into gradient-based algorithm and population-based algorithm [22–24]. The population-based algorithm simulates the natural phenomenon through random processes to search the global optimal solution in the design space and needs large number of function

evaluations. The gradient-based algorithm searches the optimal solution in special gradient direction, thus it has high efficiency with fewer iterations and faster convergence. For topology optimization of microstructure design with thousands or even more design variables, each iteration requires three dimensional FEM analysis, thus the gradient-based algorithm is more suitable for microstructure design problems. For gradient-based optimization, the availability of analytical sensitivity analysis has significant influence to the optimization efficiency. The sensitivity-based optimization algorithm can actually take longer than the zero-order optimization algorithm if the sensitivity calculation is inefficient, for example, with the finite difference method. Herein we show with the NIAH implementation recently proposed by authors, analytical sensitivity analysis is possible to enable gradient-based optimization algorithm for periodic heterogeneous beam design. Such topology optimization can be even more efficient by proper integration with commercial FEM software such as ANSYS, ABAQUS, to leverage their pre-processing and post-processing ability, expanded element types and efficient linear or eigenvalue matrix solver. To integrate the newly developed structural optimization algorithm with the existing commercial software, one key problem is sensitivity calculation, which most commercial software lack. Such demand drove the development of semianalytic sensitivity method [25] and its many extensions. Zhang et al. [26] developed a practical tool to deal with sizing sensitivity analysis with ABAQUS. Jeong et al. [27] proposed a novel approach using indirect calculation of internal finite element information for a stress-based topology optimization problem. Henrichsen et al. [28] performed optimum stiffness design of laminated composite structures using the commercially available programs ANSYS and MATLAB. This paper develops a topology optimization method for the microstructure design of periodic beam structure under the available material volume. The FEM element density of the unit cell is the topology design variable. The objective is to design extreme or specified effective stiffness using NIAH method efficiently. Our analytically derived sensitivities of the effective stiffness is proportional to the element strain energy or mutual strain energy corresponding to the generalized unit strain fields, extracted using commercial software and indirect calculation of internal finite element information. Section 2 focuses on the mathematical formulation of microstructure design problem for periodic heterogeneous beam structure. In section 3, we introduce the FEM implementation of asymptotic homogenization theory for periodic beam structure. Section 4 describes in detail sensitivity analysis and the integration of evaluation of the objective function, sensitivity

73

S. Yi et al. / Computers and Structures 172 (2016) 71–80

assumed to be a function of the element density given by the SIMP interpolation scheme as:

Ee ¼ qpe E0 ;

qmin 6 qe 6 1

where E0 is the Young’s modulus of solid material, and p is the penalization power. An alternative interpolation scheme can also be used, including Modified SIMP and RAMP [30]:

Periodic heterogeneous beam with effective stiffness D H Unit cell

Ee ¼ Emin þ qpe ðE0  Emin Þ;

Modified SIMP : Ee ¼

RAMP : Fig. 2. Microstructure design of periodic heterogeneous beam structure.

calculation and the optimization procedure with commercial software. Combining efficient calculation of objective function and sensitivity and the modeling flexibility in commercial software, this optimization approach has significantly lower computational cost compared to traditional approach and is extendable to periodic truss beam structures (bar cross-section optimization). Section 5 demonstrates the optimization method in specific examples of the optimization method and its extension to periodic truss beam structures, with future direction summarized in final section. 2. Mathematical formulation of the optimization problem

q ¼ ðq1 ; q2 ; . . . ; qN ÞT f ðDHij ðqÞÞ i; j ¼ 1; . . . ; 4 N X

qe v e =V   1 6 0

qe

1 þ pð1  qe Þ

E0 ;

0 6 qe 6 1

0 6 qe 6 1

ð3Þ ð4Þ

where Emin is a small number representing the stiffness of void material, which is used to avoid singularity of the stiffness matrix and is usually set Emin ¼ 1e  9. In the formulation (1), DH is the 4  4 matrix which is defined as the effective stiffness of beam structure in the macro scale between the generalized axial force F, bending moments My, Mz in two outaxial directions and twisting moment T and the generalized axial extension ex , bending curvatures jy, jz in two out-axial directions, torsion jyz, i.e.,

3 2 D11 F 7 6 6 6 M y 7 6 D21 7¼6 6 7 6 6 4 M z 5 4 D31 2

T

The optimization problem considered in the present work is to find optimal topology of periodic beam structure with extreme or specified effective stiffness in macro by optimizing its micro structural topology of unit cell under available material volume constraint. Since the effective stiffness of periodic beam structure is obtained by analysis of its unit cell, which will be shown in Section 3, the optimization problem is formulated on the unit cell scale as follows. By using NIAH method, the unit cell can be modeled by any type of elements and modeling techniques available in commercial FEM software. However, for sake of clarity, hereafter we use 3D solid element to model the unit cell and assume the design domain of the unit cell is discretized into N finite elements, see Fig. 2. The artificial material density qe ¼ 0 or 1; e ¼ 1; . . . ; N of each element is design variable and represents a void or solid in the unit cell, respectively. In this way, the material density distribution of qe ; e ¼ 1; . . . ; N decides the micro-structural topology of the unit cell. In order to make the solution process computationally tractable, the so called solid isotropic material with penalization (SIMP) model proposed in [29] is applied. The optimization problem described above can be formulated as an optimal material distribution problem.

8 > find > > > > > min > > > > < s:t: > > > > > > > > > > :

ð2Þ

D41

D14

32

ex 3

D12

D13

D22

D23

D32

D33

76 D24 76 jy 76 76 D34 54 jz

D42

D43

D44

7 7 7 7 5

ð5Þ

jyz

DH can be calculated based on the solution of the four unit cell 3D elasticity problems, which will be detailed in Section 3. The material volume constraint V  is imposed on the material of the unit cell, see ((1)b). In the formulation ((1)c), Dlower is the minimum value of the main stiffness to avoid the optimized result singularity, which is considered only when the structural topology has obvious separation. The objective function f in ((1)a) is a function of the effective stiffness DH , and its definition varies with different problems. For example, for the maximum torsion stiffness problem in the example 5.2, f is defined as:

f ¼ DH44 ðqe Þ=Dconst

ð6Þ

where Dconst is the normalized constant, which can be specified or chosen as the torsion stiffness of the initial design to make the objective function f non-dimensional. In addition, for specified stiffness design, the objective function f can be defined as the regularized 2-norm form:

f ¼

X

wij

DSij –0

DHij DSij

!2 1

þ

X

2

wij ðDHij  DSij Þ

ð7Þ

DSij ¼0

or the D-function form [31]:

ðaÞ ðbÞ

ð1Þ

f ¼

X

  X  2 wij DHij lnðDHij =DSij Þ  DHij þ DSij þ wij DHij  DSij

DSij –0

ð8Þ

DSij ¼0

e¼1

; i ¼ 1; . . . ; 4 DHii ðqÞ P Dlower ii 0 < qmin 6 qe 6 1; e ¼ 1; 2; . . . ; N

ðcÞ ðdÞ

where qe ; e ¼ 1; . . . ; N is the element artificial density, ve is the element volume. The subscript e refers to the eth element. V⁄ is the given amount of material. N is the total number of elements in the design domain. qmin is the lower bound of element density variable to avoid the possible singularity of structural stiffness matrix. Note that the design variable qe ¼ 0 or 1; e ¼ 1; . . . ; N is relaxed in ((1)d) in order to avoid solution of 0–1 programming. To obtain 0–1 design, the artificial Young’s modulus of eth element Ee is

where DSij is the specified target stiffness, and wij is the weight coefficient. The second term in (7) and (8) suppresses the coupling terms, or ensures a symmetric microstructure when the unit cell is symmetric one. For specified stiffness design, the volume constraint in ((1)b) may not be active if the target stiffness is chosen too small. On the other hand, the optimization problem (1) may not have a feasible design when the target stiffness is too large in comparison with the given material volume. As a result, iteration algorithm for formulation (1) often becomes difficult to converge. An alternative formulation for specified stiffness design is as follows:

74

S. Yi et al. / Computers and Structures 172 (2016) 71–80

: q ¼ ðq1 ; q2 ; . . . ; qN ÞT N X min : f ðqÞ ¼ V10 qe v e find

ðaÞ

e¼1

s:t:

: 1  DHij =DSij 6 0 i; j ¼ 1; . . . ; 4

ðbÞ

: DHii ðqÞ P Dlower ; ii

ðcÞ

i ¼ 1; . . . ; 4

: 0 < qmin 6 qe 6 1;

e ¼ 1; 2; . . . ; N

ð9Þ

ðdÞ

Here V 0 denotes the material volume when all density design variables are 1. In the new formulation, we design the material distribution of the microstructure under specified stiffness constraint for minimum material volume. 3. AH theory for periodic beam structures and implementation by NIAH method In this section, as the base of microstructure design optimization, the FEM implementation of the asymptotic homogenization method for periodic beam structures are introduced briefly, including the AH theory, FEM formulation and implementation by NIAH method. 3.1. Asymptotic homogenization theory for heterogeneous periodic beam structures Here the basic assumption of the theory and the unit cell problems are shown briefly. For detailed derivation of the theory, the interested readers may refer to Kolpakov and Kalamkarov [18– 20] and Yi et al. [21]. The theory concerns the homogenization of structures having one large global dimension in comparison with the others. As illustrated in Fig. 3, a domain of periodic structure obtained by repeating a certain small periodic cell eY along the Ox1 axis is considered. Here, the parameter e is the characteristic dimension of the periodic cell, assumed small and approaching zero ðe ! 0Þ. The overall beam structure domain in macroscopic coordinates is defined as Xe with the surface boundary @ Xe ¼ Se [ Su , while the definition of the unit cell domain in microscopic coordinates is Y with periodic boundary x and non-periodic boundary S. The asymptotic homogenization method starts from perturbing the structural displacement and stress to the asymptotic series in terms of the small parameter e, substituting them into the three dimensional equilibrium equations, and grouping the terms multiplied by the same power of e. Through the above perturbation procedure, the characteristic displacement and stress fields are obtained from solving four characteristic equations of extension, flexure and torsion in the unit cell domain, which is called the unit cell problems and can be formulated as unified form:

8 @bij > ¼0 > @yj > > > < bij Nj ¼ 0   > bij Nj  ¼ bij Nj  > > xþ x > > :u j ¼ u j i xþ

in Y on S

ð10Þ

on x on x

i x

Here, the repeated indices denote the summation convention, Nj denotes jth component of unit outward normal on boundary S and ui denotes ith component of displacement field. The characteristic stress fields bij have specific expressions for different characteristic equations:

@U 11 m þ C ij11 @yl

ð11Þ

@V am1  ya C ij11 @yl

ð12Þ

Extension bij ¼ C ijml a

Flexure bij ¼ C ijml 

Torsion bij ¼ C ijml

@W 11 m @yl

þ eab yb C ija1

ð13Þ

where the Latin indices range from 1 to 3, and the Greek ones range from 2 to 3. U 11 ; V a1 ; W 11 are the characteristic displacement fields substituted for u, and the symbol eab is defined as e23 ¼ 1; e32 ¼ 1; e22 ¼ e33 ¼ 0. Neglecting the influence of transverse shear and warping, the periodic beam structure is equivalent to an Euler–Bernoulli beam model in macro with the effective stiffness matrix as: 2

6 6 D21 D22 D23 6 D ¼6 6 D31 D32 D33 4 H

D41 D42 D43

3

2

hb11 i 7 6 6 hy2 b11 i D24 7 7 6 7¼6 6 D34 7 hy3 b11 i 5 6 4 D44 helk yk bl1 i

D11 D12 D13 D14

2

3

hb11 i

3



hb11 i

hb11 i

7  hy2 b11 i 7 7 7 7 2 3  hy3 b11 i hy3 b11 i hy3 b11 i 7 5 2 3  helk yk bl1 i helk yk bl1 i helk yk bl1 i ð14Þ 2

hy2 b11 i

3

hy2 b11 i

Here D11 is the stretching stiffness in the periodic direction, D22 ; D33 the bending stiffness, D44 the torsional stiffness and the off-diagonal terms are corresponding coupling stiffness. The angular bracket denotes averaging over the volume Y of the unit cell:

hwi ¼

1 jYj

Z

Y

wdy1 dy2 dy3

3.2. FEM formulation and implementation in commercial software The asymptotic homogenization theory of periodic beam structures seems difficult to implement numerically. Note that the form of Eqs. (10)–(13) is similar to the equilibrium equations and constitutive equations in the elasticity theory. According to the principle of virtual work, for arbitrary general virtual displacement v, the equivalent integral form of the Eq. (10) is expressed as:

Z

Y

v i bij;j dY 

Z

S

v i bij Nj dS 

Z

xþ

v i bij Nj dS 

Z

x

v i bij Nj dS ¼ 0

ð15Þ

The last two terms on opposite boundary x cancel each other out due to the periodicity of the boundary displacements and stresses, and then making integration by parts within Y domain, the equivalent integral weak form is obtained as:

Z

Y

v i;j bij dY ¼ 0

ð16Þ

Discretizing the unit cell and applying the variation procedure of elasticity in finite element method, the following finite element formulation of Eqs. (11)–(13) and (16) with the periodic displacement boundary condition can be drawn:

Kvi ¼ Fi Fig. 3. The heterogeneous beam structure and its unit cell.

ð17Þ

75

S. Yi et al. / Computers and Structures 172 (2016) 71–80

where

Z

K¼

Z

BT CBdY

ð18Þ

BT Cei0 dY

ð19Þ

v0node

Y

Fi ¼ Y

where K is the global stiffness matrix of the unit cell, B is the finite element strain–displacement matrix and C is the constitutive stiffness matrix of material. Fi are nodal force vectors corresponding to the generalized unit strain fields ei0 , which are defined as:

e¼

8 e11 9 > > > > > > > > > > > > > e 22 > > > > > > > > > > > > = 33

e10 ¼

> > > c12 > > > > > > > > > > > > > > > c > > 23 > > > > > > : ;

c31

e ¼ 4 0

8 0 > > > > > > > 0 > > > > > > <0

8 9 1> > > > > > > > > > > > > 0> > > > > > > > > > > > <0> = > > > 0> > > > > > > > > > > > > > > 0 > > > > > > > : > ; 0

e20 ¼

8 9 y2 > > > > > > > > > > > > > 0 > > > > > > > > > > > > <0 > = > > 0 > > > > > > > 0 > > > > : 0

> > > > > > > > > > > > > ;

e30 ¼

8 9 y3 > > > > > > > > > > > > > 0 > > > > > > > > > > > > <0 > = > > 0 > > > > > > > 0 > > > > : 0

9 > > > > > > > > > > > > > =

> > > > > > > > > > > > > ; ð20Þ

Y

To avoid assembling or extracting the global stiffness matrix K, we can further apply the characteristic displacement field v j on every node in commercial software, run one static analysis, and

Dij ¼

i T

j 0

j

i; j ¼ 1; . . . ; 4

ð21Þ

The Eqs. (17) and (21) are governing equations, which relate the microstructure design variables of the unit cell to the effective stiffness of equivalent Euler–Bernoulli beam in macro. The above finite element formulation remains difficult to solve because the force vector and strain energy in Eqs. (19) and (21) are formulated by the element-related matrix, such as constitutive matrix and strain–displacement matrix, which compels us to find the corresponding formulations and write the code for different element types. One novelty of the NIAH method is to obtain the force vector and strain energy by the displacement fields that corresponds to the generalized unit strain and thus be able to implement the FEM formulation by using the commercial software as a black box. For this purpose, the generalized unit strain fields e0 are represented by the corresponding displacement fields v0 ,

e0 ¼ Bv0

ð22Þ

and rewrite the forceZvector in (19) as: Z

Fi ¼

ð25Þ

ð26Þ

So Eq. (25) can be rewritten as:

ðe  e Þ Cðe  e ÞdY; i 0

1 i T ðv  vi Þ Kðv0j  v j Þ i; j ¼ 1; . . . ; 4 jYj 0

P j ¼ Kv j

Then the effective stiffness matrix DHij can be obtained as follows:

Z

ð24Þ

get the corresponding nodal reaction force P j , which means:

> > > y3 > > > > > > > > > > > > > > > 0 > > > > > > > > : ; y2

1 ¼ jYj

9 > > > =

For beam, plate and shell elements, which can be used to model the unit cell if there is need, the corresponding nodal displacement fields can be found in [21]. After getting the characteristic nodal displacement fields vi by solving Eq. (17), we rewrite the stiffness matrix (21) in the mutual strain energy form calculated by the initial nodal displacement fields vi0 and the characteristic nodal displacement fields vi :

Dij ¼

In practical numerical calculation in commercial software, the periodic boundary condition can be applied by coupling the relevant nodes in the opposite faces of the unit cell in the periodic direction.

DHij

v30node

8 9 8 9 x> xy > > > > > > > > > > > < = < = ¼ v v10node ¼ 0 v20node ¼ x2 =2 > > > > > > > > > > > > > > > > : > : : > ; ; ; 0 0 w 8 8 9 9 xz > 0 > > > > > > > > > > > < < = = 4 ¼ 0 v0node ¼ xz > > > > > > > > > > > : : 2 > ; ; x =2 xy 8 u > > > <

BT CBvi0 dX ¼

BT CBdXvi0 ¼ Kvi0

ð23Þ

In commercial software, we can apply nodal displacement vector on the unit cell and calculate nodal reaction force from static analysis without dependence of element type and modeling technique, as long as the stiffness matrix K is available. For solid element, the nodal displacement fields v0 corresponding to the generalized unit strain field e0 in Eq. (20) can be constructed as:

1 i T j ðv  vi Þ ðf  P j Þ jYj 0

ð27Þ

Note that all quantities in Eq. (27) are either given or outputs directly from the commercial software, thus it is easy to calculate effective stiffness of the periodic beam structure by using the commercial software as a black box. 4. Sensitivity analysis and implementation by commercial software In order to solve the optimization problem using gradient-based optimization algorithm, the sensitivity analysis of the objective function needs to be further considered. The sensitivity of the effective stiffness with respect to the design variable qe is derived analytically as:

@DHij 1 ¼ @ qe jYj

Z

T

ðei0  ei Þ

Y

1  jYj

Z  Y

@ ei @ qe

@C j ðe  e j ÞdY @ qe 0

T

Cðe0j  e j ÞdY 

1 jYj

Z Y

T

ðei0  ei Þ C



 @e j dY ð28Þ @ qe

Note that the equivalent integral weak form Eq. (16) can be expressed in the FE-notation as:

1 jYj

Z

Y

eðv ÞT Cðei0  ei ÞdY ¼ 0

ð29Þ

where eðv Þ is the arbitrary strain field within the unit cell domain, thus the last two terms in Eq. (28) are zero. Since the i-th design variable qi is defined in the i-th element, the sensitivity formulation can be transformed from the unit cell domain Y to the element domain Ye:

@DHij 1 ¼ @ qe jYj

Z

1 ¼ jY e j

Y

T

ðei0  ei Þ

Z Ye

@C j ðe  e j ÞdY @ qe 0 T @Ce ðEe Þ

ðei0  ei Þ

@ qe

ðe0j  e j ÞdY e

ð30Þ

76

S. Yi et al. / Computers and Structures 172 (2016) 71–80

Here Ce is the elastic constitutive matrix of the element e. If the interpolation scheme is utilized, the sensitivity formulation can be further derived as:

@DHij ¼ aW ije @ qe

ð31Þ

where

Z W ije ¼

Ye

1 i T ðe  ei Þ Ce ðe0j  e j ÞdY e ; 2 0

i; j ¼ 1; . . . ; 4

ð32Þ

When the interpolation schemes in (2)–(4) are used, the parameter

a can be written as: 2n a¼ qe jYj

SIMP :

ð33Þ

2nqn1 ðE0  Emin Þ e Ee jYj 2ð1 þ pÞ a¼ qe ½1 þ pð1  qe ÞjYj

Modified SIMP : RAMP :

a¼

ð34Þ ð35Þ

We can see from Eqs. (31)–(35), no matter what interpolation scheme is used, the sensitivity can always be written as the product of W ije and the corresponding parametric term a. The parametric term can be easily calculated, so the next problem is how to obtain W ije . For i ¼ j, W iie denotes the strain energy of the element e corresponding to the generalized unit strain case of ðei0  ei Þ, which can be extracted directly from the solution data by commercial software. For i–j, W ije is the element mutual strain energy of the general-

ized unit strain cases ðei0  ei Þ and ðe0j  e j Þ, and cannot be obtained directly from commercial software, so another approach need to be find to calculate it. Let us define a new generalized strain case eiþj as the sum of the unit strain cases ei and e j :

eiþj ¼ ei þ e j

ð36Þ

where

ei ¼ ei0  ei The corresponding element strain energy of the strain case be presented as:

ð37Þ

e

iþj

can

Z

Z 1 iþj T 1 i T ðe Þ Ce ðeiþj ÞdY e ¼ ðe þ e j Þ Ce ðei þ e j ÞdY e 2 2 Ye Ye  Z  1 i T i 1 j T T ðe Þ Ce e þ ðe Þ Ce e j þ ðei Þ Ce e j dY e ¼ 2 2 Ye

Periodic heterogeneous beam-like structures, such as composite corrugated skins, metallic corrugated core sandwich panels, beams with periodically variable cross-sections and graded corrugated truss core composite sandwich beams are widely used in industries. Many high-rise building, crane arm, bridge deck slab have nearly periodic configuration along its length and are often simplified as beam in the primitive design stage. For such structures, their effective stiffness, including bending stiffness and twisting stiffness are the most concerned mechanical characteristics, other stiffnesses can also be considered. In this section, we present three numerical examples. The first one concerns bending stiffness topology optimization, in which two optimization formulations are examined and discussed, the second example deals with maximum torsional stiffness design of a beam with square holes, and the last example deals with size optimization of a truss beam. 5.1. Microstructure design with two specified bending stiffness The first example is to design the microstructure of periodic beam with variable square cross section, with two bending stiffness targets. The design domain of the unit cell consists of two alternate cubes of L1 = 20 mm, L2 = 10 mm, shown in Fig. 5, and the material property parameters are E0 = 200 GPa, m = 0.3. The modified SIMP in Eq. (3) with the penalty parameter p = 3 is applied. The prescribed material volume fraction of the design domain V  is 20%. Two cases are tested. In case I, the target bending stiffnesses are chosen as D22 = 100 N m2, D33 = 160 N m2 with the weight coefficients w22 ¼ w33 ¼ 1. In case II, D22 = 160 N m2, D33 = 200 N m2. To avoid the checker-board problem and get clear black-white topology result, the volume preserving nonlinear density filter based on Heaviside functions [34] is applied, which is

q~ ¼

g½ebð1q=gÞ  ð1  q=gÞeb  06q6g ð1  gÞ½1  ebðqgÞ=ð1gÞ þ ðq  gÞeb =ð1  gÞ þ g g < q 6 1 ð40Þ

ð38Þ

So the element mutual strain energy W ije can be calculated as follows: ii jj W ije ¼ ðW iþj e  W e  W e Þ=2

5. Numerical examples

(

W eiþj ¼

¼ W iie þ W jje þ 2W ije

avoids the extraction of element stiffness matrix. We verified equations (32)–(39) and the above algorithm by finite difference method for a number of examples. The gradient-based optimization procedure is straightforward with the objective function and sensitivity available from commercial software, as shown in Fig. 4. Here the Method of Moving Asymptotes (MMA) algorithm [32,33] is used as the optimizer, which proves to be stable and suitable for large scale topology optimization.

ð39Þ

Based on Eqs. (32)–(39), all the element strain energy and mutual strain energy can be calculated under the corresponding strain cases. Following the essence of the NIAH method, the elecaused by the strain case eiþj can also be ment strain energy W iþj e implemented by nodal displacement fields. We can just construct i þ v  j , where v  i ¼ vi0  vi , as  iþj ¼ v the nodal displacement field v input to the finite element model of the unit cell. From static stress analysis, the strain energy of every element and the corresponding mutual strain energy W ije can be calculated by Eq. (39). This analysis step is merely a matrix multiplication with vector and thus

~ are the density before and after the nonlinear denwhere q and q sity filter. The parameter g is chosen to preserve the structural material volume before and after the nonlinear density filter. Fig. 6 illustrates the iteration history of the formulation (1) for D22, D33 and material volume fraction. After early oscillation, iteration stabilizes in the later stage, with D22, D33 converging to the target stiffness approximately, but the actual amount of material is less than material constraints (20%), which means the material volume constraint is not active in this example. Optimal structure needs less material to meet the target stiffness, causing the oscillations of the volume constraint and the two stiffnesses, resulting in the difficulty of optimal convergence. Fig. 7 illustrates the iteration history using the alternative formulation (9), in which the material volume fraction is objective function, and D22, D33 are constrained. The two effective bending stiffnesses approach the target value steady in the optimization process, except for a small oscillation when the parameter b

77

S. Yi et al. / Computers and Structures 172 (2016) 71–80

Fig. 4. The flowchart of microstructures design by commercial software.

240

L1

D22

220

D

Bending Stiffness (N m2)

L2 Fig. 5. The design domain of the beam and its unit cell.

200

1

180

0.9

D

22

140

0.7

D

33

120

Vol

0.6

100

0.5

80

0.4

60

0.3

40

0.2

20

0.1

0

0

200

400

600

800

Volume Fraction

0.8

2

Bending Stiffness (N m )

160

Vol

1

180

0.9

160

0.8

140

0.7

120

0.6

100

0.5

80

0.4

60

0.3

40

0.2

20

0.1

0

Volume Fraction

33

200

0 0

100

200

300

400

500

Iteration Fig. 7. Iteration history of optimization based on new formulation (D22 = 100 N m2, D33 = 160 N m2).

0 1000

Iteration Fig. 6. The iteration history of D22, D33 and volume fraction for the case D22 = 100 N m2, D33 = 160 N m2.

changes in the nonlinear density filter. The material volume fraction as the objective function decreases gradually in the entire optimization process, and finally turns to be 16.76% at the end of the optimization, which is smaller than the value of original formulation 20%. This means that material consumption reaches a minimum by optimization technology. The optimal topology result of periodic beam by this new formulation is shown in Fig. 8. Due to the use of nonlinear density filter, clear 1/0 material density distribution is formed, and we delete the elements of zero density value in order to display the internal topology more clearly. We can see that the configurations in two directions are different. Holes appear in the horizontal direction, which is the direction of the weak bending stiffness. In the

Fig. 8. The material distribution for the case of D22 = 100 N m2, D33 = 160 N m2.

direction of the strong stiffness, the thickness of the wall on both sides appeared to be strengthened, which is physically meaningful. For case II with D22 = 160 N m2 and D33 = 200 N m2, the iteration history and topology result of the formulation (9) are shown in Fig. 9 and Fig. 10. Similar to Case I, bending stiffness constraints are satisfied accurately, and the objective function declines gradually in the iteration. The final material volume fraction is 24.11%, slightly higher than 16.76% of the previous case, because of the stronger bending stiffness goal as constraints.

78

S. Yi et al. / Computers and Structures 172 (2016) 71–80 240

1

180

0.9

160

0.8

140

0.7 D

120

22

D

100

33

Vol

80

0.6 0.5 0.4 0.3

40

0.2

20

0.1

0.8

600 D44 Vol

500

0.7

400

0.6

300

0.5

200

0.4

100

0.3

0

0 0

50

100

150

200

250

300

350

400

450

500

Iteration Fig. 9. Iteration history of optimization after specified stiffness improved (D22 = 160 N m2, D33 = 200 N m2).

Volume Fraction

2

60

Fig. 12. Topology design of periodic beam for maximum twisting stiffness.

Torsion Stiffness D44 ( N m )

200

Volume Fraction

Bending Stiffness (N m2)

220

0.2

0 0

50

100

150

200

250

300

Iteration Fig. 13. Iteration history of torsion stiffness and volume fraction.

Fig. 10. The material distribution for the case of D22 = 160 N m2, D33 = 200 N m2.

The optimization result in Fig. 10 shows that reinforced ribs appear in the direction of strong bending stiffness compared with the previous case, and in the direction of weak bending stiffness the hole gets smaller and the wall gets thicker, therefore increasing bending stiffness in both directions. This example implies that for the specified stiffness optimization, the formulation using the specified stiffness as constraints and minimizing the material volume achieves more stable optimization iterations and better results. The initial formulation using fixed material volume as constraint to optimize for specified stiffness, causes oscillation and convergence difficulty, and can lead to material waste.

5.2. Maximum torsion stiffness design This example is the extreme stiffness design for maximum torsion stiffness. The beam in macro has periodicity in x direction, and the design domain of the unit cell is a cube with a square hole in y direction, see Fig. 11. The size of the cube is 20 mm 20 mm  20 mm, and the hole is 10 mm  20 mm  10 mm. The material of the unit cell is isotropic with E0 = 200 GPa, v ¼ 0:3. The prescribed volume fraction of the design domain is 40%, and

the torsion stiffness obtained in the first iteration is chosen as the constant Dconst. The optimal topology of the maximum torsion stiffness is shown in Fig. 12. From the result, we can see that a non-trivial new force transmission path is built since the traditional trivial path that material distributed in the exterior has been blocked by the given square hole, and in the new path the material reduced in the interior along the axis and the corners of the unit cell. From the iteration history in Fig. 13, it can be seen that the torsion stiffness improves rapidly with the iteration. This conceptual topology optimization in Fig. 12 can be further processed for a practical design. 5.3. An extension to periodic truss beam Periodic truss beam structures are widely used in satellite aerial, high-rise buildings and other engineering areas. For this kind of structures, the cross section areas of bars are chosen as design variables, and the mathematic formulation of the optimal unit cell design is expressed as:

8 > > > find > > > min > > > > < s:t: > > > > > > > > > > :

: A ¼ ðA1 ; A2 ; . . . ; AN ÞT : f ðDHij ðAÞÞ i; j ¼ 1; . . . ; 4 :

N X Ae Le =V   1 6 0

ð41Þ

e¼1

: DHii ðAÞ > Dlower ii

i ¼ 1; 2; . . . ; 4

: Amin 6 Ae 6 Amax

e ¼ 1; 2; . . . ; N

where Ae and Le are the cross section area and the length of the bar e, Amin and Amax are the lower and upper bound of cross section area, V  is the prescribed material volume of the unit cell. The sensitivity of the effective stiffness DHij to the design variable Ae can be written as:

Fig. 11. The design domain of the beam and its unit cell.

@DHij ¼ aW ije ; @Ae

where a ¼

2 Ae jYj

ð42Þ

79

S. Yi et al. / Computers and Structures 172 (2016) 71–80

Fig. 14. The periodic truss beam structure and its unit cell.

Table 1 Effective stiffness results of truss beam microstructure design.

11

2

D11 (10 N m ) D44 (1011 N m2) D14 (1011 N m2)

Initial

Min D14

Max D14

D14 to 0

1.0130 1.1952 0.3586

1.1276 1.6202 1.1984

0.2544 1.5000 0.2016

0.7872 0.7352 <106

It can be seen that sensitivity information is also related to the strain energy or mutual strain energy of the bar element, and can be obtained from commercial software as shown in the flowchart of Fig. 4. Here one example of extreme stiffness design of the truss beam structure is presented. The beam has elliptical cross section, which rotates its orientation as it moves along the axis in Fig. 14. The geometry of the unit cell is L ¼ 8 m, a ¼ 2 m,b ¼ 1:5 m. The initial cross section area A0 is 0.05 m2, and the lower and upper bound Amin and Amax are 0.005 m2 and 0.5 m2. The prescribed material volume V  equals to the initial material volume 16.7014 m3 and keeps constant in the process of optimization. The material properties of the bars are E0 ¼ 200 GPa;v ¼ 0:3. Because of the particular configuration of this beam structure, for the initial structure with all cross sections of the bars being uniform, the extension-twisting coupling effect exists, and the effective stiffness of the beam calculated are shown in Table 1. By applying (41), we can design the cross sectional areas of bars in the unit cell and optimize the extension-twisting coupling

(a) Min D14

stiffness D14 to obtain maximum or minimum coupling stiffness, or even specify the coupling stiffness D14 = 0. The corresponding objective functions can be simply defined as:

Maximum D14 :

f ¼ D14

Minimum D14 : f ¼ D14 Specify D14 to 0 : f ¼ D214 For the case of specifying D14 to 0, the inequality of volume constraint is not an active constraint. In other words, we can obtain a set of designs with D14 = 0 and different material volume, their cross sectional area are proportional to each other. To make the result unique and comparable, the equality of volume constraint is imposed and the mathematic optimization formulation is modified as:

8 > find > > > > > min > > > > < s:t: > > > > > > > > > > :

: A ¼ ðA1 ; A2 ; . . . ; AN ÞT : f ðDHij ðAÞÞ i; j ¼ 1; . . . ; 4 :

N X Ae Le =V   1 ¼ 0

ð43Þ

e¼1

: DHii ðAÞ > Dlower ii : Amin 6 Ae 6 Amax

i ¼ 1; 2; . . . ; 4 e ¼ 1; 2; . . . ; N

The three optimal results for minimum D14, maximum D14 and D14 = 0 of effective stiffness by the truss microstructure design are shown in Table 1, and the truss cross section areas distributions of three cases are presented in Fig. 15. Through material redistribution of the bars, the extension-twisting coupling stiffness varies

(b) Max D14

(c) Specify D14 to 0

Fig. 15. Optimal truss beam structure for (a) minimum extension-twisting stiffness, (b) maximum extension-twisting stiffness, (c) zero extension-twisting stiffness.

80

S. Yi et al. / Computers and Structures 172 (2016) 71–80

greatly, even changes the sign of D14. The negative D14 in the 3rd row, Table 1 means the coupling rotation is opposite when you stretch this beam. Conversely, introducing the equality of volume constraint then allows a high-precision optimization of truss beam with no coupling effect. This example proves that our method can optimize the coupling stiffness, and the sensitivity analysis by commercial software is effective. 6. Conclusion This paper presents an optimization procedure for the microstructure design of periodic beam structure, which aims at its extreme stiffness or specified stiffness under available volume. The mathematical formulation of design problem is established based on design variables of element density of the unit cell, and the objective function related to the beam stiffness matrix can be calculated by the NIAH method efficiently. The sensitivities of the macro effective stiffness are obtained from the element strain energy or mutual strain energy, which can be obtained from the solution of the unit cell problem under the corresponding generalized strain fields and extracted from commercial software directly. The entire topology optimization procedure can be integrated with the commercial software as demonstrated in this paper. With efficient calculation of objective function and sensitivity and the modeling flexibility in commercial software, the efficiency of optimization is significantly improved and the overall coding and computational workloads are reduced. Two typical examples illustrate the validity of proposed optimization method, and the optimization formulations for specified stiffness design are discussed. An extension to the optimization of periodic truss beam structures is presented with different coupling stiffness achieved in the same optimization framework. Acknowledgements This research was supported by National Natural Science Foundation (11332004, 91216201) and the 973 Program (No. 2014CB049000). The financial supports are gratefully acknowledged. We are thankful for Krister Svanberg for his MMA program made freely available for research purposes. References [1] Liang CC, Yang MF, Wu PW. Optimum design of metallic corrugated core sandwich panels subjected to blast loads. Ocean Eng 2001;28:825–61. [2] Shaw AD, Dayyani I, Friswell MI. Optimisation of composite corrugated skins for buckling in morphing aircraft. Compos Struct 2015;119:227–37. [3] Zheng T, Ji T. Equivalent representations of beams with periodically variable cross-sections. Eng Struct 2011;33:706–19. [4] Xu G, Yang F, Zeng T, Cheng S, Wang Z. Bending behavior of graded corrugated truss core composite sandwich beams. Compos Struct 2016;138:342–51. [5] Sigmund O. Design of material structures using topology optimization. Lyngby: Technical University of Denmark; 1994. [6] Guest JK, Prévost JH. Optimizing multifunctional materials: design of microstructures for maximized stiffness and fluid permeability. Int J Solids Struct 2006;43:7028–47.

[7] Guest JK, Prévost JH. Design of maximum permeability material structures. Comput Methods Appl Mech Eng 2007;196:1006–17. [8] Xu S, Cheng G. Optimum material design of minimum structural compliance under seepage constraint. Struct Multidisc Optim 2010;41:575–87. [9] Torquato S, Hyun S, Donev A. Multifunctional composites: optimizing microstructures for simultaneous transport of heat and electricity. Phys Rev Lett 2002;89:1–4. [10] Kim YY, Kim TS. Topology optimization of beam cross sections. Int J Solids Struct 2000;37:477–93. [11] Liu X, Hu G. A continuum micromechanical theory of overall plasticity for particulate composites including particle size effect. Int J Plast 2005;21:777–99. [12] Blasques JP, Stolpe M. Multi-material topology optimization of laminated composite beam cross sections. Compos Struct 2012;94:3278–89. [13] Papanicolau G, Bensoussan A, Lions JL. Asymptotic analysis for periodic structures. Amsterdam: North Holland Publ.; 1978. [14] Sanchez-Palencia E, Zaoui A. Homogenization techniques for composite media. Berlin: Springer-Verlag; 1987. [15] Bakhvalov NS, Panasenko GP. Homogenization: averaging processes in periodic media: mathematical problems in the mechanics of composite materials. Kluwer Academic Publishers; 1989. [16] Cheng G, Cai Y, Xu L. Novel implementation of homogenization method to predict effective properties of periodic materials. Acta Mech Sin 2013;29:550–6. [17] Cai Y, Xu L, Cheng G. Novel numerical implementation of asymptotic homogenization method for periodic plate structures. Int J Solids Struct 2014;51:284–92. [18] Kolpakov AG. Calculation of the characteristics of thin elastic rods with a periodic structure. J Appl Math Mech 1991;55:358–65. [19] Kolpakov AG. Stressed composite structures: homogenized models for thinwalled non-homogeneous structures with initial stresses. Berlin, New York: Springer-Verlag; 2004. [20] Kalamkarov AL, Kolpakov AG. Analysis, design and optimization of composite structures. Chichester, New York: John Wiley & Sons; 1997. [21] Yi S, Cheng G, Xu L. FEM formulation of homogenization method for effective properties of periodic heterogeneous beam and size effect of basic cell in thickness direction. Comput Struct 2015;156:1–11. [22] Goldberg DE. Genetic algorithms in search, optimization, and machine learning. New York: Addison-Wesley; 1989. [23] Dorigo M, Gambardella LM. Ant colony system: a cooperative learning approach to the traveling salesman problem. IEEE Trans Evol Comput 1997;1:53–66. [24] Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220:671–80. [25] Cheng G, Liu Y. A new computation scheme for sensitivity analysis. Eng Optim 1987;12:219–34. [26] Zhang WH, Domaszewski M, Bassir H. Developments of sizing sensitivity analysis with the ABAQUS code. Struct Optim 1999;17:219–25. [27] Jeong SH, Park SH, Choi D-H, Yoon GH. Toward a stress-based topology optimization procedure with indirect calculation of internal finite element information. Comput Math Appl 2013;66:1065–81. [28] Henrichsen SR, Lindgaard E, Lund E. Free material stiffness design of laminated composite structures using commercial finite element analysis codes. Struct Multidisc Optim 2014. on line. [29] Bendsøe MP. Optimal shape design as a material distribution problem. Struct Optim 1989;1:193–202. [30] Stolpe M, Svanberg K. An alternative interpolation scheme for minimum compliance topology optimization. Struct Multidisc Optim 2001;22:116–24. [31] Qiu G, Li X. Design of materials with prescribed elastic properties using Dfunctions. Chinese J Solid Mech 2008;29:250–5. [32] Svanberg K. A class of globally convergent optimization methods based on conservative convex sepatable approximations. Soc Indust Appl Math 2002;12:555–73. [33] Svanberg K. Method of moving asymptotes – a new method for structural optimization. Int J Numer Methods Eng 1987;24:359–73. [34] Xu S, Cai Y, Cheng G. Volume preserving nonlinear density filter based on heaviside functions. Struct Multidisc Optim 2010;41:495–505.