Exact sensitivity analysis of stresses and lightweight design of Timoshenko composite beams

Composite Structures 143 (2016) 272–286

Contents lists available at ScienceDirect

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

Exact sensitivity analysis of stresses and lightweight design of Timoshenko composite beams Qimao Liu ⇑ Department of Civil Engineering, Aalto University, Espoo 02150, Finland

a r t i c l e

i n f o

Article history: Available online 16 February 2016 Keywords: Composite beam Lightweight design Exact sensitivity Failure Optimization

a b s t r a c t The paper describes the novel optimization techniques for lightweight design of composite beams. The optimization model is to find width and depth of composite beams to minimize the mass of beams under the stiffness, strength and delamination failure constraints. The exact formulae for displacements, stresses and their sensitivities with respect to width and depth are derived using Timoshenko continuous beam theory. The analytical stiffness, strength and delamination failure functions, and their gradients are obtained using the exact expressions of displacements and stresses, and their sensitivities. The mass and its gradient are also expressed analytically. The standard gradient-based nonlinear programming algorithms are employed to solve lightweight design problems of composite beams. The lightweight designs of composite beams are performed using the proposed optimization techniques. Three standard gradient-based optimization methods (SQP, interior-point and active-set) using exact derivatives converge to the same lightweight design. However, gradient-based algorithms using finite difference sensitivities may not lead to optimal lightweight designs. It is necessary to develop exact sensitivity analysis method instead of the difference methods for gradient-based algorithms. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Composite structures are gaining more and more successful applications in aerospace, aircraft, automobile, train, naval, and defence industries because of their high performance characteristics, such as high strength-to-weight ratio, high stiffnessto-weight ratio, superior fatigue properties and high corrosion resistance [1,2]. For example, 50% of Boeing’s 787 Dreamliner has been manufactured of advanced composites, specifically epoxy and carbon fiber [3]. As the research, knowledge and confidence on the composite structures increase, the composite structures are gradually acting as the main load-carrying components, not only as the secondary load-carrying components in the engineering application. The beams are the major transversal load-carrying member in the engineering structure systems. For example, aircraft wing, helicopter blade, wind turbine blade, robot arm, and space antenna are typical composite beam structures. The floor beams of Boeing’s 787 are made of composite material. It is the first commercial airplane to use composite floor beams. From the viewpoint of the mechanics models of a composite laminated beam, the composite laminated beams can be classified ⇑ Tel.: +358 (0) 503502879; fax: +358 (0) 947023758. E-mail address: [email protected] http://dx.doi.org/10.1016/j.compstruct.2016.02.028 0263-8223/Ó 2016 Elsevier Ltd. All rights reserved.

into Euler–Bernoulli beam model [4,5], Euler–Bernoulli beam model with torsional rotation [6,7], Timoshenko beam model [8] and Timoshenko beam model with torsional rotation [9,10]. Euler–Bernoulli beam model is only applied to the composite laminated beams with a big ratio of span to thickness. However, Timoshenko beam model is applied to the composite laminated beams with a small or big ratio of span to thickness. From the viewpoint of the number of degree of freedom of a composite laminated beam, the composite laminated beams can be treated as continuous composite beam models (infinite degree of freedom) and discrete models (finite degree of freedom). The continuous models can achieve the exact solution of the composite beams. The discrete models (finite element models) can only obtain the approximate solution of the beams. The finite element models are believed to be easier implemented than the continuous composite beam models. It is often tedious and difficult to derive the exact solution of a continuous composite beam model. Many researchers have developed the analysis and optimal design methods for the composite beam structures using discrete models. For example, Blasques and Stolpe [11] performed the maximum stiffness and minimum weight optimization of laminated composite beams using finite element approach. The fiber orientations and layer thicknesses are design variables. Cardoso et al. [12,13] used finite element technique to deal with

273

Q. Liu / Composite Structures 143 (2016) 272–286

design sensitivity and optimal design of composite thin-walled laminated beams using torsion-bending beam model. Liu [14,15] developed the analytical sensitivity analysis method for the composite structures based on the finite element method and lightweight design method using the analytical sensitivity. Neto et al. [16] performed the sensitivity analysis and optimal design of composite beam structures using finite element solver FEAP. Sedaghati et al. [17] developed a finite element model to study the mechanical and electrical behavior of laminated composite beam with piezoelectric actuators and a design optimization methodology by combining the finite element model and the sequential quadratic programming technique. Hamdaoui et al. [18] investigated an optimal design approach for choosing the most suitable material for high damping and low mass for a sandwich beam. Valido and Cardoso [19] implemented the optimal design of the various geometrically nonlinear composite laminate beam structures, which is based on finite element analysis and sensitivity analysis model. Kim et al. [20] designed and manufactured the hybrid glass/carbon composite bumper beam via the design optimization process combined with the impact analysis. Belingardi et al. [21] optimized beam section profile and beam curvature for crashworthiness using commercial finite element software ABAQUS. Blasques [22] developed a methodology for simultaneous topology and material optimization in optimal design of laminated composite beams with eigenfrequency constraints. Blasques and Stolpe [23] described a novel framework for simultaneous optimization of topology and laminate properties in structural design of laminated composite beam cross sections. The disadvantage of finite element analysis is that the approximate numerical solution and sensitivity obtained by the finite element models heavily depend on the mesh schedules [14]. On the other hand, the finite element approaches have much lower efficiency than the analytical approaches of the continuous composite beam models. Many researchers have developed the analytical approaches to achieve the exact responses and sensitivities for the continuous composite beam models. For example, the exact vibration frequencies of the continuous composite beams are achieved using the Euler–Bernoulli beam model [4–7] and Timoshenko beam model [8–10]. The analytical analysis methods and optimization design of the continuous composite models using the different non gradient-based algorithms (particle swarm algorithm and genetic algorithm) for the thin-walled composite box-beam helicopter rotor blades have been investigated [24–26]. Liu [5] derived the exact solutions and sensitivity of the first four frequencies using the continuous composite model and developed the gradientbased algorithm to achieve the lightweight design of the solid composite laminated beams. Lentz and Armanios [27] described a gradient-based optimization scheme for obtaining the maximum coupling in thin-walled composite beams subject to hygrothermal and frequency constraints. Roque and Martins [28] used differential evolution optimization to find the volume fraction that can maximize the first natural frequency for a functionally graded beam. The optimization methods for the composite beam structures, as mentioned above, can be classified into gradient-based and non gradient-based algorithms. The non gradient-based algorithms are also called random search algorithms. The random search algorithms can implement the optimization design without the gradient information. However, the gradient-based algorithms require the gradient to construct the searching algorithm. Therefore, the non gradient-based algorithms are easier to be carried out than the gradient-based algorithms. However, the random search algorithms have to implement a large number of structural reanalysis and may not find the optimum design if the sampling number is not big enough. Compared with the random search algorithms, the gradient-based algorithms are more efficient and can find the

optimum design (at least the local optimum) if the gradient can be computed efficiently and accurately [29–31]. Therefore, the sensitivity analysis is the key technique for the gradient-based algorithms. Today the works on the sensitivity analysis for composite laminated beams are almost limited to the sensitivity analysis of the frequency and deflection. A few work is reported on the sensitivity analysis of the stresses for composite laminated beams. The aim of this paper is to develop the novel optimization techniques for the lightweight designs of Timoshenko composite laminated beams (continuous model) subjected to the static loadings. The optimization model for the lightweight design is to find the width and depth (or layer thickness) of the solid composite beams to minimize the mass of the beams under the stiffness failure, strength failure and delamination failure constraints. The paper is arranged as follows. In Section 2, the exact analytical expressions of the displacements and stresses are derived using the Timoshenko composite continuous beam model. In Section 3, the exact analytical sensitivity formulae of deflections and stresses are achieved by direct differentiation. In Section 4, the main failure criteria for the composite laminated beams are analysed and employed. In Section 5, the lightweight design optimization model is formulated. In Section 6, the analytical gradients of objective function (mass) and constraint functions are achieved by using the sensitivity information of the deflections and stresses. In Section 7, the standard gradient-based nonlinear programming algorithms coded in Matlab, i.e., Sequential Quadratic Programming (SQP), interior-point penalty algorithm and active-set algorithm, are proposed to find the lightweight designs of Timoshenko composite laminated beams. In Section 8, the lightweight designs of the composite laminated beams with different boundary conditions (pinned-pinned, fixed-fixed, fixed-free and fixed-pinned) are performed using different optimization algorithms with the exact sensitivities proposed in this paper. The finite difference sensitivities are also compared with the exact sensitivities proposed in this paper. The lightweight designs of the composite beams obtained by active set method using both exact sensitivities and finite difference sensitivities are discussed. 2. Exact analytical displacements and stresses 2.1. Exact analytical displacements Fig. 1 is a segment dx of composite laminated beam with rectangular cross section. The number of layers are denoted by ð1Þ; ð2Þ; . . . ; ðNÞ and the fiber orientations of layers are denoted by hð1Þ , hð2Þ ; . . . ; hðNÞ . The layered positions are denoted by z1 ; z2 ; . . . ; zN

