Buckling optimization of variable-stiffness composite panels based on flow field function

Accepted Manuscript Buckling Optimization of Variable-Stiffness Composite Panels based on Flow Field Function Peng Hao, Chen Liu, Xiaojie Yuan, Bo Wang, Gang Li, Tianyu Zhu, Fei Niu PII: DOI: Reference:

S0263-8223(17)32684-3 http://dx.doi.org/10.1016/j.compstruct.2017.08.081 COST 8839

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

25 July 2017 12 August 2017 17 August 2017

Please cite this article as: Hao, P., Liu, C., Yuan, X., Wang, B., Li, G., Zhu, T., Niu, F., Buckling Optimization of Variable-Stiffness Composite Panels based on Flow Field Function, Composite Structures (2017), doi: http:// dx.doi.org/10.1016/j.compstruct.2017.08.081

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Buckling Optimization of Variable-Stiffness Composite Panels based on Flow Field Function Peng Hao*, Chen Liu*, Xiaojie Yuan*, Bo Wang*#, Gang Li*, Tianyu Zhu*, Fei Niu** *

State Key Laboratory of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics,

International Research Center for Computational Mechanics, Dalian University of Technology, Dalian, 116023, China # **

Corresponding author: [email protected]

China Academy of Launch Vehicle Technology, Beijing, 100076, China

Due to the non-uniform in-plane stress distribution, variable-stiffness panel with curvilinear fiber paths is a promising structural concept for cutout reinforcement of composite structures under axial compression, due to the more diverse tailorability opportunities than simply choosing the best straight stacking sequence. However, traditional representation methods of curvilinear fiber path are usually not flexible for cutout reinforcement. In this study, the flow field function containing a uniform field and several vortex fields is utilized to represent the fiber path due to its inherent non-intersect and orthotropic features, and a bi-level optimization framework of variable-stiffness panels considering manufacturing constraints is then proposed. A typical rectangular composite panel with multiple cutouts is established to demonstrate the advantage of proposed framework by comparison with other fiber path functions. Results indicate that the flow fiber path only needs few variables to finely represent the fiber path, which can provide satisfying and manufacturable fiber paths by combination use of curvature constraint.

Keywords: variable-stiffness panels; cutout; buckling; fiber path; optimization; flow field function

I. .... Introduction Buckling is the main failure mode for composite thin-walled panels under axial load or combined load [1]-[9]. Cutouts are widely existed in various branches of thin-walled structures to accommodate the need of easy access,

1

inspection, electric lines, especially for launch vehicles, aircrafts, etc. When axial compression load is applied, buckling is usually the governing failure mechanism for these types of structures, and the buckling resistance would be reduced due to the presence of cutouts. Up to now, a large number of research activities have been devoted to investigate the influence of cutouts on the buckling behavior and load-carrying capacity of thin-walled structures [10]-[15]. It is found that the stress distribution is disturbed due to the material discontinuity caused by cutout, and an undesired stress concentration at the cutout vicinity is usually observed. When multiple cutouts are involved, local stress distribution would become extremely non-uniform [16]. Therefore, local stiffness tailoring near the cutout should be performed to improve the load-carrying mechanism. According to the literature review [16][17], much of the effort focuses on the increase of ply thickness near the cutout, however, this may lead to the fact that more loads are attracted to the region near the cutout, which has an unfavorable effect on the load-carrying capacity. Engels and Becker [18] found that the elliptical reinforcement can alleviate the stress concentration around the circular cutout. Besides, Nagendra et al. [19] performed the optimization of composite stiffened panels with a centrally located cutout. Compared to constant-stiffness designs (straight fiber path), a superior structural performance (e.g. buckling load, tension property) can be achieved for variable-stiffness designs, where fiber path is curvilinear and thus inplane stiffness is spatially varying throughout the structure [20]-[23]. In the previous studies [24], numerical results were provided to demonstrate that an improvement of buckling load ranging from 35% to 67% can be achieved by using a variable-stiffness design. The performance gain is attributed to the redistribution of in-plane loads to relatively stiff regions, and then resist buckling in critical regions. With the development of advanced manufacturing techniques for composite materials, Automated Fiber Placement (AFP) method makes it possible to fabricate laminates consisting of layers with curvilinear fibers [25]. As a significant insight into the cutout reinforcement of composite structures, Huang and Haftka [26] performed the optimization of fiber orientations near the cutout based on a piecewise bi-linear interpolation. For cutout reinforcement, variable-stiffness design is a good choice for tailoring local stiffness and global stiffness, which can improve the loading path to enhance the load-carrying efficiency. In the recent works by the authors [27][28], curvilinear stiffeners have been employed to reinforce the cutout, and an efficient bi-level optimization framework was established to achieve a high efficiency. As another type of variable-stiffness concept, functionally graded material (FGM) is also very promising, which can provide

2

varying mechanics properties along thickness direction, and the recent contributions regarding FGM can be found in [29]-[39]. For the design of variable-stiffness panels, a very important requirement is to guarantee the continuity and manufacturability of the distribution of fiber angles. Thus, the parameterization of fiber path and shift criterion within ply is of great significance. Liu and Haftka [40] introduced a new measure of continuity by distinguishing the composition continuity and stacking sequence continuity, and then a composite aircraft wing was designed based on this continuity constraint. As the simplest function for fiber path, linear variation fiber orientations were investigated by Gürdal et al. [41] in early times, which have already been widely used in the design of variable-stiffness panels. Besides, quadratic, cubical functions [42] and trigonometric functions [43] were also employed to describe the fiber paths. To enhance the design space, the fiber angles were then represented as cubic Bezier curves [44] and Lagrangian polynomials [45]. The principal stress trajectories were used to describe the family of fiber path by Zhu et al. [46]. Although the above representation method for fiber path definition is convenient for manufacturing, the potential of load-carrying capacity for variable-stiffness panels is greatly restricted. In order to ensure the continuity of fibers and keep a high level of design flexibility, a specific scalar function corresponding to the cocurrent and equipotential lines in the flow field was used to define a set of fiber angles within one ply [47]. The physical essence of stream and potential functions ensures that each line within one ply is non-intersect, which are particularly suitable for the description of fiber path. In addition, it should be noted that the stream and potential functions are a pair of conjugate functions, which means that the other function can be determined if one is defined. This feature is beneficial for the definition of orthogonal plies, because a pair of cocurrent and equipotential lines passing through arbitrary point are orthogonal to each other, which is also conducive for increasing the load-carrying capacity due to the mechanical balance. The result of tension experiment indicates that the strength of curvilinear laminate based on flow field is improved significantly. However, it should be noted that a pair of cocurrent and equipotential lines may not always lead to the optimum fiber path, but it is indeed suitable for substituting a series of manufacturing constraints (e.g. balance constraint, adjacent ply constraint, etc.). For the optimization of variable-stiffness panels, another challenging problem lies in the conflict between design space and optimization efficiency. Obviously, with the increase of formulation complexity for fiber path and shift criterion, the design space can be enlarged, however, the number of design variables also grows rapidly, and thus both the convergence rate and computational cost become very difficult to guarantee. Moreover, once the ply

