Multi-objective free-form optimization for shape and thickness of shell structures with composite materials

Composite Structures 135 (2016) 262–275

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Multi-objective free-form optimization for shape and thickness of shell structures with composite materials Kenichi Ikeya a, Masatoshi Shimoda b,⇑, Jin-Xing Shi b a b

Graduate School of Advanced Science and Technology, Toyota Technological Institute, 2-12-1 Hisakata, Tenpaku-ku, Nagoya, Aichi 468-8511, Japan Department of Advanced Science and Technology, Toyota Technological Institute, 2-12-1 Hisakata, Tenpaku-ku, Nagoya, Aichi 468-8511, Japan

a r t i c l e

i n f o

Article history: Available online 21 September 2015 Keywords: Composite material Free-form Shape optimization Shell Thickness optimization Orthotropic material

a b s t r a c t In this paper, we present a two-phase optimization method for designing the shape and thickness of shell structures consisting of orthotropic materials. We consider a multi-objective in terms of the compliances under multi boundary conditions, and use the weighted sum compliance as the objective functional and minimize it under the volume and the state equation constraints. The 1st phase is shape optimization, in which a shell structure is varied in the out-of-plane direction to the surface to create its optimal shape. In the 2nd phase, thickness optimization is implemented after shape optimization to decrease the compliance further. A free-form shape and thickness optimization problem is formulated in a distributed-parameter system based on the variational method. The shape and thickness sensitivities are theoretically derived and applied to the H1 gradient method for shape and size optimization. The optimal multi-objective free-form of a shell structure with an orthotropic material can be determined using the proposal method, and the influence of the orthotropic angle to the optimal shape and thickness distribution is investigated in detail. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Shell structures are widely used in various industrial products. From an economic point of view, weight reduction is strictly required in the structural design of cars, aircrafts and so on. The usage of composite materials in shell structures is one of the solutions to meet the requirement since they have higher material performances than metals. In especial, orthotropic materials can be used for making specific stiff directions of shell structures. Moreover, with design optimization, mechanical properties can be significantly improved. In design optimization of shell structures with orthotropic materials, three structural optimization methods, topology optimization, shape optimization and size optimization, are commonly considered among designers. Most of the previous works focused on the topology optimization of shell structures with composite materials. Matteo and Pierre [1] studied topology optimization for minimizing the weight of structures under compliance and stress constraints. The investigation of material influence has also been demonstrated by several authors. Tenek and Hagiwara [2] investigated the isotropic and anisotropic influence of plates for ⇑ Corresponding author. Tel.: +81 52 809 1782. E-mail address: [email protected] (M. Shimoda). http://dx.doi.org/10.1016/j.compstruct.2015.09.011 0263-8223/Ó 2015 Elsevier Ltd. All rights reserved.

minimizing their weight. Blaspues and Stolpe [3,4] treated orthotropic beams with several material distribution patterns to find the influence of material characteristics. Furthermore, size optimization for composite materials has also been reported by some scholars. Li and Yoshihiro [5] designed optimal laminated shallow shell structures under variable material orientations. Multi-objective optimization of composite laminate was proposed by Nik et al. [6]. It was found that by allowing the fiber orientation perpendicular to the loading direction, both stiffness and buckling behavior under volume constraint can be improved and the Pareto optimal design space can be found using the genetic algorithm. Zhu et al. [7] presented a multi-objective optimization combining weight, thermal expansion and buckling load by varying ply thickness of CFRP beams using the genetic algorithm. Shape optimization of shell structures, including parametric and non-parametric methods, is also an effective means in optimal design for shell structures under the material determination. Most of the previous works of shape optimization design of shell structures are classified as parametric methods, in which a shell structure is parameterized using CAD parameters such as height and thickness, parametric surface or design elements in advance. The parametric method is effective for reducing the number of design variables and seldom causes a jagged surface problem. However, designers need considerable knowledge and experience

263

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

of shape parameterization. The parametric methods have been studied not only for designing simple formed shell structures but also for arbitrarily formed shell structures. Ramm et al. [8] introduced the Bezier function and used its control points as the design variables. Rao and Hinton [9] used the Coons patch in shape optimization. Ugail and Wilson [10] obtained a smooth shell surface with the partial differential equation method. In contrast, non-parametric techniques for designing free-form shell structures were also proposed previously. Daoud et al. [11] proposed a non-parametric method with an artificial filtering technique for maintaining smoothness. Leiva [12] proposed a non-parametric method based on the basis vector method, but without smoothing. Some other works using non-parametric methods that only dealt with shape sensitivity analysis of shell structures were also carried out [13–15]. The size optimization involving thickness optimization methods have been developed by several authors. For example, Tomás and Martí [16] introduced shape and size optimization of concrete shell structures used control points as the design variables. Velea at el. [17] considered thickness of the FRP structures by a combined method. The above two methods, are classified as parametric method which needs designers’ considerable knowledge and experience of thickness parameterization. On the other hand, non-parametric methods for thickness optimization have been also developed by some scholars. Arnout at el. [18] proposed a nonparametric method for thickness optimization considering stress with Kreisselmeier–Steinhauser function. Czarnecki and Lewinski [19] proposed a non-parametric method based on a potential method. Both of the above two works used the finite dimensional approach, however, could not assure smoothness of structures. The free-form optimization method for shells is one of the nonparametric methods for arbitrarily formed shell structures that can determine the optimal smooth and natural free-form without causing jagged surfaces and without requiring shape parameterization. This method was proposed based on the traction method or H1 gradient method [20–23]. However, there has seldom study of shape optimization for shell structures consisting of orthotropic materials. In our previous work [23], we developed a free-form optimization method for determining free-form optimal shell structures with isotropic material. In the present work, the method is extended to a shell structure with an orthotropic material, and the influence of the orthotropic angle is investigated. Moreover, in this study, a non-parametric method for thickness distribution based on a gradient method is newly developed introducing Poisson’s equation both to reduce the objective functional and maintain thickness smoothness. That is, the shape and thickness optimization method is integrated to obtain higher stiffness of shell structures, or to eliminate waste of the material. The key point of the integration of the two-phase optimization of shell structures is determining the shape at first, and reducing unnecessary thickness for lighting-weight subsequently. In addition, shell structures with multi boundary conditions are considered for actual applications. The term of multi boundary conditions here refer to more than two boundary conditions that act independently on a shell structure. In the present work, we