z

b

z N +1

zN

z2 z1

y

h 2

θ (N )

zk +1 zI

( N −1)

(N)

θ

( N − 1)

θ (k )

(k ) (I ) ( 2)

θ

(I )

θ ( 2)

θ (1)

(1) dx

Fig. 1. Composite laminated beam (segment dx).

x h 2

274

Q. Liu / Composite Structures 143 (2016) 272–286

q ( x)

and zNþ1 . The width and thickness of the rectangular cross section of the beams are denoted by b and h, respectively. Based on the first-order shear deformation theory (also called Timoshenko beam theory), the assumed displacement fields for the composite laminated beam are

uðx; zÞ ¼ u0 ðxÞ þ z/ðxÞ

ð1Þ

wðx; zÞ ¼ w0 ðxÞ

ð2Þ

where u0 ðxÞ and w0 ðxÞ are the displacements of the points on the middle plane along the x axis and z axis, respectively. /ðxÞ is the rotation of the normal to the middle plane about the y axis. uðx; zÞ and wðx; zÞ are the in-plane displacement and deflection of the composite beam. The strain and displacement relationships are given by the following equations.

ex ¼

du0 ðxÞ d/ðxÞ þz dx dx

cxz ¼ /ðxÞ þ

dw0 ðxÞ dx

ð3Þ

du0 ðxÞ d/ðxÞ þ bB11 Nx ¼ bA11 dx dx du0 ðxÞ d/ðxÞ þ bD11 dx dx
 dw0 ðxÞ Q x ¼ bKA55 /ðxÞ þ dx

M x ¼ bB11

N X  ðkÞ ðzkþ1  zk Þ Q 11

dx

dNx þ pðxÞ ¼ 0 dx

ð12Þ

dQ x þ qðxÞ ¼ 0 dx

ð13Þ

dMx  Q x þ mðxÞ ¼ 0 dx

ð14Þ

By substituting Eqs. (5)–(7) into Eqs. (12)–(14), we have 2

d u0 ðxÞ

d u0 ðxÞ 2

dx

ð6Þ

A55

2

dx


 dw0 ðxÞ þ mðxÞ ¼ 0  bKA55 /ðxÞ þ dx

2

d u0 ðxÞ 2

2

þ bB11

d /ðxÞ

¼0 2 dx dx " # 2 d/ðxÞ d w0 ðxÞ bKA55 þ þ qðxÞ ¼ 0 2 dx dx

bA11

2

d u0 ðxÞ

bB11

2

dx

2

þ bD11

d /ðxÞ 2

dx


 dw0 ðxÞ  bKA55 /ðxÞ þ ¼0 dx

ð18Þ

ð19Þ

ð20Þ

Eq. (18) can be rewritten as

ð8Þ

N 1X  ðkÞ ðz2  z2 Þ Q k 2 k¼1 11 kþ1

ð9Þ

N 1X  ðkÞ ðz3  z3 Þ ¼ Q k 3 k¼1 11 kþ1

ð10Þ

55

d /ðxÞ

ð16Þ

When only the transversal loading qðxÞ is applied to the composite laminated beam, Eqs. (15)–(17) can be simplified as

2

d u0 ðxÞ 2

N X  ðkÞ ðzkþ1  zk Þ ¼ Q

2

þ bD11

ð15Þ

ð17Þ

dx

D11

2

þ bB11

d /ðxÞ

þ pðxÞ ¼ 0 2 2 dx dx " # 2 d/ðxÞ d w0 ðxÞ þ qðxÞ ¼ 0 bKA55 þ 2 dx dx

bA11

2

k¼1

B11 ¼

Qx + dQx

Fig. 2. Free-body diagram.

bB11

ð7Þ

N x + dN x

m ( x)

ð5Þ

where Nx ,Mx and Q x are the normal force along the x axis, bending moment about the y axis and shear force along the z axis on the cross section of the composite laminated beam, respectively. The subscript x indicates the plane of action of the internal forces. b is the width of the beam. K is the shear correction factor, taken as 5/6 to account for the parabolic variation of transverse shear stress. The extensional stiffness A11 , bending stiffness D11 , bendingextensional coupling stiffness B11 and extensional stiffness A55 of the composite laminate are given by

A11 ¼

M x + dM x p ( x)

Nx

ð4Þ

where ex is normal strain and cxz is shear strain. All the other strain components (i.e., ey , ez , cxy and cyz ) are zero. The integrated constitutive equations can be written as

Qx

Mx

B11 d /ðxÞ A11 dx2

ð21Þ

By substituting Eq. (21) into Eq. (20), we have

D11 

ð11Þ

2

¼

! 2
 B211 d /ðxÞ dw0 ðxÞ ¼0  KA55 /ðxÞ þ 2 dx A11 dx

ð22Þ

Differentiating Eq. (22) with respect to x leads to

! " # 3 2 B211 d /ðxÞ d/ðxÞ d w0 ðxÞ þ b D11   bKA55 ¼0 3 2 dx A11 dx dx

ð23Þ

k¼1

where zk and zkþ1 are the position coordinates of the bottom and top surfaces of the kth lamina, as shown in Fig. 1. The stiffness coeffi ðkÞ of the kth lamina in the laminate coordinate system cients Q ij

can be computed by using Eq. (60) in Section 2.2. When a segment dx of composite laminated beam is subjected to the distributed normal force pðxÞ, transversal force qðxÞ and moment mðxÞ, as shown in Fig. 2. The equilibrium equations of a segment dx are as follows.

By substituting Eq. (19) into Eq. (23), we obtain 3

d /ðxÞ 3

dx

¼

A11 bðB211  A11 D11 Þ

qðxÞ

ð24Þ

Integrating Eq. (24) leads to 2

d /ðxÞ 2

dx

¼

A11 bðB211  A11 D11 Þ


Z

 qðxÞdx

ð25Þ

275

Q. Liu / Composite Structures 143 (2016) 272–286

R

where ½ qðxÞdx produces one indefinite integration constant C 1 . Thereafter the other formulae have the same integration constant R (i.e., C 1 ) when qðxÞdx is in a square bracket.

 u0 ðxÞ ¼ B

We integrate Eq. (25) and have

d/ðxÞ A11 ¼ dx bðB211  A11 D11 Þ


ZZ

 qðxÞdxdx

ð26Þ

RR

 where qðxÞdxdx produces two indefinite integration constants C 1 and C 2 . Thereafter the other formulae have the same integration RR constants (i.e., C 1 and C 2 ) when qðxÞdxdx is in a square bracket. Further integrating Eq. (26) leads to

/ðxÞ ¼


ZZZ

A11 bðB211  A11 D11 Þ

 qðxÞdxdxdx

ð27Þ

RRR where ½ qðxÞdxdxdx produces three indefinite integration constant C 1 , C 2 and C 3 . Thereafter the other formulae have the same RRR integration constants (i.e., C 1 , C 2 and C 3 ) when ½ qðxÞdxdxdxis in a square bracket. By substituting Eq. (25) into Eq. (21) and integrating, we have 2