3

number is relatively large, the optimization efficiency would be extremely low, and the optimization process would even be terminated prematurely. In this case, surrogate model is recommended to substitute the expensive finite element analysis, and a detailed comparison of different surrogate models was reported by Nik et al. [48] for the optimum design of variable-stiffness composite panels. Besides, cellular automata was adopted for the optimum design of fiber paths by Setoodeh et al. [49], where local rules are used to update both field and design variables to satisfy equilibrium and optimality conditions. Moreover, a bi-level framework was developed for the optimum design of variable-stiffness panels, and laminations parameters were used as design variables to limit the number of variables without sacrificing generality [50][51]. Peeters et al. [52] performed a multi-level optimization of variablestiffness laminates. In their works, the norm of the gradient of the fiber angle distribution is constrained to ensure manufacturability. Besides, Huang et al. [53] developed a new reanalysis assisted method to facilitate the optimization of variable-stiffness panels, and the contour lines of functions are employed to represent the fiber path. A multi-step optimization framework was established by Rouhi et al. [54], and its basic idea is to narrow the design space according to the optimum results in previous steps. Kiyono et al. proposed a novel fiber optimization method based on the normal distribution function, and fiber continuity is guaranteed by using a spatial filter [55]. Peeters and Abdalla [56] found that a significant improvement of buckling load can be achieved by combining ply drop and fiber angle optimizations, which provides some new insights for the design of variable-stiffness panels. Recently, Hao et al. [57] established the isogeomtric buckling analysis framework for variable-stiffness panels, and both high computational efficiency and prediction accuracy can be obtained compared to traditional FEA, which is promising for the efficient optimization of variable-stiffness panels. Montemurro and Catapano [58] presented a multi-scale two-level optimization framework for variable-stiffness panels, and all types of constraints are considered in the first level, which can be regarded as the improvement of the work by Wu et al. [50]. Additionally, a streamline analogy was used to represent the relationship between ply thickness and fiber angle variation, and the optimization of maximizing surface smoothness was performed to avoid the overlap of tow-placed courses [59]. In this study, variable-stiffness panels based on flow field is employed to meet the requirement of cutout reinforcement, since the stress distribution is highly non-uniform for the panel with cutout, and stiffness tailoring is significant for improving the stress distribution and load-carrying efficiency. The cocurrent and equipotential lines are used to parameterize each pair of adjacent plies. The paper is organized as follows. In Section 2, the fiber path based on flow field function as well as other typical representation functions of fiber paths are introduced. In Section

4

3, the optimization formulation for variable-stiffness panels based on flow field function is established, moreover, a bi-level optimization framework is proposed to accelerate the convergence rate. In Section 4, illustrative example is used to demonstrate the efficiency of optimization formulation for cutout reinforcement, and the influence of manufacturing constraint is investigated.

II. ... Representation of curvilinear fiber path 1.

Linear variation function According to the previous study [41], the in-plane analysis of rectangular panels under uniaxial compression

has been studied, which are balanced symmetric laminates using a linear variation of a pair of fiber angles [±θ(x)]s. The linear variation along the x axis, which is assumed to be centered along the length l of the panel, can be defined as follows

θ ( x) =

2(T1 − T0 ) x + T0 l

(1)

where T0 is the fiber orientation angle at the panel center, i.e. x = 0. T1 is the fiber orientation angle at the panel ends, i.e. x = ±l/2. By using such a simple function, closed-form solutions for the in-plane stress distribution and displacements of rectangular panels under uniform applied edge displacements are even possible. A generalization of this idea was to allow the direction of fiber orientation angle variation to be rotated with respect to the coordinate direction x, rather than limiting it to be along the x axis or the y axis. Hence, starting from an arbitrary reference point A with a fiber orientation angle T0 and moving along a direction x′ that is oriented by an angle φ from the coordinate axis x, the fiber orientation angle is assumed to reach a value T1 at a characteristic distance h from the reference point. Based on the linear variation of the fiber orientation angle between the points A and B, the equation for the fiber orientation angle along this reference path takes the following form

θ ( x′) = φ + (T1 − T0 )

x′ + T0 h

(2)

The fiber reference path associated with this unidirectional fiber orientation variation given in Eq. (2) is shown by the solid gray line in Fig. 1 (a). Although the fiber orientation angle is varying along a single axis, x′, assuming that the fiber orientation at other points in the domain is obtained by shifting this basic path in a direction

5

perpendicular to the x′ axis, the actual fiber orientation variation in the x–y plane is now a function of two coordinate directions, i.e. θ = θ (x, y). For traditional straight fiber path, each layer is represented by only a single angle. On the other hand, for variable-stiffness layer, three angles and a characteristic distance are used to represent a single layer. Assuming the characteristic distance to be associated with a geometric property of the part, the representation of a single curvilinear layer may be specified by φ. Three parameters may seem to be too restrictive for generating laminates with a wide range of tailoring possibilities. However, even if we chose the angles φ, T0 , and T1 from a small set of discrete orientations, a variety of laminate constructions are possible. For example, one can construct laminates with different combinations of the angle φ, laminates with ±φ angles, and stacks of layers. An example of a pair of laminates are shown in Fig. 1. (b).

(a) Definition of reference path;

(b) Variable-stiffness panel based on linear variation function.

Fig. 1 Fiber path based on linear variation function.

2.

Cubic polynomial function The fiber angles can also be represented in terms of cubic polynomials, which may represent the fiber contour

in the xy-coordinate system as

y ( x ) = ax 3 + bx 2 + cx + d

(3)

where a, b, c and d are constant. In order to make it easier to represent highly steep fiber contours everywhere in the xy-coordinate system, an auxiliary coordinate system shown in Fig. 2. (a) is utilized. Then, the fiber contours can be represented by

6

s ( r ) = ar 3 + br 2 + cr + d

(4)

where r = xcosα + ysinα. Therefore, for an arbitrary point in the xy-system, the fiber angle θ can be defined as

θ = arctan(3ar 2 + 2br + c ) + α

(5)

Since fibers are assumed to be straight within an element, the fiber angle for an element with mid-point coordinates (xm,ym), can be represented by