use the weighted sum compliance under multi boundary conditions as the objective functional. The compliance minimization problem is formulated in a distributed-parameter shape and thickness optimization system. The sensitivity functions, also called the shape and thickness gradient functions or the optimal conditions, are theoretically derived using the material derivative method and the adjoint method. The derived shape gradient functions are applied to the proposed two-phase optimization method. In the following sections, the governing equation for shell structures, the formulation of the compliance minimization problem, the two-phase optimization method and the H1 gradient methods for shape and thickness optimization are described in detail. At last, numerical analysis results for several design problems are presented to demonstrate the effectiveness and practical utility of our solution. In addition, the influence of the orthotropic angle is investigated by changing the material direction of the initial shape. 2. Governing equation for a shell structure with an orthotropic material As shown in Fig. 1(a) and Eqs. (1)–(3), we consider a shell structure owning an initial bounded domain X
R3 with the boundary @ X, mid-surface A (with the boundary @A), side surface S and thickness h. It is assumed for simplicity that a shell structure occupying a bounded domain is a set of infinitesimal flat surfaces. The notation dA expresses a small area.






t t 2 2



X ¼ fðx1 ; x2 ; x3 Þ 2R3 jðx1 ; x2 Þ 2 A
R2 ; x3 2  ;


t t 2 2

ð1Þ



X¼A  ;

ð2Þ


 t t S ¼ @A   ; 2 2

ð3Þ

It is assumed that the mapping of the local coordinate system ðx1 ; x2 ; 0Þ, which gives the position of the mid-surface of the shell structure, to the global coordinate system ðX 1 ; X 2 ; X 3 Þ, i.e., / : ðx1 ; x2 ; 0Þ 2 R3 #ðX 1 ; X 2 ; X 3 Þ 2 R3 , is piecewise smooth. The Mindlin–Reissner plate theory is applied concerning bending of plates. The coupling of the membrane stiffness and bending stiffness is ignored for simplicity. Using the sign convention in Fig. 1 (b), the displacements expressed by the local coordinates u ¼ fui gi¼1;2;3 are considered by dividing them into the displacements in the in-plane direction fua ga¼1;2 and in the out-of-plane direction u3 , shown as:

ua ðx1 ; x2 ; x3 Þ ¼ u0a ðx1 ; x2 Þ  x3 ha ðx1 ; x2 Þ;

ð4Þ

u3 ðx1 ; x2 ; x3 Þ ¼ wðx1 ; x2 Þ:

ð5Þ

The notations fu0a ga¼1;2 , w and fha ga¼1;2 express the in-plane displacements, out-of-plane displacement and rotational angles of the mid-surface of the shell structure, respectively.

Fig. 1. Shell as a set of infinitesimal flat surfaces. (a) Geometry of shell and global coordinates and (b) Local coordinates and DOF of flat surface.

264

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

When NB boundary conditions are independently applied to a shell structure, the weak form of the nth state equation with

3. Domain and material variation for designing free-form shell structure with an orthotropic material

ðnÞ

respect to ðu0 ; wðnÞ ; hðnÞ Þ 2 U ðnÞ ; ðn ¼ 1; . . . ; NBÞ can be expressed as Eq. (6) by substituting Eqs. (4) and (5) into the variational equation of the three-dimensional linear elastic body, considering ðnÞ rðnÞ 33 ¼ 0 and eliminating e33 . Eq. (7) can also be expressed as Eq. (8)

by considering the relationship of

R

X ðÞdX

¼

R R t=2 A

t=2

ðÞdzdA.

ðnÞ  ðnÞ  ðnÞ  ðnÞ ðnÞ  ðnÞ ðnÞ aððu0 ; wðnÞ ; hðnÞ Þ; ðu 0 ; w ; h ÞÞ ¼ lððu0 ; w ; h ÞÞ;  ðnÞ ; w  ðnÞ ; hðnÞ Þ 2 U ðnÞ ; ðuðnÞ ; wðnÞ ; hðnÞ Þ 2 U ðnÞ ; n ¼ 1; . . . ; NB; 8ð u 0

0

ð6Þ where the energy bilinear form að; Þ and the linear form lðÞ for the ðnÞ

nth state variables ðu0 ; wðnÞ ; hðnÞ Þ are respectively defined as:

0 ; w  ðnÞ ; hðnÞ ÞÞ ¼ aððu0 ; wðnÞ ; hðnÞ Þ; ðu ðnÞ

ðnÞ

Z

ðnÞ

X

ðnÞ

ðnÞ

 0c;d fEabcd ðu0a;b  x3 ha;b Þðu

ðnÞ    x3 hc;d Þ þ ESab ðwðnÞ ;a  ha Þðw;b  hb ÞgdX; ðnÞ

ðnÞ

ðnÞ

ð7Þ


Z  S ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ  a;b þ eM   u þ ke c c dA ; ¼ feBabcd jc;d j u abcd 0c;d 0a;b ab b a

ð8Þ

A

0 ; w  ðnÞ ; hðnÞ ÞÞ ¼ lððu ðnÞ

Z

 M a ha þ Q

@AðnÞ

ðnÞ

 ðnÞ Þds; w

ð9Þ

where the notations fu0a ga¼1;2 , w and fha ga¼1;2 express the in-plane displacements, out-of-plane displacement and rotational angles of the mid-surface of the shell structure, respectively. In this paper, the tensor subscript notation uses Einstein’s summation convention and a partial differential notation for the spatial coordinates ðÞ;i ¼ @ðÞ=@xi . The notation ðÞ denotes a variation. Loads relative to the local coordinate system ðx1 ; x2 ; 0Þ are defined as: f ðnÞ ¼ ðnÞ

ðnÞ

mðnÞ ¼ fma ga¼1;2 ,

ff a ga¼1;2 , ðnÞ

(n)

fMa ga¼1;2 , Q

ðnÞ

and tb

¼

(n)

q ,

ðnÞ ftbi gi¼1;2;3 ,

N

ðnÞ

ðnÞ

¼ fN a ga¼1;2 ,

M

ðnÞ

¼

which denote in-plane loads,

out-of-plane moments, out-of-plane load on the sub-domains in A, in-plane loads, bending moments, shearing force on the subboundaries in @A and body forces on A, respectively. In addition, fEabcd ga;b;c;d¼1;2 and fESab ga;b¼1;2 express an orthotropic plane stress elastic tensor and an orthotropic elastic tensor with respect to the shearing, respectively. feBabcd ga;b;c;d¼1;2 , feSab ga;b¼1;2 and feM abcd ga;b;c;d¼1;2 express orthotropic elastic tensors with respect to bending, shear and membrane component, respectively. The constant k expresses a shear correction factor (assuming k = 5/6). ðnÞ

of domain variation. Then, X s  T s ðXÞ; Xs  TðXÞ. If T 1 exists and by thinking of s as time, the design velocity field (i.e., the amount of domain variation per iteration) can be given as:

VðX s Þ 

@X s @T s 1 ¼ ðT ðX s ÞÞ: @s @s s

ð13Þ

When the mapping T s ðXÞ is assumed to be regular enough in the neighbor of s, it can be expanded using the Taylor series around s. Then, by ignoring higher-order terms, the infinitesimal domain variation can be given as:

T sþDs ðXÞ ¼ T s ðXÞ þ DsVðXÞ;

ðnÞ ðnÞ ðnÞ  ðnÞ ðnÞ  0a  mðnÞ ðf a u a ha þ q w ÞdA A Z Z ðnÞ ðnÞ 0a þ b3ðnÞ w  ðnÞ ÞdA þ  ðnÞ tðba u ðNðnÞ þ a u0a ds ðnÞ

AðnÞ ðnÞ ðnÞ

Domain variation using the velocity analysis is introduced briefly for formulating the free-form optimization problem. A detailed explanation of this method has been expressed in the references [24–26]. As shown in Fig. 2, we consider that a shell structure undergoes a domain variation in the out-of-plane direction to the surface such that the initial bounded domain X, mid-surface A and side surface S become Xs ð XðsÞÞ, As and Ss, respectively. It is assumed that the thickness h keeps a constant with respect to the domain variation. This domain variation is expressed by a mapping T : X ! X s ðXÞ; X 2 X. The subscript s indicates the iteration history

ð14Þ

The optimal design velocity field VðXÞ is determined by the H1 gradient method for shell structures which will be explained later. The orthotropic material with E1 and E2 is assumed to move with the surface with regard to the shape variation as shown in Fig. 2. The material directions are not design variables to be determined, but are automatically determined according to the shape variation. Note that, with our method, the material directions vary smoothly rather than suddenly in the process of creating free-form surface. Therefore, the problem of sudden change of the strain and stress can be avoided, and thus the shear locking problem can be ignored. 4. Multi-objective free-form optimization for shape-thickness design considered orthotropic material In this study, with the aim of maximizing the stiffness of a shell structure under multi boundary conditions, the compliance vector ð1Þ

ð2Þ

ðNBÞ

flððu0 ; wð1Þ ; hð1Þ ÞÞ; lððu0 ; wð2Þ ; hð2Þ ÞÞ; . . . ; lððu0 ; wðNBÞ ; hðNBÞ ÞÞg is used as an index of structural stiffness under multi boundary conditions. This objective functional is scalarized by the weighted sum method given by: ðnÞ NB X lððu0 ; wðnÞ ; hðnÞ ÞÞ cðnÞ ; ðnÞ linit n¼1

ðnÞ

The notations fjab ga;b¼1;2 and fca ga¼1;2 express the curvatures and the transverse shear strains which are defined as:

1 2

ðnÞ ðnÞ jðnÞ ab  ðha;b þ hb;a Þ;

ð10Þ

ðnÞ ðnÞ cðnÞ a  w;a  ha :

ð11Þ

It should be noted that U ðnÞ in Eq. (6) is given by: ðnÞ

ðnÞ

ðnÞ

ðnÞ

5

U ðnÞ ¼ fðu01 ; u02 ; wðnÞ ; h1 ; h2 Þ 2 ðH1 ðAÞÞ j satisfy the given Dirichlet conditions on each subboundaryg; where H1 is the Sobolev space of order 1.

ð12Þ Fig. 2. Shape variation in the out-of plane by V.

ð15Þ

265

K. Ikeya et al. / Composite Structures 135 (2016) 262–275 NB X cðnÞ ¼ 1;

ðnÞ

ð16Þ

þ Hf a

n¼1 ðnÞ

where linit indicates the compliance of the nth boundary condition of the initial shape, which is used for normalizing the compliance. The notation c(n) indicates the weighting coefficient of the nth boundary condition, and has the relationship given by Eq. (16). Letting the volume and the state equations in Eq. (6) be the constraint conditions and utilizing the weighted sum compliance in Eq. (15) as the objective functional, a distributed-parameter shape optimization problem for finding the optimal design velocity field V (or As ) and the optimal thickness distribution can be formulated as:

Given A; t

ð17Þ

find Vðor As Þ; ts

ð18Þ

that minimizes

NB X

ðnÞ

cðnÞ

lððu0 ; wðnÞ ; hðnÞ ÞÞ ðnÞ

linit

n¼1


Z  ^ tdA 6 M subject to M ¼

;

ð19Þ

and Eq: ð6Þ

ð20Þ

A

^ denote the volume and its constraint value, where M and M respectively.  ðnÞ  ðnÞ ðnÞ Letting ðu 0 ; w ; h Þ and K denote the Lagrange multipliers for the nth state equation and the volume constraint, respectively, the Lagrange functional L associated with this problem can be expressed as:

 0ð1Þ ; w  ð1Þ ; hð1Þ Þ; ðu

ð1Þ ; ðu0 ; wð1Þ ; hð1Þ Þ;

LðX

 ðNBÞ  ðNBÞ ; hðNBÞ Þ; KÞ ¼ ;w ðu 0 þ

ð2Þ ðu0 ; wð2Þ ; hð2Þ Þ; . . . ;

cðnÞ ðnÞ linit

cðnÞ

ðnÞ

þ Hq þ

hGt ; t it ¼

cðnÞ

ðnÞ Hb3

Z

0

ðnÞ

linit ðnÞ

linit

0

Gt t dA ¼ A

þ

! ðnÞ

@t

 HmðnÞ a !

ðnÞ

w

ðnÞ

w

þw

!

cðnÞ

ðnÞ u ðnÞ 0a linit

 ðnÞ þu 0a

)

Z (X NB @EM abcd @t

ð23Þ

 0a;b u0c;d u

@ESab c cb @t a

!

) þ K t 0 dA;

ð24Þ

where C H is the suitably smooth function space that satisfies the constraints of domain variation. In the derivation of Eq. (24), the relationship of ðV  ntop Þntop ¼ ðV  nbtm Þnbtm is assumed. As shown in Fig. 3, the notations ntop ; nbtm denote unit outward normal vectors at the top and bottom surface of a virtual domain in the thickness direction, respectively, and a unit normal vector at the midsurface nmid ð nÞ ¼ ntop ¼ nbtm is assumed. The notation H denotes twice of the mean curvature of the surface. Eqs. (23) and (24) express the effects of the area variation by V and the thickness variation, respectively. The term with respect to the volume variation is included in Eqs. (23) and (24) as KtHn  V and Kt 0 , separately. The optimal conditions of the Lagrange functional L with ðnÞ