d u0 ðxÞ 2

dx

¼


Z

B11 bðB211  A11 D11 Þ

du0 ðxÞ B11 ¼ dx bðB211  A11 D11 Þ u0 ðxÞ ¼ 

B11 bðB211

 A11 D11 Þ

 qðxÞdx

Z
Z

ZZ
Z

ð28Þ

 qðxÞdx dx

ð29Þ

 qðxÞdx dxdx

ð30Þ

The integral term in Eq. (30) produces three indefinite integration constants, i.e., C 1 , C 4 and C 5 , where C 1 is the same indefinite integration constant as that of Eq. (25), and produced by R ½ qðxÞdx . By substituting Eq. (26) into Eq. (19) and integrating, we obtain 2

d w0 ðxÞ 2

dx

1 A11 ¼ qðxÞ  2 bKA55 bðB11  A11 D11 Þ

 /ðxÞ ¼ A


ZZZ


ZZ

 qðxÞdxdx

ð31Þ


ZZ ) Z ( dw0 ðxÞ 1 A11 ¼  qðxÞ  qðxÞdxdx dx dx bKA55 bðB211  A11 D11 Þ

 qðxÞdxdxdx

ZZ
Z

ZZ  w0 ðxÞ ¼

ð37Þ

 qðxÞdx dxdx

  CqðxÞ A


ZZ

 qðxÞdxdx dxdx

ð38Þ

ð39Þ

The integral expressions of displacements and rotation in Eqs. (37)–(39) can be applied to any integrable loading forms. In this paper, we take the uniform distributed loading as an example to find the analytical expressions of the displacements and rotation. When qðxÞ ¼ q0 , Eqs. (37)–(39) lead to the following explicit formulae.

  q 0 x3 þ 1 C 1 x 2 þ C 2 x þ C 3 /ðxÞ ¼ A 2 6

ð40Þ

  q 0 x3 þ 1 C 1 x 2 þ C 4 x þ C 5 u0 ðxÞ ¼ B 2 6

ð41Þ

  1 x3  C  2 x2 þ C 6 x þ C 7  q0 x4  1 AC  q0 þ 1 AC w0 ðxÞ ¼ A 6 24 2 2

ð42Þ

The indefinite integration constants, i.e., C 1 , C 2 , C 3 , C 4 , C 5 , C 6 and C 7 , can be determined using the boundary conditions of the composite beams. The internal forces also need to be used to express the boundary conditions. The internal forces can be computed as follows. By substituting Eqs. (41) and (40) into Eq. (5), we have



 q0  bA11 B   bA11 B  q0 x2 þ bB11 A  C1x Nx ¼ bB11 A 2 3  2  bA11 BC  4 þ bB11 AC

ð43Þ

By substituting Eqs. (41) and (40) into Eq. (6), we obtain



 q0  bB11 B   bB11 B  q0 x2 þ bD11 A  C1x Mx ¼ bD11 A 2 3  2  bB11 BC  4 þ bD11 AC

ð44Þ

By substituting Eqs. (42) and (40) into Eq. (7), we have

 3 þ bKA55 C 6  bKA55 x þ bKA55 AC Q x ¼ Cq 0

ð45Þ

ð32Þ w0 ðxÞ ¼

ZZ ( 

1 A11 qðxÞ  bKA55 bðB211  A11 D11 Þ


ZZ

) qðxÞdxdx dxdx ð33Þ

The four indefinite integration constants of Eq. (33) are C 1 , C 2 , C 6 and C 7 , where C 1 and C 2 are the same indefinite integration conRR  stants as those of Eq. (26), and produced by qðxÞdxdx . We define the following three parameters so that the formulae become more concise and convenient.

¼ A ¼ B ¼ C

A11 bðB211  A11 D11 Þ B11 bðB211  A11 D11 Þ 1 bKA55

ð34Þ

ð35Þ

ð36Þ

The displacements, i.e., Eqs. (30) and (33), and rotation, i.e., Eq. (27) can be written as

2.2. Exact analytical stresses The nonzero strain components can also be achieved. By substituting Eqs. (41) and (40) into Eq. (3), we have

q 

 0 2  q 0 x2 þ C 1 x þ C 2 x þ C 1 x þ C 4 þ zA 3 2

ex ¼ B

ð46Þ

By substituting Eqs. (42) and (40) into Eq. (4), we obtain

 q0 3 1  1 x2  q0 x3  1 AC x þ C 1 x2 þ C 2 x þ C 3  A 2 2 6 6  2 x þ C 6  þ AC  Cq

cxz ¼ A

0

ð47Þ

The material coordinate system (123) and beam coordinate system/laminate coordinate system (xyz) are shown in Fig. 3. The angle between the x axis and 1 axis is hðkÞ (counter clockwise). The 1 axis is along the direction of the fiber. In this paper, the notations of stress and strain are the same as those in Jones’ text book [1]. Therefore, the plane stress components in the material coordinate system (123) are

276

Q. Liu / Composite Structures 143 (2016) 272–286

 ðkÞ ¼ Q cos4 hðkÞ þ 2ðQ þ 2Q Þsin2 hðkÞ cos2 hðkÞ þ Q sin4 hðkÞ Q 11 12 66 22 11

y 1

z=3

2

 ðkÞ ¼ ðQ þ Q  4Q Þsin2 hðkÞ cos2 hðkÞ þ Q ðsin4 hðkÞ þ cos4 hðkÞ Þ Q 11 22 66 12 12 x

θ (k )

0

 ðkÞ ¼ Q sin4 hðkÞ þ 2ðQ þ 2Q Þsin2 hðkÞ cos2 hðkÞ þ Q cos4 hðkÞ Q 11 12 66 22 22  ðkÞ ¼ ðQ  Q  2Q ÞsinhðkÞ cos3 hðkÞ Q 11 12 66 16 þðQ 12  Q 22 þ 2Q 66 Þsin hðkÞ cos hðkÞ 3

Fig. 3. The material coordinate system (123) and laminate coordinate system (xyz).

0

r1

1

0

B C ðkÞ  ðkÞ B @ r2 A ¼ T Q @

s12

 ðkÞ ¼ ðQ  Q  2Q Þsin3 hðkÞ coshðkÞ Q 11 12 66 26

1

ex ey C A; zkþ1 6 z 6 zk cxy

ð48Þ

þðQ 12  Q 22 þ 2Q 66 ÞsinhðkÞ cos3 hðkÞ  ðkÞ ¼ ðQ þ Q  2Q  2Q Þsin2 hðkÞ cos2 hðkÞ Q 11 22 12 66 66 þQ 66 ðsin hðkÞ þ cos4 hðkÞ Þ 4

ey = 0, cxy = 0, and

where

2

cos2 hðkÞ

6 T ðkÞ ¼ 4

sin hðkÞ

sin hðkÞ

cos2 hðkÞ

2

 sin h

ðkÞ

cos h

ðkÞ

sin h

ðkÞ

cos h

3

 ðkÞ ¼ Q cos2 hðkÞ þ Q sin2 hðkÞ Q 44 55 44

7 2 sin hðkÞ cos hðkÞ 5

 ðkÞ ¼ ðQ  Q Þ coshðkÞ sinhðkÞ Q 55 44 45

2 sin hðkÞ cos hðkÞ

2

ðkÞ

ðkÞ

cos2 h

2

 sin h

ðkÞ

ð49Þ 2

 ðkÞ Q 12  ðkÞ Q 22  ðkÞ Q 26

 ðkÞ Q 6 11 6Q  ðkÞ 4 12  ðkÞ Q 16

 ðkÞ ¼ Q

 ðkÞ Q 16  ðkÞ Q 26  ðkÞ Q 66

3 7 7 5

ð50Þ

1

r1 C rp1 ¼ B @ r2 A s12 0

ex

ð51Þ

1

C epx ¼ B @0A

ð52Þ

0 We rewrite Eq. (48) as

rp1 ¼ T ðkÞ Q ðkÞ epx ; zkþ1 6 z 6 zk

ð53Þ

The shear stress components in the material coordinate systems









s23  ðkÞ cyz ; zkþ1 6 z 6 zk ¼ T sðkÞ Q s cxz s13

where cyz = 0 and

"

T ðkÞ s ¼

sin hðkÞ cos h

"  ðkÞ ¼ Q s

 ðkÞ Q 44  ðkÞ Q 45

ðkÞ

cos hðkÞ

ð54Þ

cs ¼





s23 s13 0

where the stiffness coefficients Q ij of the lamina in material coordinate system are

E1 m12 E2 ; Q 12 ¼ ; 1  m12 m21 1  m12 m21 ¼ G12 ; Q 44 ¼ G23 ; Q 55 ¼ G13

#

 sin hðkÞ

where E1 and E2 are the modulus of elasticity, m12 and m21 are the Poisson’s ratio. G12 , G23 and G13 are the shear modulus, of the composite lamina. The indefinite integration constants are derived by using the different boundary conditions (pinned-pinned, fixed-fixed, fixedfree and fixed-pined), i.e., (1) Pinned-pinned The boundary conditions, i.e., u0 ð0Þ ¼ 0, w0 ð0Þ ¼ 0, N x ð0Þ ¼ 0, Mx ð0Þ ¼ 0, w0 ðLÞ ¼ 0, M x ðLÞ ¼ 0, Q x ð2LÞ ¼ 0, are used to determine the indefinite integration constants. The seven indefinite integration constants are

C1 ¼

 L   3D11 AÞq ð2B11 B 0 ;   6ðD11 A  B11 BÞ C 5 ¼ 0;

C6 ¼

C 2 ¼ 0;

C3 ¼

  B11 B  3D11 A 3   q0 L ; 72ðD11 A  B11 BÞ

   3D11 A B11 B  3 1   Aq0 L þ 2 Cq0 L; 72ðD11 A  B11 BÞ

C7 ¼ 0 ð62Þ

The boundary conditions, i.e., u0 ð0Þ ¼ 0, w0 ð0Þ ¼ 0, /ð0Þ ¼ 0, u0 ðLÞ ¼ 0, w0 ðLÞ ¼ 0, /ðLÞ ¼ 0, u0 2L ¼ 0, are employed to determine the indefinite integration constants. The seven indefinite integration constants are

ð57Þ

1 1 q L2 ; C 1 ¼  q0 L; C 2 ¼ 2 12 0 1 C 6 ¼ Cq L; C 7 ¼ 0 2 0

ð58Þ

C 3 ¼ 0;

C4 ¼

1 q L2 ; 12 0

C 5 ¼ 0; ð63Þ

(3) Fixed-pinned

So that we can rewrite Eq. (54) as

ð59Þ

 ðkÞ of the kth lamina in the laminate coorThe stiffness coefficients Q ij dinate system in Eqs. (50) and (56) are

ð61Þ

ð56Þ



ss ¼ T sðkÞ Q sðkÞ cs ; zkþ1 6 z 6 zk

E2 1  m12 m21

(2) Fixed-fixed

55

cxz

Q 22 ¼

ð55Þ

# ðkÞ

 Q 45  ðkÞ Q

Q 66

C 4 ¼ 0;

We also denote the shear stress and shear strain as

ss ¼

ð60Þ

Q 11 ¼

We use the following notations to depict the plane stress and strain in the material coordinate systems.

0

 ðkÞ ¼ Q sin2 hðkÞ þ Q cos2 hðkÞ Q 44 55 55

The boundary conditions, i.e., u0 ð0Þ ¼ 0, w0 ð0Þ ¼ 0, /ð0Þ ¼ 0, Þ ¼ 0, are used to determine w0 ðLÞ ¼ 0, N x ðLÞ ¼ 0, M x ðLÞ ¼ 0, Q x ð5L 8 the indefinite integration constants. The seven indefinite integration constants are

277

Q. Liu / Composite Structures 143 (2016) 272–286

 L2  3Cq  L2 þ 3Cq   5Aq Aq 0 0 0 0 ¼ ; C ; C 3 ¼ 0; 2   8AL 8A 2   7Aq0 L þ 9Cq0 5 ; C 5 ¼ 0; C 6 ¼ Cq L; C 7 ¼ 0 C4 ¼  8 0 24A

By substituting Eqs. (67) and (68) into Eqs. (8)–(11), we have

C1 ¼

ð64Þ

(4) Cantilever/Fixed-free

A11 ¼

 11 Lbq  3Lq  2B11 BLbq   2 bq 3AD 3L2 q0  B11 BL 0 0 0 0 ; C2 ¼   11 Þ ;   6bðB11 B  AD11 Þ 6bðB11 B  AD  11 L2 bq 3L2 q0  AD 0  C 3 ¼ 0; C 4 ¼  11 Þ ; C 5 ¼ 0; C 6 ¼ Cq0 L; C 7 ¼ 0   AD 6bðB11 B

C1 ¼

B11

D11

N 1X  ðkÞ Q ¼ 3 k¼1 11

ð65Þ

3.1. Option of design variables

ð66Þ

where d1 ¼ b and d2 ¼ h. 3.2. First derivatives of A11 , B11 , D11 , and A55 The deflection of the beam, i.e., Eq. (42), the strain components, i.e., Eqs. (46) and (47), and the indefinite integration constants, i.e., Eqs. (62)–(65), show it is necessary to obtain the first derivatives of A11 , B11 , D11 , and A55 so that the first derivatives of the deflection and stresses can be achieved. In this paper, the thickness of each lamina are supposed to be the same. The position coordinates of the bottom and top surfaces of the kth lamina can be expressed used the depth of the beam, i.e.,

h h 2k  N  2 h zk ¼  þ ðk  1Þ ¼ 2 N 2N h h 2k  N zkþ1 ¼  þ ðk þ 1  1Þ ¼ h 2 N 2N

" 3  3 # 2k  N 2k  N  2 h  h 2N 2N

ð71Þ

 N N X X  ðkÞ 2k  N h  2k  N  2 h ¼ h  ðkÞ Q Q 55 55 2N 2N N k¼1 k¼1

ð72Þ

The first derivatives of the stiffness coefficients A11 , B11 , D11 , and A55 with respect to design variables are obtained by direct differentiation.

ð73Þ

k¼1

The design variables of the composite laminated beam can be the fiber orientations, width and depth (thickness of layer) of the beam, or even the fiber volume fractions. Since the Timoshenko beam model in this paper is one dimensional beam, it is not necessary to optimize the fiber orientations because all the fiber orientations will go to zero to obtain the maximum bending stiffness. In this paper, the stiffness failure, strength failure and delamination failure are included in the optimization mathematical model to achieve the lightweight design of the composite laminated beams. As pointed out in the author’s previous work [5], the strength of the lamina will change when the fiber volume fraction of the lamina (layer) varies in the optimization process. Today the strength of the lamina is often determined by using experiment methods. It is still difficult to accurately predict the strength of the lamina using any mechanics theories [1], although the moduli of elasticity can be accurately calculated using micromechanics theory. Therefore, when the strength failure and delamination failure are included in the optimization mathematical model, it is unpractical to choose the fiber volume fractions as the design variables. Therefore, the width and depth of the beam are chosen as the design variables to achieve the lightweight designs of the composite laminated beams, i.e., T

ð70Þ

8 when di ¼ b > 0; @A11 < X ¼ 1 N  ðkÞ > @di Q 11 ; when di ¼ h :N

3. Exact analytical sensitivity of deflections and stresses

d ¼ ½b h

A55 ¼

ð69Þ

" 2  2 # 2k  N 2k  N  2 h  h 2N 2N

N 1X  ðkÞ ¼ Q 2 k¼1 11

The boundary conditions, i.e., u0 ð0Þ ¼ 0, w0 ð0Þ ¼ 0, /ð0Þ ¼ 0, Mx ð0Þ ¼  12 q0 L2 , N x ðLÞ ¼ 0, M x ðLÞ ¼ 0, Q x ðLÞ ¼ 0, are employed to determine the indefinite integration constants. The seven indefinite integration constants are

 N N X X  ðkÞ 2k  N  2k  N  2 h ¼ h  ðkÞ Q Q 11 11 2N 2N N k¼1 k¼1

ð67Þ

8 when di ¼ b > 0; @B11 < X N h i



¼  ðkÞ 2kN 2  2kN2 2 h; when di ¼ h > @di Q : 11 2N 2N k¼1

8 when di ¼ b > 0; @D11 < X N h i



¼  ðkÞ 2kN 3  2kN2 3 h2 ; when di ¼ h > Q @di : 11 2N 2N

ð75Þ

k¼1

8 when di ¼ b > 0; @A55 < X N ¼ 1  ðkÞ ; when di ¼ h >N Q @di : 55

ð76Þ

k¼1

3.3. First derivatives of deflections and stresses Since the first derivatives of A11 , B11 , D11 , and A55 with respect to design variables have been computed in Section 3.2, the first derivatives of deflections and stresses with respect to design variables can be obtained by direct differentiating Eqs. (42), (53) and (59). We have

   q @w0 ðxÞ @A 1 @A 0 4  @C 1 x3 ¼ C1 þ A x  @di @di 24 6 @di @di    @ C q0 1 @ A 1  @C 2 2 @C 6 @C 7 x þ  C2 þ A xþ þ @di 2 2 @di 2 @di @di @di  ðkÞ @ rp1 @T ðkÞ  ðkÞ @Q @e ¼ epx þ T ðkÞ Q ðkÞ px ; Q epx þ T ðkÞ @di @di @di @di

ð77Þ

zkþ1 6 z 6 zk ð78Þ

where the element

@ ex @di

in

@ epx @di

is

   q @ ex @B 0 2  @C 1 x þ @C 4 ¼ x þ C1x þ C4  B @di 3 @di @di @di    q @A @C @C 2 1 0 2  xþ þz x þ C 1 x þ C 2 þ zA @di 2 @di @di @e

The elements @dyi and

ð68Þ

ð74Þ

@ cxy @di

in

@ epx @di

ð79Þ ðkÞ

 ðkÞ

Q are equal to zero. @T@di and @@d can i

be obtained by direct differentiating Eqs. (49) and (50) with the design variable di .

278

Q. Liu / Composite Structures 143 (2016) 272–286

s @ ss @T sðkÞ  ðkÞ @Q @c ¼ c þ T sðkÞ Q sðkÞ s ; Q c þ T sðkÞ @di @di s s @di s @di ðkÞ

r1 r2 r2 s2 X c  X t Y  Yt  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ 2 þ 12 þ r1 þ c r2 ¼ 1 Xt Xc Xt Xc YtYc X t X c Y t Y c Y t Y c S2 r21

zkþ1 6 z 6 zk ð80Þ

where the element

@ cxz @di

in

@ cs @di

is

 q @ cxz @ A 1 0 3 ¼ x þ C 1 x2 þ C 2 x þ C 3 @di 6 2 @di   q  1 @C @C 2 @C 3 @A 1 @A 1 2 0 3  þA  x þ xþ C 1 x2 x  2 @di @di 6 2 @di @di @di    1  @C 1 2 @C @A  @C 2 x þ @C 6  A x  q0 þ C2 þ A 2 @di @di @di @di @di The element

@ cyz @di

in

@ cs @di

is equal to zero.

ðkÞ

@T s @di

and

ðkÞ

 @Q s @di

ð81Þ can be

obtained by direct differentiating Eqs. (55) and (56) with the design variable di . The analytical expressions of

 @B   @A , , @C @di @di @di

and

@C j @di

(j ¼ 1; 2;    7) in

Eqs. (77)–(81) can be derived using Eqs. (73)–(76).

where X t and X c are the tensile strength and compressive strength along the 1 axis of the material coordinate system, respectively. Y t and Y c are the tensile strength and compressive strength along the 2 axis of the material coordinate system, respectively. S is the shear strength on the plane 102 of the material coordinate system. The stress state (r1 , r2 and s12 ) are depicted in the material coordinate system 123, as shown in Fig. 3. 4.3. Stiffness failure criterion In the serviceability, if the deflection of the composite laminated beam is greater than the serviceability limit. The stiffness failure is supposed to happen. Therefore, the deflection of the composite laminated beam is not allowed to be greater than the serviceability limit. We have the following stiffness failure criterion.

 w0 6 w 4. Failure criteria of a composite laminated beam

Delamination, the separation of two adjacent layers in the composite laminates, is one of the most critical failure modes in composite laminated structures. Debonding between the adjacent layers depends on the stresses acting on that interface, i.e., the normal stress component r3 and the two shear stresses s12 and s23 . The researchers have developed many delamination failure criteria to predict the delamination failure of the composite laminated structures. In this paper, the delamination failure criterion by Yeh and Kim [32] is adopted to predict when the delamination will occur. The tensile delamination occurs if

r3

2

 þ

ZT

s13

2

 þ

S13

s23 S23

2 P 1;

when

r3 > 0

ð82Þ

where S12 and S23 are the shear strength on the planes 103 and 203 of the material coordinate system 123. Z T is the tensile strength along the 3 axis of the material coordinate system.The shear delamination occurs if



s13 S13

2

 þ

s23 S23

2 P 1;

when

r3 < 0

ð85Þ

 is the serviceability limit on the deflection of the composwhere w ite laminated beams.

4.1. Delamination failure criterion



ð84Þ

ð83Þ

In this paper, we let S13 ¼ S23 . For a composite laminated beam subjected to the transversal loading, the normal stress component on the plane 102 of the material coordinate system is less than zero, i.e., r3 < 0, therefore, the failure is shear delamination, not tensile delamination. 4.2. Tsai-Wu strength failure criterion Strength failure is also one of the most critical failure modes in composite laminated structures. Over the last five decades, there have been continuous efforts in developing strength failure criteria for composite laminated structures. Today, a large number of lamina failure criteria and laminate failure analysis methods have been developed, such as Tsai-Hill, Hoffman and Tsai-Wu failure criteria. In this paper, Tsai-Wu failure criterion is used to keep strength failure away from the design since Tsai-Wu strength failure criterion gives prediction that range from acceptable to excellent when it is compared with test results. The Tsai-Wu strength failure criterion has the following form [33].

5. Lightweight design optimization model Considering the stiffness failure criterion, strength failure criterion and delamination failure criterion, we have the lightweight design optimization model of composite beams as follows. T

Find d ¼ ½b h Minimize WðdÞ  r2 r2 s2 r1 r2 ffi Y c Y t t Subject to g j ¼ X t X1 c  pffiffiffiffiffiffiffiffiffiffiffiffiffi r þ r 1 <0 þ Y t Y2 c þ S122 þ XXctX 1 2 Xc Yt Yc fj¼



s213 S213

s223



Xt Xc Y t Y c

j

þ S2 1 < 0 23

j

0 6 0 r ¼ w0 ðaLÞ  w  b6b6b  h6h6h ð86Þ where WðdÞ is the mass of a composite laminated beam. aL is the location (x coordinate) where the deflection of the beam is monitored. a is different for different boundary conditions, i.e., PP: a ¼ 1=2, FF: a ¼ 1=2, FP: a ¼ 505=873, CL: a ¼ 1. b and b are the  are the lower and upper limits on the width of the beam. h and h lower and upper limits on the depth of the beam, respectively. In this paper, g, f and r are called strength failure function, delamination failure function and stiffness failure function, respectively. The subscript (j ¼ 1; 2; . . . ; N m :) means the jth monitored point of the strength and delamination, and Nm is the total number of the monitored points. It is noted that the analytical expressions of the strength failure function g, delamination failure function f and stiffness failure function r are naturally obtained by just substituting the analytical stresses (r1 , r2 , s12 , s23 and s13 ) and deflection w0 ðaLÞ in their right hand sides. When the composite laminated beams are subjected to the uniform distributed loading, the maximum of the shear force Q x is at x ¼ 0 or x ¼ L for PP, FF, FP and CL. The maximum of the bending moment M x is at x ¼ L=2 for PP and FF, at x ¼ 5L=8 for FP, at x ¼ 0 for CL. Therefore, we have the monitored points for strength failure and delamination failure, as shown in Fig. 4 for PP and FF, Fig. 5 for FP and Fig. 6 for CL. Double circles means two monitored points for the adjacent layers since their fiber orientations are not the same.

279

Q. Liu / Composite Structures 143 (2016) 272–286 Table 2 Material properties of lamina.

L 4

L 8

L 8

L 8

L 8

L 4

Fig. 4. Monitored points (Nm = 36) for strength failure and delamination failure (PP and FF).

Property

T300/5208

E1 (GPa) E2 ¼ E3 (GPa) G12 ¼ G13 (GPa) G23 (GPa) m12 ¼ m13

m23 q (kg/m3)

136.00 9.80 4.70 5.20 0.28 0.15 1540

X t (MPa) X c (MPa) Y t (MPa) Y c (MPa) S (MPa) S13 ¼ S23 (MPa)

1550 1090 59 59 75 75

Table 3 Deflection and sensitivity of deflection at monitored points (initial design).

L 4

L 8

L 8

L 8

L 8

L 4

Fig. 5. Monitored points (Nm = 36) for strength failure and delamination failure (FP).

Beam type (location)

PP (x ¼ L=2)

FF (x ¼ L=2)

FP (x ¼ 505L=873)

CL (x ¼ L)

w0 (m) @w0 =@b @w0 =@h

0.0145 0.0484 0.0863

0.0038 0.0126 0.0191

0.0068 0.0228 0.0375

0.1334 0.4445 0.8153

The first derivatives of the mass with respect to the design variables are

L 4

L 8

L 8

L 8

L 8

L 4

@W ¼ hLq @b

ð88Þ

@W ¼ bLq @h

ð89Þ

We have the analytical gradient of the objective functions as

½rWT ¼ Fig. 6. Monitored points (Nm = 12) for strength failure and delamination failure (CL).

 @W

 @W T @h

@b

ð90Þ

6.2. Analytical gradients of constraint functions The first derivatives of the strength failure function with respect to design variables can be achieved by direct differentiation.

Table 1 Main frame for nonlinear programming algorithms in Matlab. %% Choose optimization algorithm using gradient % ‘SQP’ or ‘interior-point’ or ‘active-set’ options = optimset(‘GradObj’, ‘on’, ‘GradConstr’, ‘on’, ‘Algorithm’, ‘sqp’); %% Implement optimization algorithms  @mycon, options); d = fmincon(@myfun, d , [], [], [], [], d, d, 0

% d0 is the initial design.

 @g 2r1 @ r1 1 @ r1 @ r2 2r2 @ r2 þ ¼  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ r1 @di X t X c @di @di @di Y t Y c @di Xt Xc Y t Y c 2s12 @ s12 X c  X t @ r1 Y c  Y t @ r2 þ 2 þ þ @di X t X c @di Y t Y c @di S

ð91Þ

We have the analytical gradient of strength failure function as

 are the lower and upper limits on the design variables. % d and d %% Computing objective and its gradient [mass, gradient_mass] = myfun(d); %% Computing constraints and their gradients [constraint, gradient_constraints] = mycon(d);

½rgT ¼

 @g

@b

 @g T @h

ð92Þ

The first derivatives of the delamination failure function with respect to design variables can also be achieved by direct differentiation.

6. Analytical gradients of objective function and constraint functions 6.1. Analytical gradient of objective function

@f 2s13 @ s13 2s23 @ s23 ¼ 2  2 @di S12 @di S23 @di

ð93Þ

The analytical gradient of delamination failure function is T

½rf  ¼

h

@f @b

@f @h

iT

ð94Þ

The mass of the composite laminated beam can be calculated as

W ¼ bhLq where q is the density of the composite material.

ð87Þ

7. Nonlinear programming algorithms Since the analytical expressions and analytical gradients of the objective function and constrained functions have been achieved,

280

Q. Liu / Composite Structures 143 (2016) 272–286

Table 4 Normal stresses and sensitivities of normal stresses (initial design). Beam type (location)

PP (x ¼ L=2, z ¼ h=2)

FF (x ¼ L=2, z ¼ h=2)

FP (x ¼ 5L=8, z ¼ h=2)

CL (x ¼ 0, z ¼ h=2)

r1 (MPa) @ r1 =@b (MPa.m1) @ r1 =@h (MPa.m1) r2 (MPa) @ r2 =@b (MPa.m1) @ r2 =@h (MPa.m1)

88.13 293.77 550.81 1.78 5.93 11.11

29.38 97.92 183.60 0.59 1.98 3.70

45.18 150.61 277.74 0.91 3.04 5.60

352.52 1175.10 2203.30 7.11 23.71 44.45

Table 5 Shear stresses and sensitivities of shear stresses (initial design) (at the bottom of the 4th layer). Beam type (location)

PP (x = 0, z = 0)

FF (x = 0, z = 0)

FP (x = 0, z = 0)

CL (x = 0, z = 0)

s12 (MPa) @ s12 =@b (MPa.m1) @ s12 =@h (MPa.m1) s23 (MPa) @ s23 =@b (MPa.m1) @ s23 =@h (MPa.m1) s13 (MPa) @ s13 =@b (m1) @ s13 =@h (m1)

0 0 0 2.01 6.71 4.20 2.23 7.43 4.64

0 0 0 2.01 6.71 4.20 2.23 7.43 4.65

0 0 0 2.52 8.39 5.25 2.79 9.29 5.80

0 0 0 4.03 13.43 8.39 4.46 14.86 9.29

Fig. 7. Pinned-pinned beam (initial design).

Q. Liu / Composite Structures 143 (2016) 272–286

281

Fig. 8. Fixed-fixed beam (initial design).

any standard gradient-based nonlinear programming algorithms can be used to solve the lightweight design optimization model, i.e., Eq. (86). In this paper, the standard gradient-based nonlinear programming algorithms coded in Matlab, i.e., Sequential Quadratic Programming (SQP), interior-point penalty algorithm and active-set algorithm, are used to find the lightweight designs of Timoshenko composite laminated beams. The main frame to implement the nonlinear programming algorithms in Matlab is shown in Table 1. 8. Examples The span of the composite laminated beams (PP, FF, FP and CL) are L = 7.2 m. The composite beams are subjected to the uniform distributed loading with intensity q0 = 105 N/m. The beams have the layered number N = 8. The material properties of the composite lamina is shown in Table 2. The given fiber orientations are [0/90/45/45]s in the examples. The initial design of the composite laminated beams are b = 0.3 m and h = 0.48 m (the thickness of each layer is 0.06 m), mass W = 1597 kg. The geometric design  = 2 m, h = 0.2 m and h  = 2 m. The serviceability space is b = 0.1 m, b

 = 0.01 m in Eq. (86). All the computations limit of the deflection w are implemented using the same computer (Processor: Intel(R) Core(TM) i5-2320 CPU @ 3.00 GHz, RAM: 8 GB). 8.1. Sensitivity of deflection and stresses The deflection and its first derivatives with respect to the width and depth of the composite laminated beam can be calculated using the exact analytical formulae Eqs. (42) and (77), respectively. The maximum deflections and their first derivatives with respect to design variables are shown in Table 3. The stresses in the material coordinate system can be computed using the exact analytical formulae Eqs. (48) and (54). The first derivatives of the stresses in the material coordinate system can be achieved using the exact analytical formulae Eqs. (78) and (80). The normal and shear stresses, their first derivatives with respect to design variables at the certain monitored points are listed in Tables 4 and 5. Compared with the stress, the deflection has more intuitively for the readers to see the deformation behavior of the composite laminated beams. All the deflections and their first derivatives with respect to the width and depth of the initial design with different

282

Q. Liu / Composite Structures 143 (2016) 272–286

Fig. 9. Fixed-pinned beam (initial design).

boundary conditions are plotted in Figs. 7–10. Since the allowable  = 0.01 m, we can see in Table 1 and Figs. 7(a) and 10 deflection is w (a) that the initial beams with pinned-pinned and fixed-free boundary conditions are infeasible design. It is noted that the initial design is not necessary to be feasible design for the three optimization methods, i.e., sequential quadratic programming, interior-point penalty algorithm and active-set algorithm, in this paper. The sensitivity information in Figs. 7(b), 8(b), 9(b) and 10 (b) indicates that, for the initial design, the depth is more sensitive to the deflection than the width of the beam. 8.2. Lightweight designs of composite beams In the section, the lightweight designs of Timoshenko composite laminated beams with pinned-pinned, fixed-fixed, fixed-pinned and fixed-free boundary conditions are achieved by using the optimization techniques proposed in this paper. The three standard gradient-based optimization methods coded in Matlab, i.e., Sequential Quadratic Programming method (SQP), interior-point penalty method (interior-point) and active set method (activeset), are employed to achieve the lightweight designs (optimal designs) of the composite laminated beams. The lightweight

designs are shown in Table 6 by SQP, Table 7 by interior-point and Table 8 by active set. The computational effort of the three optimization methods are given in Tables 6–8. The computational time to calculate the objective function, constrained functions and their gradients (one time) is shown in Table 9. The lightweight designs of the composite laminated beams with PP, FF, FP and CL boundary conditions in Tables 6–8 indicate that the three gradient-based optimization methods converge to the same lightweight design. Therefore, any one of the three gradient-based optimization methods can be employed to find the lightweight designs of the composite laminated beams. However, the iterations and function counts of the three optimization methods are different. The SQP use the least iteration and function count to achieve the lightweight designs. However, the CPU time depends not only on the iteration and function count, but also on the computational efficiency of objective function, constrained functions and their gradients. For example, the computational time of objective function, constrained functions and their gradients (one time) for FF beam, as shown in Table 9, is much less than those of PP, FP and CL beams. Even though the iteration and function count to achieve the lightweight designs of FF beam are

283

Q. Liu / Composite Structures 143 (2016) 272–286

Fig. 10. Fixed-free beam (initial design).

Table 6 Lightweight designs of Timoshenko composite laminated beams (SQP). Beam type

PP

FF

FP

CL

Optimal design

Mass W (kg) Width b (m) Depth h (m) Layer thickness h=N (m)

909.2634 0.1000 0.8200 0.1025

560.7427 0.1000 0.5057 0.0632

706.5145 0.1000 0.6372 0.0796

2065 0.1000 1.8623 0.2328

Computational effort

Iteration and function count CPU time (s)

7 7 45

8 9 6

6 7 30

9 11 32

Table 7 Lightweight designs of Timoshenko composite laminated beams (interior-point). Beam type

PP

FF

FP

CL

Optimal design

Mass W (kg) Width b (m) Depth h (m) Layer thickness h=N (m)

909.2634 0.1000 0.8200 0.1025

560.7427 0.1000 0.5057 0.0632

706.5145 0.1000 0.6372 0.0796

2065 0.1000 1.8623 0.2328

Computational effort

Iteration and function count CPU time (s)

11 17 60

13 32 11

12 16 38

11 13 21

284

Q. Liu / Composite Structures 143 (2016) 272–286

Table 8 Lightweight designs of Timoshenko composite laminated beams (active-set). Beam type

PP

FF

FP

CL

Optimal design

Mass W (kg) Width b (m) Depth h (m) Layer thickness h=N (m)

909.2634 0.1000 0.8200 0.1025

560.7427 0.1000 0.5057 0.0632

706.5145 0.1000 0.6372 0.0796

2065 0.1000 1.8623 0.2328

Computational effort

Iteration and function count CPU time (s)

7 15 52

24 163 56

18 74 174

9 17 27

Table 9 Computational time for objective function, constrained functions and their gradients. Beam type

PP

FF

FP

CL

Computational time (s)

3.5

0.3

2.4

1.6

lightweight designs, the stiffness failure is completely active constraints (equal to zero), the strength failure and delamination failure are not really active (not equal to zero). This numerical  increases, it is still results mean if the allowable deflection w possible to obtain the lighter weight designs of composite beams. 8.3. Comparison with central difference method

Table 10 Deflection and sensitivity of deflection (lightweight design). Beam type (location)

PP (x ¼ L=2)

FF (x ¼ L=2)

FP (x ¼ 505L=873)

CL (x ¼ L)

w0 (m) @w0 =@b @w0 =@h

0.0100 0.1000 0.0319

0.0100 0.1000 0.0470

0.0100 0.1000 0.0382

0.0100 0.1000 0.0125

greater than those of PP, FP and CL beams (as shown in Table 7 more obviously), the CPU time to achieve the lightweight design of FF beam is far less than those of PP, FP and CL beams. This indicate that the high efficient sensitivity and gradient computation methods can greatly improve the efficiency of the optimization methods. It is very important to the lightweight design of the large scale composite beam structure system. Therefore, the computational efficiency of the sensitivity and gradient is often called the bottle neck for the gradient-based optimization methods in structural optimization. The analytical sensitivity and gradient methods are often more efficient than the numerical sensitivity and gradient analysis methods. However, it is often tedious and of heavy work to derive the exact analytical expression formulae of sensitivity and gradient. It should be noted that all the exact analytical expression formulae of sensitivity and gradient for the composite laminated beams in this paper are derived by using the Symbolic Math Toolbox of Matlab R2015a. The maximum deflection and its first derivatives with respect to the width and depth of the lightweight design are listed in Table 10. The deflections, as shown in Table 10, indicate that the maximum deflection of the lightweight designs is equal to the allowable deflection. The sensitivity of the deflection in Table 10 indicates that, for the lightweight designs, the width is more sensitive to the deflection than the depth of the composite beams. That’s because the width of the lightweight design beams has reached the lower limit of the width. The numerical value/scope of the constrained functions of lightweight designs at all the monitored points are listed in Table 11. The numerical value/scope in Table 11 shows that, for the

In this section, the Central Difference Method (CDM) is employed to compute the first derivatives of the responses of Timoshenko composite laminated beams with respect to the width and depth. The central difference method has higher accuracy to obtain the sensitivity than both forward and backward difference methods. The sensitivity obtained by the central difference method is used in the gradient-based optimization methods in this paper. Due to space limitation, the sensitivity analysis and lightweight design of the composite beam with fixed-fixed boundary condition using CDM and active set method are implemented in this paper. The sensitivities of deflection with respect to width and depth are shown in Figs. 11 and 12, respectively. The lightweight designs of Timoshenko composite laminated beam are listed in Table 12. The results in Figs. 11 and 12 indicate that the sensitivities of deflection with respect to width and depth obtained by the central difference method are going to the exact solution by the proposed analytical method in this paper while the difference step size is getting smaller and smaller. The lightweight designs using exact derivatives and central difference derivatives (Table 12) show that the active set method using difference derivatives may not achieve the optimal design (the mass the optimum is 675.1210 kg using difference derivatives, 560.7427 kg using exact derivatives), although the computational efficiency of active set method using difference derivatives is equivalent to that of active set method using exact derivatives (CPU times are about 56 s). Therefore, it is necessary to develop the exact sensitivity analysis method instead of the difference methods for the gradient-based optimization algorithms. 9. Conclusion This paper developed the lightweight design methods for the composite laminated beams subjected to static loading. It should be noted that the proposed lightweight design methods are only applied to the continuous design variables. When the design variables are discrete, the evolutionary algorithms such as

Table 11 Numerical value/scope of constrained functions of lightweight designs at monitored points. Beam type

PP

FF

FP

CL

Stiffness failure function r Strength failure function g j

0.00000 0.97168 1.02101 0.99455 1.00000

0.00647 0.99137 1.01482 0.98193 1.00361

0.00000 0.97297 1.02987 0.98591 1.00000

0.00000 0.99433 1.01693 0.97867 0.99799

Delamination failure function f j

285

Q. Liu / Composite Structures 143 (2016) 272–286

Fig. 11. Sensitivity of deflection with respect to width.

Fig. 12. Sensitivity of deflection with respect to depth.

Table 12 Lightweight designs using exact derivatives and difference derivatives (active-set). Exact derivatives

Difference derivatives (db ¼ b=5; dh ¼ h=5)

Optimal design of composite beam (FF)

Mass W (kg) Width b (m) Depth h (m) Layer thickness (m)

560.7427 0.1000 0.5057 0.0632

675.1210 0.1371 0.4442 0.0555

Computational effort

Iteration and function count CPU time (s)

24 163 56

17 201 56

Sensitivity analysis using

Genetic Algorithms and Swarm Optimization Algorithms can be employed to solve optimization problems with discrete design variables. The main contributions and conclusions of this work are as follows. (1) The exact analytical formulae of displacements and stresses of the composite laminated beams with pinned-pinned, fixed-fixed, fixed-pinned and fixed-free boundary conditions are derived using Timoshenko beam theory (continuous beam model).

(2) The exact analytical sensitivity formulae of deflection and stresses, with respect to the width and depth (thickness of layer), of the composite laminated beams are derived using direct differentiation. (3) The lightweight design optimization model considering stiffness failure criterion, strength failure criterion and delamination failure criterion is proposed in this paper. (4) The analytical expression formulae for the mass, stiffness failure function, strength failure functions, delamination failure functions and their analytical gradients are derived

286

Q. Liu / Composite Structures 143 (2016) 272–286

using the exact analytical deflection, stresses and their analytical sensitivities with respect to the width and depth (thickness of layer) of the beams. (5) Three standard gradient-based optimization methods, i.e., sequential quadratic programming algorithm, interiorpoint penalty algorithm and active-set algorithm, coded in Matlab, are proposed to the lightweight designs of the composite laminated beams. The lightweight designs of the composite laminated beams with pinned-pinned, fixed-fixed, fixed-pinned and fixed-free boundary conditions are performed. (6) The three standard gradient-based optimization methods (SQP, interior-point and active-set) converge to the same lightweight design of the composite laminated beams. Any one of the three gradient-based optimization methods can be employed to find the lightweight designs of the composite laminated beams. The efficiency of the three optimization methods depends not only on the iteration and function count, but also on the computational efficiency of objective function, constrained functions and their gradients. The high efficient sensitivity and gradient computation methods can greatly improve the efficiency of the gradient-based optimization methods. (7) The optimization algorithms using the difference sensitivities may not lead to the optimal lightweight designs of the composite beams. It is necessary to develop the exact sensitivity analysis method instead of the difference methods for the gradient-based optimization algorithms.

Acknowledgements The author would like to thank the reviewers for their constructive comments. References [1] Jones RM. Mechanics of composite materials. 2nd ed. Philadelphia (PA): Taylor & Francis; 1999. [2] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. 2nd ed. Washington (DC): CRC Press; 2003. [3] Hale J. Boeing 787 from the ground up. boeing.com/commercial/aeromagazine 2006; 24(4): 17-23. . [4] Vinson JR, Sierakowski RL. Behavior of structures composed of composite materials. 2nd ed. Secaucus, NJ, USA: Kluwer Academic Publishers; 2002. [5] Liu Q. Analytical sensitivity analysis of eigenvalues and lightweight design of composite laminated beams. Compos Struct 2015;134:918–26. [6] Dancila DS, Armanios EA. The influence of coupling on the free vibration of anisotropic thin-walled closed-section beams. Int J Solids Struct 1998;35:3105–19. [7] Banerjee JR. Explicit analytical expressions for frequency equation and mode shape of composite beams. Int J Solids Struct 2001;38:2415–26.

[8] Chandrashekhara K, Krishnamurthy K, Roy S. Free vibration of composite beams including rotary inertia and shear deformation. Compos Struct 1990;14:269–79. [9] Chandrashekhara K, Krishnamurthy K, Wu WZB. Analytical solution to vibration of generally layered composite beams. J Sound Vib 1992;159 (1):85–99. [10] Banerjee JR. Frequency equation and mode shape formulae for composite Timoshenko beams. Compos Struct 2001;51:381–8. [11] Blasques JP, Stolpe M. Maximum stiffness and minimum weight optimization of laminated composite beams using continuous fibre angle. Struct Multidiscip Optim 2011;43:573–88. [12] Cardoso JB, Sousa LG, Castro JA, Valido AJ. Optimal design of laminated composite beam structures. Struct Multidiscip Optim 2001;24:205–11. [13] Cardoso JB, Valido AJ. Cross-section optimal design of composite laminated thin-walled beams. Comput Struct 2011;89:1069–76. [14] Liu Q. Analytical sensitivity analysis of frequencies and modes for composite laminated structures. Int J Mech Sci 2015;90:258–77. [15] Liu Q, Paavola J. Lightweight design of composite laminated structures with frequency constraint. Compos Struct 2015. http://dx.doi.org/10.1016/ j.compstruct.2015.08.116. [16] Netoa MA, Yu W, Leal RP. Generalized Timoshenko modelling of composite beam structures: sensitivity analysis and optimal design. Eng Optim 2008;40 (10):891–906. [17] Sedaghati R, Zabihollah A, Ahari M. Sensitivity analysis and optimal design of smart piezolaminated composite beams. AIAA J 2006;44(12):2987–96. [18] Hamdaoui M, Jrad M, Daya EM. Optimal design of frequency dependent threelayered rectangular composite beams for low mass and high damping. Compos Struct 2015;120:174–82. [19] Valido AJ, Cardoso JB. Geometrically nonlinear composite beam structures: optimization design. Eng Optim 2003;35(5):553–60. [20] Kim DH, Kim HG, Kim HS. Design optimization and manufacture of hybrid glass/carbon fiber reinforced composite bumper beam for automobile vehicle. Compos Struct 2015;131:742–52. [21] Belingardi G, Beyene AT, Koricho EG. Geometrical optimization of bumper beam profile made of pultruded composite by numerical simulation. Compos Struct 2013;102:217–25. [22] Blasques JP. Multi-material topology optimization of laminated composite beams with eigenfrequency constraints. Compos Struct 2014;111:45–55. [23] Blasques JP, Stolpe M. Multi-material topology optimization of laminated composite beam cross sections. Compos Struct 2012;94:3278–89. [24] Suresh S, Sujit PB, Rao AK. Particle swarm optimization approach for multiobjective composite box-beam design. Compos Struct 2007;81:598–605. [25] Murugan MS, Suresh S, Ganguli R, Mani V. Target vector optimization of composite box beam using real-coded genetic algorithm: a decomposition approach. Struct Multidiscip Optim 2007;33:131–46. [26] Kathiravan R, Ganguli R. Strength design of composite beam using gradient and particle swarm optimization. Compos Struct 2007;81:471–9. [27] Lentz WK, Armanios E. Optimum coupling in thin-walled, closed-section composite beams. J Aerosp Eng 1998;11:81–9. [28] Roque CMC, Martins PALS. Differential evolution for optimization of functionally graded beams. Compos Struct 2015;133:1191–7. [29] Liu Q, Paavola J, Zhang Jing. Shape and cross-section optimization of plane truss structures subjected to earthquake excitation using gradient and Hessian matrix calculations. Mech Adv Mater Struct 2016;23(2):156–69. [30] Ganguli R, Chopra I. Aeroelastic optimization of a helicopter rotor with composite coupling. J Aircr 1995;32(6):1326–34. [31] Ganguli R, Chopra I. Aeroelastic optimization of a helicopter rotor with twocell composite blades. AIAA J 1996;34(4):835–41. [32] Yeh HY, Kim CH. The Yeh-Stratton criterion for composite materials. J Compos Mater 1994;28(10):926–39. [33] Tsai SW, Wu EM. A general theory of strength for anisotropic materials. J Compos Mater 1971;5(1):58–80.