θ = arctan 3a ( xm cos α + ym sin α ) 2 + 2b ( xm cos α + ym sin α ) + c  + α

(6)

where a, b, c and α are the design variables. And then by the way of translation, the above function can be obtained through the entire design of the fiber path, as shown in Fig. 2 (b).

(a) Main and auxiliary coordinate systems; (b) Variable-stiffness panel based on cubic polynomial function. Fig. 2 Fiber path based on cubic polynomial function.

3.

Contour lines of cubic functions According to Ref. [53], the contour lines of cubic functions were used to represent fiber paths. Also, the cubic

functions can be expanded to other general functions. For example, the fiber paths shown in Fig. 3 are generated by the path functions z = x+a1y+a2xy+a3x2+a4y2+a5x2y+a6xy2+a7x3+a8y3. In order to analyze the mechanical properties, the fiber orientation angles need to be calculated according to path functions. Assume that the path function is

z = f ( x, y ) The contour line passing through an arbitrary point (x0, y0) can be expressed as 7

(7)

f ( x, y ) = x + a1 y + a2 xy + a3 x 2 + a4 y 2 + a5 x 2 y + a6 xy 2 + a7 x3 + a8 y 3

(8)

Therefore, the fiber orientation angle at (x, y) can be defined as

f x ( x, y )  arctan(− f ( x, y ) )  y θ ( x, y ) =  π   2

f y ( x, y ) ≠ 0 (9)

f y ( x, y ) = 0

Fig. 3 Variable-stiffness panel based on contour lines of cubic function. It should be noticed that the constant term of path function is not necessary to calculate fiber orientation angle. According to Eq. (8), we can obtain the fiber angle function as follow

 1 + a2 y + 2a3 x + 2a5 xy + a6 y 2 + 3a7 x 2 arctan( − )   a1 + a2 x + 2a4 y + a5 x 2 + 2a6 xy + 3a8 y 2 θ ( x, y ) =  π   2

4.

f y ( x, y ) ≠ 0 (10)

f y ( x, y ) = 0

Flow field function According to the works by Niu et al. [47], a basic potential flow can be decomposed into uniform linear flow

field, point source, point sink and point vortex. In this study, the flow field function is composed of one uniform linear flow field and two point vortexes based on superposition theorem. The flow function of uniform linear flow field can be expressed as

ψ 0 = U ⋅ ( y cos θ 0 − x sin θ 0 ) 8

(11)

where θ0 is the direction of the linear flow with respect to the x-axis of the local coordinate system, U is the velocity of linear flow. The flow function of point vortex is written as

ψi = −

Γi ⋅ ln ri + C1i 2π

(12)

where

ri =

2

( x − xi ) + ( y − yi )

2

(13)

where i expressed as ith point vortex. Γi is defined as the intensities of ith point vortex, and its sign indicates the point vortex flow direction. (xi, yi) is the coordinates of ith point vortexes. C1i is the regulating parameter of the vortex field function. Take the combination of uniform flow and two point vortexes as an example (i.e. i = 1, 2) as Fig. 4 (a)(b), the new flow function can be written as n

ψ = ∑ψ i = ψ 0 +ψ 1 +ψ 2 = U ⋅ ( y cos θ0 − x sin θ0 ) − 0

Γ1 Γ ⋅ ln r1 − 2 ⋅ ln r2 + C1 2π 2π

(14)

In order to simplify the function, several assumptions are made as follows

U = 1; C1 = C11 + C12 = 0 k = −Γ1 (2π ) = Γ 2 (2π ) x ∈ [ −l / 2 l / 2 ] ; y ∈ [ − w / 2 w / 2 ]

(15)

x1 = 0; y1 = 0.6w; x2 = 0; y2 = −0.6 w where k is the relative intensity of the point vortex. In the subsequent optimization problem, ki is substituted by Γi /2π to control the point vortex intensity. l is the length of panel, and w is the width of panel. Then, Eq. (14) can be given as

ψ = y cos θ 0 − x sin θ 0 − k ln r1 + k ln r2  2 2  r1 = ( x − x1 ) + ( y − y1 )  2 2 r2 = ( x − x1 ) + ( y − y1 ) 

(16)

According to the Cauchy-Riemann law, the conjugate potential function corresponding to the flow function can be expressed as

φ = y sin θ 0 + x cos θ 0 + k ⋅ arctan ( y − y1 , x − x1 ) − k ⋅ arctan ( y − y2 , x − x2 ) + C2

9

(17)

To pass through the origin (0, 0), the potential value has to be set as zero.

C2 = −k ⋅ [ arctan(−0.6w, 0) − arctan(0.6 w, 0)] = kπ

(18)

The fiber angle of flow function in (x, y) can be given as

 ∂ψ ∂ψ  ,   ∂x ∂y 

θ1 = arctan  −

where

∂ψ ∂x

and

∂ψ ∂x

(19)

are the partial derivatives of x and y, respectively. According to Eq. (15), we can obtain

k⋅x k⋅x  ∂ψ  ∂x = − sin θ0 − r 2 + r 2  1 2   ∂ψ = cos θ − k ( y − 0.6 w ) + k ( y + 0.6w ) 0  ∂y r12 r22

(20)

Similarly, the fiber angle of potential function in (x, y) can be obtained as

 ∂φ ∂φ  ,   ∂x ∂y 

θ 2 = arctan  −

(21)

According to Eq. (18), we can obtain the partial derivatives of x and y as follow

k ( y − 0.6 w ) k ( y + 0.6 w )  ∂φ +  = cos θ 0 − r12 r22  ∂x   ∂φ = sin θ + k ⋅ x − k ⋅ x 0  ∂y r12 r22

(22)

According to the potential flow theory, the cocurrent and equipotential lines passing through arbitrary point on the flow field are orthogonal to each other. For the flow field function, the fiber angle in a single layer is represented by⊥<θ0, k>, as shown in Fig. 4, in which the symbol ⊥ stands for a pair of layers described by cocurrent and equipotential lines. Obviously, the design flexibility is remarkably increased compared to linear variation fiber orientations and other representation methods, because local improvement of fiber path can be implemented by moving the point vortex or altering its intensity. Fig. 4 shows two examples using two flow field function based on two point vortexes and four point vortices to describe a layer respectively. In the FEM, the fiber orientation angle is discretized based on elements. Specifically, the fiber orientation angle in an element is treated as constant, and calculated using the centroid of the element.

10

y

x

(a) Fiber path based on flow function;

(b) Fiber path based on potential function; (with two point vortexes)

(c) Flow potential combination fiber path of rectangle panel with cutout (with four point vortexes). Fig. 4 Fiber path based on flow field function.