respect to the state variables ðu0 ; wðnÞ ; hðnÞ Þ, the adjoint variables  ðnÞ  ðnÞ ðnÞ ðu 0 ; w ; h Þ and K are derived as:

  ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ ðnÞ  00 ; w  00 ; w  0 ; h0 Þ ¼ lððu  0 ; h0 ÞÞ; a ðu0 ; wðnÞ ; hðnÞ Þ; ðu  00 ; w  0 ; h0 Þ 2 U ðnÞ ; ðu0ðnÞ ; wðnÞ ; hðnÞ Þ 2 U ðnÞ ; 8ðu ðnÞ

ðnÞ NB X lððu0 ; wðnÞ ; hðnÞ ÞÞ cðnÞ ðnÞ linit n¼1

ðnÞ hðnÞ a þ ha

þ KtH n  VdA;

þw

hc;d ha;b þ j0

ðnÞ

linit

þ Hba

 ðnÞ

n¼1

cðnÞ

ðnÞ

 ðnÞ

!#

A

@EBabcd

ðnÞ

 0a u0a þ u

!

ðnÞ

ðnÞ

n ¼ 1;    ; NB; ð25Þ

NB X  ðnÞ  0ðnÞ ; w  ðnÞ ; hðnÞ Þ  aððu0ðnÞ ; wðnÞ ; hðnÞ Þ; ðu  ðnÞ ðnÞ flððu 0 ; w ; h ÞÞg

ðnÞ

n¼1

^ þ KðM  MÞ

ð21Þ

ðnÞ

ðnÞ

0 0 0 ðnÞ ðnÞ ðnÞ ðnÞ lððu0 ; w ; h ÞÞ  ðnÞ  ðnÞ ðnÞ aððu00 ; w0 ; h0 Þ; ðu ; 0 ; w ; h ÞÞ ¼ c ðnÞ linit ðnÞ ðnÞ ðnÞ  ðnÞ ; w  ðnÞ ; hðnÞ Þ 2 U ðnÞ ; n ¼ 1; . . . ; NB; 8ðu0 ; w0 ; h0 Þ 2 U ðnÞ ; ðu 0

0

ð26Þ

For the sake of simplicity, it is assumed that the sub-boundaries ðnÞ

ðnÞ

ðnÞ

on which the non-zero boundary forces N ; Q ; M do not vary (i.e., V ¼ 0) and the surface forces f ðnÞ ; mðnÞ ; qðnÞ and the body forces ðnÞ

tb do not vary with regard to both time and space (i.e., _ ðnÞ ¼ tb_ ðnÞ ¼ 0; q_ ðnÞ ¼ 0). It is also assumed that the direction f_ ðnÞ ¼ m of the orthotropic material is fixed on its local coordinate system (E_ ¼ 0). Then, using the design velocity field V and t0 to represent the amount of domain variation and thickness variation, the material derivative L_ [1,25] of the Lagrange functional L involving the thickness variation can be derived as:

L_ ¼

NB X

ðnÞ

c

ðnÞ

ðnÞ

ðnÞ

linit

n¼1



ðnÞ

lððu00 ; w0 ; h0 ÞÞ

^ ¼ 0; KðM  MÞ

^ 6 0; K P 0; MM

ð27Þ ðnÞ

The state equations in Eq. (25) for ðu0 ; wðnÞ ; hðnÞ Þ and the adjoint  ðnÞ  ðnÞ ðnÞ equations in Eq. (26) for ðu 0 ; w ; h Þ can be solved using a standard commercial FEM code. Considering the following quasi self-adjoint relationships ðnÞ

c ðnÞ ðnÞ ðnÞ  ðnÞ  ðnÞ ðnÞ ðu 0 ; w ; h Þ ¼ ðnÞ ðu0 ; w ; h Þ; linit

n ¼ 1; . . . ; NB ðnÞ

obtained from Eqs. (25) and (26), and substituting ðu0 ; wðnÞ ; hðnÞ Þ and K determined by these optimality conditions into Eq. (22),

NB X ðnÞ ðnÞ ðnÞ  00 ; w  0 ; h0 Þ þ flðu n¼1

ðnÞ ðnÞ ðnÞ ðnÞ  00 ; w  0 ; h0 ÞÞ aððu0 ; wðnÞ ; hðnÞ Þ; ðu

^ 0 ; w  ðnÞ ; hðnÞ ÞÞg þ K0 ðM  MÞ  aððu00 ; w0 ; h0 Þ; ðu ðnÞ

ðnÞ

ðnÞ

þ hGA n; ViS þ hGt ; t 0 it ;

ðnÞ

V 2 CH:

ð28Þ

ð22Þ