III. .... Optimization framework of variable-stiffness panels considering manufacturing constraints As already mentioned in the abstract, the AFP technology allows the designer to find more flexible fiber paths compared to straight fiber paths, and this significantly increases the load-carrying potential and design space of composite components. Since the fiber path based on flow field function is highly flexible and non-intersect naturally, moreover, it can be described by only limited governing variables, which can reduce the dimension of optimization problem. However, new types of constraints may be raised in the manufacturing process of variablestiffness panels, which should be integrated in the design optimization to guarantee the manufacturability of the

11

optimum design. In this section, two typical manufacturing constraints are introduced for the AFP process, which are described mathematically as follow The path function can be assumed as

z = f ( x, y )

(23)

The fiber path passing through an arbitrary point (x0, y0) can be expressed as

f ( x , y ) = f ( x0 , y 0 )

(24)

Therefore, the equation of the fiber path can be defined as

F ( x, y ) = f ( x, y ) − f ( x0 , y0 ) = 0

(25)

We assume y′ as follow

y′ = −

Fx ( x, y ) Fy ( x, y )

(26)

where, the subscript x and y indicates partial derivative. Once the curvature of an innermost tow of a course becomes too severe, the compressed side would exhibit local buckling or wrinkle modes [52][53]. For a variable-stiffness panel, the curvature of fiber path at arbitrary point can be expressed as

cur ( x, y ) =

y ′′ 1 + y′

3 2 2

where y′ can be calculated by Eq.(27), and y″ is the derivative of y′ on x .

Fig. 5 The selected points for constructing constraint.

12

(27)

Due to the fact that the influence of the flow field function on the curvature of the fiber path at the boundary is larger than the inner part, we only select the appropriate set of points in the finite region at the boundary to establish the constraints, as shown in Fig. 5.

Fig. 6 The proposed optimization framework of variable-stiffness panels based on flow field function. On this basis, a bi-level optimization framework is proposed to find the optimum design of variable-stiffness panels, and the basic idea is to split the original problem into several small-set problems, as shown in Fig. 6. In the first level, fiber path is imposed to be straight (i.e. uniform field of flow function), and fiber orientation angles are optimized with the objective of maximizing the buckling load, and a quasi-optimum design can be obtained. The optimization problem in the first level can be formulated as Maximum: λ l

u

Subject to: X i ≤ X i ≤ X i ,

i = 1, 2, 3, 4

(28)

where λ is the buckling factor, Xi is the ith design variable, X il and X iu are the lower and upper bounds of the ith design variable. In the second level, based on the optimum design in the first level, the flow field function containing a uniform field and several vortex fields is utilized to represent the fiber path. The optimization problem in the second level can be formulated as 13

Maximum: λ

cur ( x, y ) ≤ cm   l u Subject to:  X i ≤ X i ≤ X i , i = 1, 2, … , n  ( x, y ) ∈ Ω 

(29)

where Ω is the design domain, cm is the upper limit of fiber curvature.

IV. .... Optimization of composite panel with multiple cutouts 1. Model description In the following sections, a rectangle panel included two circular cutouts is investigated, as shown in Fig. 7, whose geometric dimensions are selected according the previous works by NASA [60][61]. To be specific, the length is 609.6 mm, and the width is 711.2mm. The center coordinate of large cutout is (-196, -86) mm, and the diameter is 107.5 mm. The center coordinate of small cutout is (-125, -243) mm, and the diameter is 53.8 mm. A 8ply balanced symmetric laminate is used. The lamina properties are set as E1=129.8 GPa, E2=9.2 GPa, G12=G13=5.1 GPa, G23=1.84 GPa, υ12=0.36, and the ply thickness is 0.1905 mm. For traditional method, thick rings are usually used to reinforce the stiffness loss caused by cutout. In this study, curvilinear fibers are employed to tailor the inplane stiffness, which can redistribute the loading path and improve the stiffness loss. For the purpose of illustration, the rectangular variable-stiffness panel with flow field fiber orientations are shown in Fig. 8.

Fig. 7 Sketch of composite panel with multiple cutouts.

14

Fig. 8 Sketch of composite panel with flow fiber path. 2. Loading condition and boundary condition For aircraft panels, combined compression-shear load is a common design condition, which requires more complex loading path to increase the resistance to buckling. This further highlights the load-carrying gain of variable-stiffness panels. Moreover, the combined loading condition may cause complicated mode shapes with coupling pattern for variable-stiffness panels. Also, boundary condition has a substantial influence on the buckling behavior and mode shape. Two typical boundary conditions are investigated in this study, including fully simply supported (SSSS), two loading edges simply supported and other two clamped (SSCC), as shown in Fig. 9.

(a) SSCC under combined compression-shear load; (b) SSSS under combined compression-shear load. Fig. 9 Two different boundary conditions of variable-stiffness panel.

15

3. Optimization results The proposed bi-level framework is used to conduct the optimization of variable-stiffness panels based on flow field function, and the optimization formulations in each step can be found in Section 3. Due to the constraint of symmetric laminates, only 4 plies need to be designed. In the first level, the optimization objective is to maximize the buckling factor λ, and the fiber angles of each ply are the design variable. The ranges of each variable θi are defined as follows

X i = θi ∈ [ −90°,90°]

i = 1, 2,3, 4

(30)

The stacking sequence of initial design is [0°, 45°, -45°, 90°]s. After performing linear buckling analysis, the buckling factors λ are 2.2160 and 1.3210 for the SSCC and SSSS boundary conditions, respectively. Then, the multi-island genetic algorithm (MIGA) is adopted to perform the first-level optimization. To be specific, the subpopulation size is selected as 100, and the number of islands is set as 5. The maximum generation number is set to be 20. The iteration histories for two boundary conditions are shown in Fig. 10. The optimum buckling factor λ of panel under SSCC boundary condition is 2.8928, with an improvement of 30.54%. The corresponding optimum stacking sequence is [45°, 60°, -70°, -80°]s. The optimum buckling factor λ of panel under SSSS boundary condition is 2.0297, with an improvement of 53.65%. The corresponding optimum stacking sequence is [45°, 50°, -45°, -80°]s, as shown in Fig. 11.

(a) SSCC;

(b) SSSS.

Fig. 10 Iterations history of buckling factor for two boundary conditions.

16

(a) SSCC;

(b) SSSS.

Fig. 11 Optimum stacking sequences of first-level optimization for two boundary conditions. In the second level, the representation method of fiber path is selected as the flow field function. Thus, the involved design variables include the position of each point vortex (xi, yi ), the relative intensity of the vortex field ki and the orientation of uniform field θ. To lift the design space, four point vortexes and uniform flow are used to describe the flow field. In order to better explore the advantages of four point vortexes, the point vortex is assumed to be distributed around panel. Therefore, the ranges of each design variable are defined as follows

x1 ∈ [ 304.8,914.4]

y j ∈ [ −711.2, 711.2]

x2 ∈ [ −914.4, −304.8]

y3 ∈ [ 355.6, 711.2]

xi ∈ [ −609.6, 609.6] km ∈ [ −200, 200]

y4 ∈ [ −711.2, −355.6]

i = 3, 4

j = 1, 2 (31)

θ ∈ [ −90°, 90°]

m = 1, 2, 3, 4

Besides, we normalize the ranges of each variable in the optimization, in order to obtain a smoother iteration history. 1) Optimization result of variable-stiffness panels without manufacturing constraints In the second level, the buckling optimizations of variable-stiffness panels under two boundary conditions are performed respectively. For this example, there are four point vortices and uniform flow to control the fiber angle in arbitrary domain. Similarly, only 4 plies need to be designed due to the constraint of symmetric laminates. Because of the orthogonality of the flow layer and the potential layer, a set of design variables can control two adjacent layers Ply-(1)×(2) or Ply-(3)×(4) simultaneously. Therefore, only the fiber angles in 2 plies are independent variables, and thus a total number of 26 variables are involved in the optimization. Due to the inherent feature of multimodal search problems, MIGA is used to find the optimum stacking sequence. The iteration history is also shown in Fig.

17

10. The optimum values of design variables in each ply are listed in Table 1. During the optimization process, the average CPU time of buckling analysis is about 4 min, using a PC with a CPU of Intel Xeon E5-2697 2.7 GHz and 128G RAM. Table 1 Optimum variable values of flow fiber path without manufacturing constraints. SSCC

Type

SSSS

Ply-(1)×(2)

Ply-(3)×(4)

Ply-(1)×(2)

Ply-(3)×(4)

x1

316.53

487.04

400.04

500.11

x2

-650.68

-446.61

-564.55

-437.66

x3

255.55

-379.84

-395.97

-272.57

x4

-307.36

382.35

-173.97

-64.881

y1

-488.99

-79.427

-399.13

-454.72

y2

35.367

-107.93

211.33

-661.52

y3

467.02

407.07

774.38

401.50

y4

-428.11

-390.11

-498.42

-727.46

k1

-162.17

95.006

27.726

-140.70

k2

131.57

108.85

74.308

-138.79

k3

97.789

190.80

-171.64

-190.19

k4

198.17

134.81

-197.91

186.99

θ

62.940

24.430

30.887

-8.8812

λ

3.4360

2.4504

Furthermore, the first two order buckling modes of optimum designs for two boundary conditions are listed in Table 2. It can be found that the buckling deformation mainly occurs in the vicinity of cutouts for the first-order mode. The difference of two boundary conditions leads to a distinction in the second-order mode. The buckling pattern indicates that the optimum fiber path can compensate the stiffness loss caused by cutouts, since the buckling waves are uniformly distributed on the panel in general. For the SSSS boundary condition, the first two order buckling loads only have a slight difference, but the buckling mode switches and changes significantly, which can be regarded as a symbol of high load-carrying efficiency.

18

Table 2 Buckling modes of variable-stiffness panels based on flow field function without manufacturing constraints.

Mode order

SSCC

SSSS

λ =3.4360

λ =2.4504

λ =4.2411

λ =2.4514

Legend

①

②

2) Optimization of variable-stiffness panels considering manufacturing constraints In the above optimization, it is found that the flow fiber path is more flexible, and results indicate that the loadcarrying capacity of variable-stiffness panel can be remarkably improved for a specific situation. However, the fiber curvature is exceedingly varied in some regions, and fiber enrichment can also be found. For example, when the vortex is too close to the panel, the phenomenon of fiber enrichment point can be found, as shown in Fig. 12. The existing manufacturing technologies such as ATL (Automated Tape Laying) and AFP (Automated Fiber Placement) are difficult to satisfy the requirements of this fiber path.

19

Exceedingly large fiber curvature or fiber enrichment region Fig. 12 Optimum flow fiber path without manufacturing constraints (SSCC). In order to ensure the optimum design to meet the practical manufacturing requirements, the fiber curvature constraint is introduced as a manufacturing constraint in the optimization process. The definition of fiber curvature has been given before as Eq. (27). In this study, the upper limit of fiber curvature cm is set to be 0.1, i.e. cur(x, y)≤0.1, which is identical with Ref. [53]. On this basis, the optimization of curvilinear fiber path is carried out. Taking the SSCC boundary condition for example, the iteration history of buckling factor λ is shown in Fig. 13. The optimum values of design variables are listed in Table 3.

Fig. 13 Iteration histories of buckling factor for straight and flow fiber paths (SSCC).

20

Table 3 Optimum variable values of flow fiber path (SSCC).

Type

Without manufacturing constraints

With manufacturing constraints

Ply-(1)×(2)

Ply-(3)×(4)

Ply-(1)×(2)

Ply-(3)×(4)

x1

316.53

487.04

852.73

-408.55

x2

-650.68

-446.61

451.85

-648.87

x3

255.55

-379.84

-176.89

101.33

x4

-307.36

382.35

168.26

-214.81

y1

-488.99

-79.427

-666.29

-280.11

y2

35.367

-107.93

4.1781

-16.181

y3

467.02

407.07

569.35

-658.67

y4

-428.11

-390.11

510.23

-643.28

k1

-162.17

95.006

52.067

-170.31

k2

131.57

108.85

144.88

56.510

k3

97.789

190.80

-20.291

-31.803

k4

198.17

134.81

-155.36

-74.961

θ

62.940

24.430

37.328

20.090

λ

3.4360

3.1907

After 25 iterations, the convergence of buckling factor is achieved, and the history curve exhibits a smooth manner during the optimization. By comparing the first-order mode of straight fiber path and curvilinear fiber path, the buckling factor λ of straight fiber path increases from 2.2160 to 2.8928, with an improvement of 30.54%. The buckling factor λ of curvilinear fiber path without manufacturing constraints increases from 2.2160 to 3.4360, with an improvement of 55.05%. While the buckling factor λ of curvilinear fiber path with manufacturing constraints increases from 2.2160 to 3.1907, with an improvement of 43.98%. The advantage of curvilinear fiber path can be attributed to the varying the fiber angle of governing variables, by which in-plane stiffness can be tailored to accommodate the combined compression-shear load. Furthermore, the comparison of buckling modes for different optimum designs is provided in Table 4. As is evident, two axial buckling waves can be found for the straight fiber path, and the one near the cutouts is distinctly larger than the one away from the cutouts. Due to the tailoring of in-plane stiffness distribution, the buckling waves