where


 Z (X NB   t ðnÞ ðnÞ hGA n; ViS ¼ GA V n dA ¼  Eabcd u0a;b þ ha;b 2 A A n¼1


 t ðnÞ t ðnÞ t ðnÞ ðnÞ ðnÞ  ðnÞ  h h  u  E h þ u   u a b c d 0c;d 0a;b 0c;d 2 c;d 2 a;b 2 c;d Z

Fig. 3. Virtual domain in the thickness direction.

266

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

the material derivative L_ can be expressed as the sum of the dot product of the shape gradient function GA , the design velocity field in the normal direction to the surface V n ð¼ V  nÞ, the dot product of the thickness gradient function Gt and the thickness variation field t 0 , shown as:

L_ ¼ hGA n; ViS þ hGt ; t0 it ¼

Z

Z

GA V n dA þ A

Gt t0 dA;

ð29Þ

A

The shape and thickness gradient functions GA ð¼ GA nÞ and Gt (i.e., sensitivity functions) for this problem are derived as

( " 

 NB X cðnÞ t ðnÞ t ðnÞ ðnÞ ðnÞ u0c;b þ hc;b GA ¼  ðnÞ Eabcd u0a;b þ ha;b 2 2 l n¼1
init 
  t ðnÞ t ðnÞ ðnÞ ðnÞ  Eabcd u0a;b  ha;b u0c;b  hc;b þ KtH n; 2 2

ð30Þ

( !) B NB X @EM @ESab cðnÞ @Eabcd abcd Gt ¼ cc þ hc;d ha;b þ j0 u0c;d u0a;b þ K; ðnÞ @t @t a b @t n¼1 l init

ð31Þ where H can be determined by a proposed method [1]. The derived gradient functions are applied to the H1 gradient method to determine the optimal design velocity field V and the optimal thickness variation field t0 . 5. H1 gradient method for shell structures The free-form optimization method for shells was proposed by Shimoda et al. [24], which consists of main three processes: (1) Derivation of shape gradient functions; (2) Numerical calculation of shape gradient functions; (3) The H1 gradient method for determining the optimal shape variation. The H1 gradient method is a gradient method in a Hilbert space. The original H1 gradient method was proposed by Azegami for solid structures and also called the traction method [20,26,27]. Shimoda modified the original optimization method for shell structures. In the present paper, we newly propose a H1 gradient method for determining the optimal thickness distribution and integrate it for shape optimization. The proposed method is a node-based shape and thickness optimization method that treats all nodes as design variables and does not require any design variable parameterization. 5.1. H1 gradient method for shape optimization of shell structures In this study, we apply the free-form optimization method for shells to a shell structure with orthotropic material. This makes it possible to obtain the optimal free-form of shell structures. In this method, the negative shape gradient function GA nðXÞ is applied as a distributed force to the fictitious-elastic shell structure in the normal direction to the design surface under a Robin boundary condition, i.e., an elastic support condition with a distributed

spring constant aða > 0Þ as shown in Fig. 4. The shape gradient function is not used directly to vary the shape but is replaced by a distributed force. a is introduced to control the influence range of the shape gradient function, and is called a coefficient of sensitivity influence. This approach makes it possible both to reduce the objective functional and to maintain surface smoothness, i.e., mesh regularity. The design velocity field V(X) is calculated as the displacement field obtained in this fictitious-elastic shell analysis, which is called ‘‘velocity analysis” and used to update the initial shape. This method is explained in a view of the gradient method as follows. When the state and adjoint equations are satisfied and V ¼ hV; hiT is considered, the perturbation expansion of the Lagrange functional L neglecting hGn; DSðV; hÞit can be written as

DL ¼ hGn; DSðV; hÞis þ OðjDsj2 Þ

ð32Þ

Considering the design velocities V ¼ fV i gi¼1;2;3 as a combination of the in-plane velocities fV 0a ga¼1;2 and the out-of-plane velocity V 3 , the governing equation of the velocity analysis for V ¼ fV 01 ; V 02 ; V 3 g is expressed as Eq. (32).

 0 ; w;  0 ; w;  hÞÞ þ aA hðV  nÞn; ðu  hÞis aððV 0a ; V 3 ; hÞ; ðu  0 ; w;  0 ; w;  hÞis ; ðV 0a ; V 3 ; hÞ 2 C H ; 8ðu  hÞ 2 C H ; ¼ hGn; ðu

ð33Þ

5

C H ¼ fðV 01 ; V 02 ; V 3 ; h1 ; h2 Þ 2 ðH1 ðAÞÞ satisfy Dirichlet condition for shape variation on Sg;

ð34Þ

The stiffness tensor in Eq. (33) functions as a smoother for maintaining mesh regularity, and its positive definiteness is the necessary condition for minimizing the objective functional. Substituting Eq. (33) into Eq. (32) and taking into account  0 ; w;   the positive definiteness of aððV 0a ; V 3 ; hÞ; ðu hÞÞ and   aA hðV  nÞn; ðu 0 ; w; hÞi, based on the positive definiteness of the elastic tensors Eabcd and ESab , the following relationship is obtained when Ds is sufficiently small:

DL ¼ aððV; hÞ; DsðV; hÞÞ  aA hðV  nÞn; DsðV; hÞis < 0:

ð35Þ

In problems where convexity is assured, this relationship definitely reduces the Lagrange functional in the process of updating the shell shape using the design velocity field V determined by Eq. (32). 5.2. H1 gradient method for thickness optimization of shell structures The H1 gradient method is extended to the thickness optimization in the present work. When the state equations and the adjoint equations are satisfied, the perturbation expansion of the Lagrange functional L neglecting hGt ; Dst 0 is can be written as:

DL ¼ hGt ; Dst0 it þ OðjDsj2 Þ

Fig. 4. Schematic of free-form optimization method for shell (H1 gradient method for shell).

ð36Þ

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

when Ds is sufficiently small, to obtain the optimal thickness variation field t0 , the following weak formed Poisson’s equation for t 0 is introduced.

at ðt 0 ; v Þ þ at ht 0  t 0 ; v it ¼ hGt ; Dst0 it ; where

at ðt 0 ; v Þ ¼

Z A

8v 2 C t ; t 0 2 C t

t ;i kij v ;j dA

ð37Þ

ð38Þ

where t 0 and t 0 denote the thickness variation field and the reference thickness, respectively. It is assumed that t0  t 0 > 0. The notations at and kij are equivalent to the heat transfer coefficient and the thermal conductivity tensor in the heat transfer analysis, respectively. Eq. (37) can be solved with a standard commercial FEM code. The kinematic admissible function space C t is defined as:



C t ¼ fv 2 H1 satisfy Dirichlet condition for thickness variationg

ð39Þ Substituting Eq. (37) into Eq. (36) and taking into account the voluntariness of m in Eq. (37), we obtain:

DL ’ hGt ; Dst 0 it ¼ ðat ðt0 ; v Þ þ at ht0  t0 ; v iÞ

ð40Þ

Furthermore, taking into account the positive definitiveness of at ðt0 ; v Þ > 0 and at ht0  t0 ; v i > 0 in Eq. (37), we have:

DL < 0

ð41Þ

In problems where convexity is assured, this relationship definitely reduces the Lagrange functional in the process of updating the thickness using the thickness variation field t0 determined by Eq. (37). In this method, the negative thickness gradient function Gt is applied as a distributed internal heat generation to a fictitious-elastic shell structure to the design surface. The thickness variation field t 0 is calculated as the temperature field of the Poisson’s equation and is used to update the initial thickness. Note that thickness is a discrete value since Pregregs or woven cloths are laminated in composite materials by the actual production line. In this case, determination of ply numbers, material distributions and material constants is important and optimization methods which treat these as design variables have been practically applied [17]. These optimization results can be determined directly by using parametric optimization methods, however, computational costs have to be considered at this time due to the large number of design variables. The method that the

267

final thickness could be rounded to the nearest multiple of ply thickness is applied practically. In this study, the actual numbers of ply thickness or material distributions of each layer are not directly treated as design variables and the materials of the structures are limited to the isotropic and the orthotropic materials, however, we aim to obtain the optimal and continuous results of shapes and thickness distributions of shells by using the method based on the variational method. Although it is difficult to make the optimized shapes and thickness distributions in practical production at this stage, the theoretical solutions are helpful for the designers to grasp the limitation properties of structures. Hence, the proposed optimization method and the calculated numerical examples could be important and meaningful for the engineering and industrial applications. Based on the proposed optimization method, we construct the two-phase optimization process by combining a standard commercial FEM code. Eqs. (25), (33) and (37) are solved by the FEM code and the shape gradient functions in Eqs. (30) and (31) are calculated without differentiating the stiffness matrix with respect to design variables. 5.3. Two-phase optimization method To minimize the compliance and mass of shell structures, both shape and thickness are treated as design variables in the optimization. In the present work, the shape optimization is applied firstly to shell structures composed of orthotropic materials. Then, the thickness optimization is carried out to the structures after shape convergence. The design optimization process is shown Fig. 5. The shape optimization process is summarized as: (1) Stiffness analysis by Eq. (29) and evaluation of objective and constraint functionals; (2) Calculation of shape gradient function by Eq. (30); (3) Velocity analysis; (4) Updating of the shape by using V; (5) Satisfying the volume constraint. This process is repeated until the optimal shape is obtained. The thickness optimization process following the shape optimization is summarized as: (6) Stiffness analysis by Eq. (29) and evaluation of objective and constraint functionals; (7) Calculation of thickness gradient function by Eq. (31); (8) Velocity analysis; (5) Updating of the thickness by using t0 ; (9) Satisfying the volume constraint. This process is repeated until the optimal thickness distribution is obtained. In this optimization process, the computational cost for each phase can be evaluated by the number of iterations. However, the number of iterations to the final design depends on the step size for

Fig. 5. Flowchart of the optimization process.

268

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

the shape or the thickness change (Ds in Eqs.(35) and (37)). This parameter depends on the design problem and set empirically by the designer, which is the same as the other optimization methods.

E1 ¼ 210 GPa; E2 ¼ 21 GPa, shear modulus G12 ¼ 65 GPa and Passion’s ratio mð¼ m1 ¼ m2 Þ is 0.3. 6.1. Square shell model

6. Numerical results The proposed design optimization method is applied to two design problems. Each FEM model is discretised by constant strain triangular elements. The material constants are used as:

The first design problem concerns the two-objective design of a square shell with an orthotropic material. The optimization results with isotropic material are shown for a comparison. Two-objective indicates that the objective functional is a weighted sum compliance of a shell structure under two boundary conditions as shown

Fig. 6. Boundary conditions for design problem 1. (a) Stiffness analysis (Case 1), (b) Stiffness analysis (Case 2) and (c) Velocity analysis.

Fig. 7. Optimal shapes and thickness under single loading condition of square shell problem. (a) Case 1 (c(1) = 1, c(2) = 0) and (b) Case 2 (c(1) = 0, c(2) = 1).

Fig. 8. Iteration histories. (a) Case 1 and (b) Case 2.

Fig. 9. Pareto optimal shapes and thickness under multi-loading conditions of square shell problem. (a) c(1) = 0.8, (b) c(1) = 0.5 and (c) c(1) = 0.2.

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

269

in Eq. (15). The relationship of the weighting coefficient of each loading case is shown in Eq. (16). The initial shape and problem definition are illustrated in Fig. 6. In the stiffness analysis of Case 1 shown in Fig. 6(a), three corners are supported and an upward force is applied on the fourth corner (i.e., torsional condition). In the stiffness analysis of Case 2 shown in Fig. 6(b), centers of all of the 4 edges in the square shell are simply supported and a downward force is applied at the center of the shell (i.e., bending condition). In the velocity analysis, which determines the design velocity field under shape variation constrains, it is assumed that all the edges are simply supported as shown in Fig. 6(c). The notation (SPC: 123) in the figures denotes the D.O.F. numbers constrained. The weighting coefficients c(1) of Case 1 and c(2) of Case 2 have a relationship of c(1) + c(2) = 1. The volume constraints are set as 1.05 times of the initial value in shape optimization and 1 times of the optimal value in the thickness optimization after shape convergence. Note that the volume constraint value can be arbitrarily given in our method. However, in the shape optimization (first phase) of the square shell model as shown in Fig. 6(a) with the out-of-plane variation, the initial volume of the flat shell is the minimum value. Therefore, we set the volume constraint as 1.05 times of the initial volume in the shape optimization. Fig. 10. Comparison of compliance of square shell problem.

Fig. 11. Material layout for orthotropic material of square shell problem. (a) E1 45° from x-axis (Type 1) and (b) E1 Parallel to x-axis (Type 2).

Fig. 12. Optimal shapes and thickness under single loading condition of square shell problem with orthotropic materials. (a) Type 1 and (b) Type 2.

Fig. 13. Optimal shape with the material flow of square shell problem. (a) Type 1 and (b) Type 2.

270

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

6.1.1. Isotropic material First, two-phase optimization for the isotropic model with E1 is implemented for comparison. Fig. 7(a) and (b) show the optimized shapes and thickness distributions under boundary conditions of Case 1 and Case 2, respectively. Natural bead patterns with smoothness and natural thickness distributions are created according to each load condition. An X-shaped bead pattern appears on the surface in Case 1, in which thickness is distributed as the Xshaped bead due to a large amount of load bearing and the other parts become thin. In Case 2, a central point on the surface rises and the thickness is distributed as a cross- shaped along to x and y axes. The iteration histories are shown in Fig. 8, in which the values are normalized by the initial value. The red line, blue line and green line indicate the volume ratio, and the objective functional ratio during the 1st phase (shape optimization) and the 2nd phase (thickness optimization), respectively. The objective functional ratio of each case is smoothly converged, while satisfying the volume constraint. Compared with the initial model, the compliance is reduced by 97% after the shape optimization and 98% after the thickness optimization in Case 1 (Fig. 8(a)), and reduced by 99% after the shape optimization and 99.5% after the thickness optimization in Case 2 (Fig. 8(b)).

Fig. 9(a)–(c) show the Pareto optimal shapes and thickness distributions when c(1) = 0.8, c(1) = 0.5 and c(1) = 0.2, respectively. As shown in Fig. 9(a)–(c), When c(1) decreases, thickness of the center at the square shell gradually swells and the X-shaped bead gradually disappears, while the cross-shaped thickness distribution appears gradually. It is clear that a set of Pareto optimal shapes and thickness distributions (i.e., intermediate shapes and thickness distributions) can be obtained by varying the weighting coefficients. When c(1) is large, the shape and thickness distribution are strongly influenced by Case 1. With decreasing the value of c(1), the shape and thickness distribution gradually changes to show the influence of Case 2. Furthermore, when c(1) changes from 1.0 to 0.0, the compliance ratios after two-phase optimization are shown in Fig. 10. It shows that the compliance ratio of each case involves a trade-off. 6.1.2. Orthotropic material Shell structures with orthotropic materials are optimized in this section. Two types (Type 1 and Type 2) of orthotropic materials with different material distributions are expressed in Fig. 11. Each case of material distributions is set as 45° from x axis and parallel to x axis. Young’s modulus ratio is set as E1 : E2 ¼ 10 : 1. Under the

Fig. 14. Iteration histories. (a) Type 1 and (b) Type 2.

Fig. 15. Optimal shapes and thickness under single loading condition of square shell problem. (a) Case 1(c(1) = 1, c(2) = 0) and (b) Case 2(c(1) = 0, c(2) = 1).

Fig. 16. Pareto optimal shapes and thickness under multi-loading conditions of square shell problem. (a) c(1) = 0.8, (b) c(1) = 0.5 and (c) c(1) = 0.2.

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

Fig. 17. Comparison of compliance of square shell problem.

271

loading condition of Case 1, Fig. 12(a) and (b) show the optimal shapes and thickness distributions of Type 1 and Type 2, respectively. The results show that, when considering the material distributions, the smooth shapes and thickness distributions can be also obtained based on the proposed method. For emphasizing the material distributions, we express the material flows of the results of Type 1 and Type 2 in Fig. 13(a) and (b), respectively. The iteration histories of Type 1 and Type 2 are shown in Fig. 14 (a) and (b), respectively. The objective functional ratio of each case converges smoothly while satisfying the volume constraint. The results show that the compliance is reduced by 97% after the shape optimization and 98% after the thickness optimization in Type 1, and reduced by 90% after the shape optimization and 93% after the thickness optimization. Thus, the developed twophase optimization method works well for shell structures with orthotropic materials. Fig. 15(a) and (b) show the optimal shapes and thickness distributions of the shell structure with material distribution of Type 1 in Case 1 and Case 2, in which the reductions of the compliance ratios are 98% and 99%, respectively. Moreover, when c(1) equals to 0.8, 0.5 and 0.2, the Pareto optimal shapes and thickness distributions are shown in Fig. 16(a)–(c), respectively. Fig. 17 shows the compliance ratios of the above results of Type 1 with various weighting coefficients. The compliance ratio

Fig. 18. Boundary conditions for design problem 2. (a) Stiffness analysis (Case 1), (b) Stiffness analysis (Case 2) and (c) Velocity analysis.

Fig. 19. Optimal shapes and thickness under single loading condition of T-joint model. (a) Case 1(c(1) = 1, c(2) = 0) and (b) Case 2(c(1) = 0, c(2) = 1).

272

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

Based on the two-phase optimization method for shells, the second design problem concerning a T-joint shell structure is carried out. The initial shape and the problem definition are illustrated in Fig. 18. In the stiffness analyses shown in Fig. 18 (a) and (b), left and right side edges of the lower beam of the T-joint shell structure are constrained at SPC: 123. Two coupling forces are applied as Case 1 (i.e., torsional condition) at the top edge of the pillar of the shell structure and a distributed force in y direction is applied as Case 2 (i.e., bending condition). In the velocity analysis shown in Fig. 18(c), it is assumed that the edges of right and left sides, and top and bottom surfaces are simply supported. The volume constraints of both shape and thickness are set as 1.00 times of the initial value. We illustrate the optimization results of isotropic model and orthotropic model in 6.2.1 and 6.2.2, respectively.

condition. The bottom part of the pillar is expanded and the lower beam is shrunk for satisfying the volume constraint. Most of thickness is distributed on the pillar due to the twisted force. In Case 2, the pillar of T-joint is squashed in x-axis direction and the lower beam is expanded. Thickness is mostly distributed on the pillar and the side of lower beam due to the coupling forces. Fig. 20 (a) and (b) show the iteration histories of Case 1 and Case 2, respectively. The objective functional of each case is smoothly converged under the volume constraint. In Case 1 shown as Fig. 20(a), the compliance is reduced by 47% after the shape optimization and 57% after the thickness optimization. In Case 2 shown as Fig. 20(b), the compliance is reduced by 62% after the shape optimization and 73% after the thickness optimization. When c(1) is 0.8, 0.5 and 0.2, the obtained Pareto optimal shapes and thickness distributions are shown in Fig. 21(a)–(c), respectively, which are the transitive models between the models shown in Fig. 19(a) and (b). Fig. 22 expresses the compliance ratios under each boundary conditions of Case 1 and Case 2 when c(1) equals to 1.0, 0.8, 0.5, 0.2 and 0.0. The results show that the compliance can be effectively reduced in most of the cases.

6.2.1. Isotropic material Fig. 19(a) and (b) show the optimal shapes and thickness distributions of the T-joint shell structure composed of an isotropic material with E1 under the loading conditions of Case 1 and Case 2, respectively. In both Case 1 and Case 2, the smooth shapes and thickness distribution are obtained according to each load

6.2.2. Orthotropic material The same optimization problem of the T-joint shell structure with orthotropic materials as shown in Fig. 23 is carried out to investigate the effect of material distribution. The direction of E1 is expressed by the red arrow lines, and the direction of E2 is perpendicular to E1 . Based on the two-phase optimization method, the

of shell structures with orthotropic materials of each case involves a trade-off, which is the same as isotropic materials. 6.2. T-joint shell model

Fig. 20. Iteration histories. (a) Case 1 and (b) Case 2.

Fig. 21. Pareto optimal shapes and thickness under multi-loading conditions of T-joint model. (a) c(1) = 0.8, (b) c(1) = 0.5 and (c) c(1) = 0.2.

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

Fig. 22. Comparison of compliance of T-joint model.

273

optimal shapes and thickness distributions of the T-joint orthotropic model under the load conditions of Case 1 and Case 2 are shown in Fig. 24(a) and (b), respectively. Both shapes and thickness distributions express anti-symmetry according to the material distribution, and the material flows of the two cases are shown in Fig. 25, in which the material directions are automatically determined to the shape change. From the iteration histories shown in Fig. 26(a) and (b), the objective functionals are reduced by 60% after the shape optimization and 84% after the thickness optimization in Case 1, and reduced by 65% after the shape optimization and 86% after the thickness optimization in Case 2. Fig. 27 shows the three transitive optimal shapes and thickness distributions between Case 1 (c(1) = 1.0) and Case 2 (c(1) = 0.0) when c(1) decreases from 0.8 to 0.2. When c(1) is from 1.0 to 0.0, the compliance ratios of the T-joint shell structure with orthotropic materials are shown in Fig. 28. It is showed that the compliance can be effective reduced in all of the cases, which can verify the validity of the proposed two-phase optimization method well. Note that the effect of the thickness (second phase) optimization depends on the design problem. In the square shell model of Section 6.1, the reductions of compliances are almost governed by the shape (first phase) optimization, since the membranestress structure has been already created only by the shape

Fig. 23. Material layout for orthotropic material of T-joint model.

Fig. 24. Optimal shapes and thickness under single loading condition of T-joint model. (a) Case 1(c(1) = 1, c(2) = 0) and (b) Case 2(c(1) = 0, c(2) = 1).

274

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

Fig. 25. Optimal shape with the material flow of T-joint model. (a) Case 1 and (b) Case 2.

Fig. 26. Iteration histories. (a) Case 1 and (b) Case 2.

Fig. 27. Pareto optimal shapes and thickness under multi-loading conditions of T-joint model. (a) c(1) = 0.8, (b) c(1) = 0.5 and (c) c(1) = 0.2.

K. Ikeya et al. / Composite Structures 135 (2016) 262–275

275

Acknowledgements This research was supported by grants-in-aid from the Sustainable Mechanical Systems R&D Center at the Toyota Technological Institute. References

Fig. 28. Comparison of compliance of T-joint model.

optimization. In contrast, in the T-joint shell model of Section 6.2, the compliances are considerably reduced by the thickness optimization, since the stress near the neutral axis of the box structure under bending can be increased by the thickness reduction rather than the shape change. 7. Conclusions In this work, we proposed a two-phase optimization method to solve a multi-objective problem of a shell structure consisting of orthotropic materials under multi boundary conditions. With this method, the Pareto optimal free-form shapes and thickness distributions of shell structures under multi boundary conditions and volume constraint can be designed. This optimization method is one of the gradient methods in a Hilbert space, using an elastic tensor as a positive definite tensor with a smoother and applying a distributed force proportional to the sensitivity functions to vary the shape and thickness. The formulation of a distributed parameter optimization problem was presented using weighted sum compliance as the objective functional subjected to the volume constraint. The shape and thickness gradient functions were derived and applied to the H1 gradient method. Two design examples were carried out to verify the effectiveness and practical utility of this method. The proposed method makes it possible to obtain the smooth and natural Pareto optimal shapes and thickness distributions while reducing the objective functional without shape parameterization. According to this method, a global form and/or a local bead pattern can be obtained according to the boundary conditions. Several patterns of orthotropic materials were considered in the design problems and the influence of material distributions to the optimal shapes and thickness distributions of shell structures was also investigated in detail.

[1] Matteo B, Pierre D. Topology optimization for minimum weight with compliance and stress constraints. Struct Multidisc Optim 2012;46:369–84. [2] Tenek LH, Hagiwara I. Optimization of material distribution within isotropic and anisotropic plates using homogenization. Comput Methods Appl Mesh Eng 1993;50:155–67. [3] Blasques JP, Stolpe M. Multi-material topology optimization of laminate composite beam cross sections. Compos Struct 2012;94:3278–89. [4] Blasques JP. Multi-material topology optimization of laminate composite beams with eigenfrequency constraints. Compos Struct 2014;111:45–55. [5] Li J, Narita Y. Multi-objective design for aeroelastic flutter of laminated shallow shells under variable flow angles. Compos Struct 2014;111:530–9. [6] Nik MA, Fayazbakhsh K, Pasini D, Lessard L. Surrogate-based multi-objective optimization of a composite laminate with curvilinear fibers. Compos Struct 2012;94:2306–13. [7] Zhu X, He R, Lu X, Ling X, Zhu L, Liu B. A optimization technique for the composite strut using genetic algorithms. Mater Des 2015;65:482–8. [8] Ramm E, Bletzinger KU, Reitinger R. Shape optimization of shell structures. Revue Européenne des Éléments Finis 1993;2:377–98. [9] Rao NVR, Hinton E. Analysis and optimization of prismatic plate and shell structures with curved planform – II. Shape Optim Comput Struct 1994;52 (2):341–51. [10] Ugail H, Wilson MJ. Efficient shape parameterization for automatic design optimisation using a partial differential equation formulation. Comput Struct 2003;81:2601–9. [11] Daoud F, Firl M, Bletzinger KU. Filter techniques in shape optimization with CAD-free parameterization. In: Proceedings of 6th world congress of structural and multidisciplinary optimization, Rio de Janeiro, 30 May–03 June 2005. [12] Leiva PL. Freeform optimization: A new capability to perform grid by grid shape optimization of structures. In: Proceedings of 6th China-Japan-Korea joint symposium on optimization of structural and mechanical systems, Kyoto, 22–25 June 2010. [13] Lee SS, Kwak BM. Shape sensitivity analysis of thin-shell structures. Finite Elem Anal Des 1992;10:293–305. [14] Yamazaki K, Vanderplaats GN. Design sensitivity analysis with isoparametric shell elements. Struct Multidisc Optim 1993;5:152–8. [15] Taroco E, Feijoo RA. A unified approach for shape sensitivity analysis of elastic shells. Struct Multidisc Optim 2004;27:66–79. [16] Tomás A, Martí P. Shape and size optimization of concrete shells. Eng Struct 2010;32:1650–8. [17] Velea MN, Wennhage P, Zenkert D. Multi-objective optimization of vehicle bodies made of FRP sandwich structures. Compos Struct 2014;111:75–84. [18] Arnout S, Firl M, Bletzinger KU. Parameter free shape and thickness optimization considering stress response. Struct Multidisc Optim 2012;45:801–14. [19] Czarnecki S, Lewinski T. On minimum compliance problems of thin elastic plates of varying thickness. Struct Multidisc Optim 2013;48:17–31. [20] Shimoda M, Azegami H, Sakurai T. Traction method approach to optimal shape design problems. SAE 1997 Trans J Passenger Cars 1998;106:2355–65. [21] Shimoda M, Tstuji J. Non-parametric shape optimization method for rigidity design of automotive sheet metal structures. SAE J Passenger Cars 2006;2006– 01-0584:483–92. [22] Shimoda M, Liu Y. A non-parametric free-form optimization method for controlling the stiffness of automotive sheet metal structures. SAE technical paper 2013; 2013-01-0962. [23] Shimoda M, Liu Y. A non-parametric free-form optimization method for shell structures. Struct Multidisc Optim 2014;50:409–23. [24] Sokolowski J, Zolesio JP. Introduction to shape optimization. Berlin Heidelberg: Springer; 1992. [25] Choi KK, Kim NH. Structural sensitivity analysis and optimization 1: linear systems. Springer Science & Business Media; 2006. [26] Azegami H. A solution to domain optimization problems. Trans Jpn Soc Mech Engs Ser A 1994;60:1479–86 (in Japanese). [27] Shi JX, Shimoda M. Interface shape optimization of designing functionally graded sandwich structures. Compos Struct 2015;125:88–95.