21

become more uniform for curvilinear fiber path, and the second-order buckling mode even splits into three axial waves, indicating a higher load-carrying efficiency. Table 4 Buckling modes of different optimum designs based on straight and flow fiber paths (SSCC). Mode order

Straight fiber path

Flow fiber path without manufacturing constraints

Flow fiber path with manufacturing constraints

λ =2.8928

λ =3.4360

λ =3.1907

λ =3.3191

λ =4.2411

λ =3.8380

Legend

①

②

Ply-1

Ply-2

Ply-3

Fig. 14 Optimum flow fiber path with manufacturing constraints (SSCC).

22

Ply-4

To further investigate the influence of manufacturing constraints, the optimum fiber path considering curvature constraint is shown in Fig. 14 for the SSCC boundary condition. As can be observed, the phenomenon of fiber enrichment is well eliminated, and the curvilinear fiber path is smoothly varied. In general, the arrangement of fiber satisfies the requirement of stiffness distribution to avoid stress concentration. This will allow for a better understanding of the load redistribution mechanism responsible for an increased buckling load. Moreover, the AFP technology is adequate for manufacturing this type of composite panels, and this can be attributed to the curvature constraint and inherent constraint of cocurrent and equipotential lines.

3) Comparison of other fiber path functions In this section, variable-stiffness design based on other different fiber path functions considering manufacturing constraints are compared, including linear variation function, cubic polynomial function, contour lines of cubic function. For the linear variation function, T0, T1 and φ in Eq. (2) are set as the design variables of each ply (i.e. 3 variables for each ply). For the cubic polynomial function, a, b, c and α in Eq. (6) are employed as design variables (i.e. 4 variables for each ply). For the contour lines of cubic functions, ai (i=1, 2, …, 8) in Eq. (10) are used as design variables (i.e. 8 variables for each ply). Similar to the cocurrent and equipotential lines in flow field function, a symmetrical assumption of two adjacent plies is made for the linear variation function and cubic polynomial function according to Refs. [40][44], in order to reduce the number of independent variables. As a comparison, each ply is individually designed for the contour lines of cubic function.

23

Fig. 15 Iterations history of buckling factor for different fiber path functions (SSSS). Table 5 Optimum variable values for different fiber path functions (SSSS). Linear variation function

Ply number

T0

T1

φ

Ply-(1)×(2)

19.870

15.962

25.795

Ply-(3)×(4)

45.796

89.404

-11.721

Cubic polynomial function

Ply number

a

b

c

α

Ply-(1)×(2)

6.0415E-04

2.0546E-02

4.8637E+00

-4.1272E+01

Ply-(3)×(4)

4.8844E-05

8.4733E-03

-1.9400E+00

1.7692E+01

Contour lines of cubic function

Ply number

a1

a2

a3

a4

a5

a6

a7

a8

Ply-(1)

40.901

22.474

-7.1435

-49.419

-37.149

37.904

22.161

-15.222

Ply-(2)

7.1770

29.612

40.414

-23.458

-31.212

49.854

30.732

-2.8862

Ply-(3)

-19.987

44.684

0.1030

-13.343

20.658

47.722

0.6416

21.899

Ply-(4)

48.306

-1.8395

25.639

-26.980

23.648

43.715

1.3039

23.545

Flow field function

Ply number

x1

x2

x3

x4

y1

y2

y3

Ply-(1)×(2)

386.07

-855.85

-213.51

268.50

555.25

-108.90

714.41

Ply-(3)×(4)

821.11

-354.25

411.25

308.03

134.95

138.64

573.86

24

y4

k1

k2

k3

k4

θ

Ply-(1)×(2)

-594.06

-95.451

118.76

-169.73

-80.119

37.446

Ply-(3)×(4)

-370.21

125.19

-125.96

-131.98

-149.22

-47.916

Taking the SSSS boundary condition for example, MIGA is employed to perform the optimization of variablestiffness panels based on the above fiber path functions. The algorithm parameters are identical with the previous ones. Besides, it should be stated that the curvature constraint is considered for each optimization, with a maximum allowable value of 0.1. The iteration histories of buckling factor λ are shown in Fig. 15 for different fiber path functions. It can be found that the convergence is achieved finally for each curve, and the history curve becomes very smooth after 30 iterations. By comparing the results, we can find that the performance of buckling resistance for the composite panel based on flow field function is stronger than those of other conventional fiber path functions. For each curvilinear fiber path, the buckling load is significantly higher than the one of optimum straight fiber path. The contour lines of cubic function and linear variation function yield the optimum fiber paths with similar buckling loads, and the cubic polynomial function produces the lowest buckling load. To be specific, the optimum values of design variables are listed in Table 5. Table 6 Buckling modes of different optimum designs based on various fiber path functions (SSSS). Fiber path function

Mode order ①

Mode order ②

λ =2.3391

λ =2.5269

Linear variation function

25

Legend

Cubic polynomial function

λ =2.3000

λ =2.3476

λ =2.3512

λ =2.5133

λ =2.4329

λ =2.4526

Contour lines of cubic function

Flow field function

Similarly, the buckling modes of different optimum designs based on various fiber path functions are provided in Table 6. As can be seen, the tendencies of buckling deformation for linear variation function and contour lines of cubic function are almost identical, and the buckling waves are not uniformly distributed on the panel. The ones for cubic polynomial function and flow field function are similar, while the buckling waves are more uniform for flow field function. To explain this phenomenon, the optimum fiber paths are also provided for different functions, as shown in Fig. 16. Corresponding to the highest buckling load and simultaneous buckling pattern, the fiber path described by flow field function is smooth and compatible due to its inherent feature, especially for the region near cutouts, which significantly improve the loading path and thus the buckling behavior. The AFP technology is

26

adequate for manufacturing this type of fiber path. While for other fiber path functions, although curvature constraints are considered in the optimization process, the fiber paths are still not satisfying for manufacturability, since they appear to be discontinuous and curious. For the contour lines of cubic function, even if 32 variables are involved to represent the fiber path, it still lacks the local variation capacity of fiber orientation angles to compensate the stiffness loss caused by cutouts, because all fiber paths need to be shifted from a reference path. By contrast, the flow field function only needs few variables to finely describe the fiber path (uniform linear flow field for global orientation and point vortex for local variation) and enhance the design flexibility, moreover, it includes the inherent non-intersect and orthotropic features, which can provide satisfying and manufacturable fiber paths by combination use of curvature constraint.

27

Fig. 16 Optimum fiber paths for different functions (SSSS).

28

V. .... Concluding remarks Thin-walled components with cutouts are very common in current aerospace structures, however, there is still no efficient reinforcement method to compensate for the stiffness loss caused by cutouts. Variable-stiffness panel with curvilinear fibers is a promising structural concept compared to constant-stiffness designs, especially for cutout reinforcement of composite structures under axial compression. Traditional curvilinear fiber path functions lack the local variation capacity of fiber orientation angles to compensate the stiffness loss caused by cutouts. In this study, the flow field function containing a uniform linear flow field and several vortex fields is utilized to represent the global orientation and local variation of fiber path, which can enhance the design flexibility with only few variables. Then, a bi-level optimization framework of variable-stiffness panels is developed for the flow field function, aiming to increase the design efficiency of fiber paths. A typical rectangular panel with multiple cutouts is established to demonstrate the advantage of flow field function and proposed framework. The buckling modes and fiber paths of obtained optimum designs are examined in detail. Also, the effects of boundary condition and manufacturing constraint are investigated. By comparison with other fiber path functions, including linear variation function, cubic polynomial function, contour lines of cubic function, the flow fiber path only needs few variables to finely describe the fiber path, which can provide satisfying and manufacturable fiber paths by combination use of curvature constraint.

Acknowledgments This work was supported by the National Natural Science Foundation of China (11402049 and 11372062), the National Basic Research Program of China (2014CB049000).

References [1]. Kolahchi R, Hosseini H, Esmailpour M. Differential cubature and quadrature-Bolotin methods for dynamic stability of embedded piezoelectric nanoplates based on visco-nonlocal-piezoelasticity theories. Compos Struct 2016; 157: 174-186. [2]. Meziane MAA, Abdelaziz HH, Tounsi A. An efficient and simple refined theory for buckling and free vibration of exponentially graded sandwich plates under various boundary conditions. J Sandw Struct Mater 2014; 16(3): 293-318.

29

[3]. Kolahchi R, Safari M, Esmailpour M. Dynamic stability analysis of temperature-dependent functionally graded CNTreinforced visco-plates resting on orthotropic elastomeric medium. Compos Struct 2016; 150: 255-265. [4]. Tounsi A, Houari MSA, Bessaim A. A new 3-unknowns non-polynomial plate theory for buckling and vibration of functionally graded sandwich plate. Struct Eng Mech 2016; 60(4): 547-565. [5]. Bourada M, Amara K, Tounsi A. Buckling analysis of isotropic and orthotropic plates using a novel four variable refined plate theory. Steel Compos Struct 2016; 21(6): 1287-1306. [6]. Chikh A, Tounsi A, Hebali H, et al. Thermal buckling analysis of cross-ply laminated plates using a simplified HSDT. Smart Struct Syst 2017; 19(3): 289-297. [7]. Mahi A, Tounsi A. A new hyperbolic shear deformation theory for bending and free vibration analysis of isotropic, functionally graded, sandwich and laminated composite plates. Appl Math Model 2015; 39(9): 2489-2508. [8]. Meziane MAA, Abdelaziz HH, Tounsi A. An efficient and simple refined theory for buckling and free vibration of exponentially graded sandwich plates under various boundary conditions. J Sandw Struct Mater 2014; 16(3): 293-318. [9]. Bourada M, Tounsi A, Houari MSA, et al. A new four-variable refined plate theory for thermal buckling analysis of functionally graded sandwich plates. J Sandw Struct Mater 2012; 14(1): 5-33. [10]. Starnes Jr JH. Effect of a slot on the buckling load of a cylindrical shell with a circular cutout. AIAA J 1972; 10(2): 227229. [11]. Kapania RK, Haryadi SG, Haftka RT. Global/local analysis of composite plates with cutouts. Comput Mech 1997; 19(5): 386-396. [12]. Priyadharshani S A, Prasad A M, Sundaravadivelu R. Analysis of GFRP stiffened composite plates with rectangular cutout. Compos Struct 2017; 169: 42-51. [13]. Li X, Gao W, Liu W. Post-buckling progressive damage of CFRP laminates with a large-sized elliptical cutout subjected to shear loading. Compos Struct 2015; 128: 313-321. [14]. Hao P, Wang B, Tian K, et al. Integrated optimization of hybrid-stiffness stiffened shells based on sub-panel elements. Thin-Walled Struct 2016; 103: 171-182. [15]. Dimopoulos CA, Gantes CJ. Numerical methods for the design of cylindrical steel shells with unreinforced or reinforced cutouts. Thin-Walled Struct 2015; 96: 11-28. [16]. Hao P, Wang B, Tian K, et al. Simultaneous buckling design of stiffened shells with multiple cutouts. Eng Optimiz 2017; 49(7): 1116-1132. [17]. Guo SJ. Stress concentration and buckling behaviour of shear loaded composite panels with reinforced cutouts. Compos Struct 2007; 80(1): 1-9. [18]. Engels H, Becker W. Optimization of hole reinforcements by doublers. Struct Multidiscip O 2000; 20(1): 57-66.

30

[19]. Nagendra S, Haftka RT, Gürdal Z, et al. Design of a blade-stiffened composite panel with a hole. Compos Struct 1991; 18(3): 195-219. [20]. Ghiasi H, Fayazbakhsh K, Pasini D, et al. Optimum stacking sequence design of composite materials Part II: Variable stiffness design. Compos Struct 2010; 93(1): 1-13. [21]. Rouhi M, Ghayoor H, Hoa SV, et al. Stiffness tailoring of elliptical composite cylinders for axial buckling performance. Compos Struct 2016; 150: 115-123. [22]. Samukham S, Raju G, Vyasarayani CP. Parametric instabilities of variable angle tow composite laminate under axial compression. Compos Struct 2017; 166: 229-238. [23]. Irisarri F X, Peeters DMJ, Abdalla MM. Optimisation of ply drop order in variable stiffness laminates. Compos Struct 2016; 152: 791-799. [24]. van den Brink WM, Vankan WJ, Maas R. Buckling optimized variable stiffness laminates for a composite fuselage window section. Proceedings of the 28th international congress of the aeronautical sciences, Brisbane (Australia), 2012. [25]. Lozano GG, Tiwari A, Turner C, et al. A review on design for manufacture of variable stiffness composite laminates. Proc Inst Mech Eng, B J Eng Manuf 2016; 230(6): 981-992. [26]. Huang J, Haftka RT. Optimization of fiber orientations near a hole for increased load-carrying capacity of composite laminates. Struct Multidiscip O 2005; 30(5): 335-341. [27]. Hao P, Wang B, Tian K, et al. Optimization of curvilinearly stiffened panels with single cutout concerning the collapse load. Int J Struct Stab Dy 2016; 16(07): 1550036. [28]. Hao P, Wang B, Tian K, et al. Efficient optimization of cylindrical stiffened shells with reinforced cutouts by curvilinear stiffeners. AIAA J 2016; 54(4): 1350-1363. [29]. Beldjelili Y, Tounsi A, Mahmoud SR. Hygro-thermo-mechanical bending of S-FGM plates resting on variable elastic foundations using a four-variable trigonometric plate theory. Smart Struct Syst 2016; 18(4): 755-786. [30]. Bouderba B, Houari MSA, Tounsi A, et al. Thermal stability of functionally graded sandwich plates using a simple shear deformation theory. Struct Eng Mech 2016; 58(3): 397-422. [31]. Bousahla AA, Benyoucef S, Tounsi A, et al. On thermal stability of plates with functionally graded coefficient of thermal expansion. Struct Eng Mech 2016; 60(2): 313-335. [32]. Boukhari A, Atmane HA, Tounsi A, et al. An efficient shear deformation theory for wave propagation of functionally graded material plates. Struct Eng Mech 2016; 57(5): 837-859. [33]. Houari MSA, Tounsi A, Bessaim A, et al. A new simple three-unknown sinusoidal shear deformation theory for functionally graded plates. Steel Compos Struct 2016; 22(2): 257-276.

31

[34]. Yahia SA, Atmane HA, Houari MSA, et al. Wave propagation in functionally graded plates with porosities using various higher-order shear deformation plate theories. Struct Eng Mech 2015; 53(6): 1143-1165. [35]. Attia A, Tounsi A, Bedia E A, et al. Free vibration analysis of functionally graded plates with temperature-dependent properties using various four variable refined plate theories. Smart Struct Syst 2015; 18(1): 187-212. [36]. Bourada M, Kaci A, Houari M S A, et al. A new simple shear and normal deformations theory for functionally graded beams. Steel Compos Struct 2015; 18(2): 409-423. [37]. Bouderba B, Houari MSA, Tounsi A. Thermomechanical bending response of FGM thick plates resting on WinklerPasternak elastic foundations. Steel Compos Struct 2013; 14(1): 85-104. [38]. Bellifa H, Benrahou KH, Hadji L, et al. Bending and free vibration analysis of functionally graded plates using a simple shear deformation theory and the concept the neutral surface position. J Braz Soc Mech Sci 2016; 38(1): 265-275. [39]. Besseghier A, Houari MSA, Tounsi A, et al. Free vibration analysis of embedded nanosize FG plates using a new nonlocal trigonometric shear deformation theory. Smart Struct Syst 2017; 19(6): 610-614. [40]. Liu B, Haftka RT. Composite wing structural design optimization with continuity constraints. Proceedings of 42nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA Paper 2001; 1205. [41]. Gürdal Z, Tatting BF, Wu CK. Variable stiffness composite panels: effects of stiffness variation on the in-plane and buckling response. Compos Part A-Appl S 2008; 39(5): 911-922. [42]. Muc A, Ulatowska A. Design of plates with curved fiber format. Compos Struct 2010; 92(7): 1728-1733. [43]. Blom AW, Tatting BF, Hol J MAM, et al. Fiber path definitions for elastically tailored conical shells. Compos Part B Eng 2009; 40(1): 77-84. [44]. Parnas L, Oral S, Ceyhan Ü. Optimum design of composite structures with curved fiber courses. Compos Sci Technol 2003; 63(7): 1071-82. [45]. Wu Z, Weaver PM, Raju G, et al. Buckling analysis and optimisation of variable angle tow composite plates. Thin-walled Struct 2012; 60: 163-72. [46]. Zhu Y, Liu J, Liu D, et al. Fiber path optimization based on a family of curves in composite laminate with a center hole. Compos Part B Eng 2017; 111: 91-102. [47]. Niu XJ, Yang T, Du Y, et al. Tensile properties of variable stiffness composite laminates with circular holes based on potential flow functions. Arch Appl Mech 2016; 86(9): 1551-63. [48]. Nik MA, Fayazbakhsh K, Pasini D, et al. A comparative study of metamodeling methods for the design optimization of variable stiffness composites. Compos Struct 2014; 107: 494-501. [49]. Setoodeh S, Gürdal Z, Watson LT. Design of variable-stiffness composite layers using cellular automata. Comput Method Appl M 2006; 195(9): 836-51.

32

[50]. Wu Z, Raju G, Weaver PM. Framework for the buckling optimization of variable-angle tow composite plates. AIAA J 2015; 53(12): 3788-3804. [51]. Liu D, Toropov VV. A lamination parameter-based strategy for solving an integer-continuous problem arising in composite optimization. Comput Struct 2013; 128: 170-174. [52]. Peeters DMJ, Hesse S, Abdalla MM. Stacking sequence optimisation of variable stiffness laminates with manufacturing constraints. Compos Struct 2015; 125: 596-604. [53]. Huang G, Wang H, Li G. An efficient reanalysis assisted optimization for variable-stiffness composite design by using path functions. Compos Struct 2016; 153: 409-420. [54]. Rouhi M, Ghayoor H, Hoa SV, et al. Computational efficiency and accuracy of multi-step design optimization method for variable stiffness composite structures. Thin-Walled Struct 2017; 113: 136-143. [55]. Kiyono CY, Silva ECN, Reddy JN. A novel fiber optimization method based on normal distribution function with continuously varying fiber path. Compos Struct 2017; 160: 503-515. [56]. Peeters D, Abdalla M. Optimization of ply drop locations in variable-stiffness composites. AIAA J 2016; 1760-1768. [57]. Hao P, Yuan X, Liu H, et al. Isogeometric buckling analysis of composite variable-stiffness panels. Compos Struct 2017; 165: 192-208. [58]. Montemurro M, Catapano A. On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates. Compos Struct 2017; 161: 145-159. [59]. Blom AW, Abdalla MM, Gürdal Z. Optimization of course locations in fiber-placed panels for general fiber angle distributions. Compos Sci Technol 2010; 70(4): 564-570. [60]. Havens D, Shiyekar S, Norris A, Bird RK, Kapania RK, Olliffe R. Design, optimization, and evaluation of integrally stiffened Al-2139 panel with curved stiffeners, 2011 [61]. Hao P, Wang YT, Liu XX, et al. An efficient adaptive-loop method for non-probabilistic reliability-based design optimization. Comput Method Appl M 2017; 324: 689-711.

33