- Email: [email protected]
Layer thickness optimisation in a laminated composite
Composites Part B 28B (1997) 309-317
PII: S 1 3 5 9 - 8 3 6 8 ( 9 6 ) 0 0 0 4 8 - 0
ELSEVIER
© 1997 Elsevier Science Limited Printed in Great Britain. All rights reserved 1359-8368/97/$17.00
Layer thickness optimisation in a laminated composite
P. Conti*, S. Luparello and A. Pasta Department of Mechanics and Aeronautics, University of Palermo, Viale delle Science, 90128-Palermo, Italy (Received 7 December 1994; accepted 19 February 1996) The paper describes a method to optimise the thickness balance within a composite laminate with layers oriented according to a limited set of angles. The laminate must be symmetric, balanced and loaded in-plane. The optimisation process is particularly suited to be used in conjunction with a finite element program. It provides the designer with the optimal overall engineering characteristics anda list of all the possible orientation combinations ranked with respect to their safety factor. The optimisation method is based on the first order gradient optimum search method and operates iteratively in the engineering elastic characteristics field. The cost function implemented up to now is the structure stiffness but no conceptual limitations exist for different functions. Some simple applications are listed at the end of the paper in order to verify the capabilities of the method. © 1997 Elsevier Science Limited
INTRODUCTION Most composite structures are constructed by joining together many elementary layers (pre-pregs) with different orientations with respect to an arbitrary common reference axis. Although the single layer orientation could be freely chosen between 0 ° and 360 °, in normal practice the orientations are selected among a limited number comprising at most eight angles; the aim of this work is to provide the designer with an automatic mean, based on the finite elements method (FEM), to choose the best orientation balance in a laminated composite structure with discrete values of orientations. Many optimisation methods suitable for composite materials have been developed over the years. An interesting theoretical study on optimisation of anisotropic bodies Baniciuk 1, gives a powerful conceptual framework to optimise the orientation of the fibres at every point of a structure. The method suggested cannot be easily extrapolated to provide the most convenient orientations set within a large laminated component. In earlier works Schmit and Farshi 2'3 have developed an efficient optimisation algorithm adapted to composite laminates based on the hyper sphere method (Baldur4). The method consists of the sequential search of the optimal design assuming a linearisation of both the constraint equations and the cost function (CF). At each step, the derivatives of the CF with respect to all the design variables - - in the present case, all the orientation angles - - must be * P.Conti is presently at the University of Perugia, (Istituto di Energetica, Via S. Lucia Canetola, 06125 - PERUGIA) but the work was developed during his stay at Palermo University
evaluated and the constraints linearised. As long as strain constraints are used, the method needs the evaluation of the strains in each layer of each element at each step. A large number of studies is devoted to buckling analysis and optimisation. Stroud and Anderson 5, have developed the PASCO computer code, which allows the optimisation of stiffened and unstiffened rectangular plates composed of multiple balanced angle plied symmetric laminates (MBAS laminates) with respect to buckling instability with limits on strength and stiffness. The optimisation method suggested is based on Taylor expansions of both the CF and the constraints. The extension of the method to laminates with complex shapes - - analysed by FEM - - should require the evaluation of all the derivatives at each design point and may be difficult to deal with when the CF has steep "peaks and valleys" (as often happens when the angles themselves are assumed as design variables) or when the design variables belong to a discrete set (as in the present work). Hirano 6, proposed the optimisation of fiat plates and cylindrical shells under axial compression. The most interesting aspect, from a methodological point of view, is the recourse to Powell's algorithm 7, which does not require the evaluation of the CF's derivatives. The paper yields the optimal angle, or, for many load conditions for single angleplied laminates only. Haftka and Walsh 8 and Nagendra, Haftka, G i i r d a l 9 have developed a very powerful method to optimise composite plates subjected to buckling with strain constraints. The method is based on integer programming techniques and yields the optimal stacking sequence of layers oriented according to a limited set of angles. The method is particularly suited to buckling and bending problems as it takes into account the distance of each layer from the
309
Layer thickness optimisation: P. Conti et al. mid-surface. In their present implementation, closed-form formulations of the cost function are considered. The extension to structures analysed with FEM seems to require the repeated analysis of the whole structure for all the design points selected by the integer programming algorithm. Saravanos and Chamis 1° have developed a multiobjective shape and material optimisation, which provides the best angle combination and thickness distribution. The method relies upon the choice of some control location within the structure where the thicknesses are optimised and upon the search of the angles of each single layer, which yields a design point as a compromise among all the objective functions. The use of this kind of approach is appropriate when many CF must be taken into account. In a previous work, one of the authors developed a multilevel optimisation method to choose the orientations with a finite element program, without any limit on the orientations themselves but with a limit on the thicknesses (Conti and Celia H and ContilZ). The method, however, yields orientation sets that are completely free and often are not acceptable for manufacturing reasons. The same method is used in this paper but, instead of leaving the orientation set completely free, the layer angles are chosen among a limited number and their balance is optimised. The first particular feature of the method is the design strategy, which relies upon the definition of the allowable values of the engineering elastic constants (Ell, E22, ~'~z, G t2) that can be obtained with any combination of the layer angles for MBAS laminates. An analogous design strategy has been proposed by Miki for in-plane loaded laminates Miki 14 and for flexural loaded plates Miki 15. In Miki's approach, the possible values of the elastic characteristics were deduced from the allowable values of two trigonometric functions, V~ and V2, introduced by Tsai and Pagano 16. In the authors' previous works and in the present one, the allowable characteristics of an equivalent laminate are obtained and displayed directly. This second approach is more convenient for optimisation purpose and closer to the current design practice. The second particularity is the optimisation strategy, which splits the procedure into two steps: (i) choice of the optimum overall engineering constants of the laminate and (ii) choice of an orientation set, which can provide the optimum engineering constants. A similar multilevel approach is applied for stiffened panels (Schmit and Mehrinfarl8). Most of the optimisation methods listed above require the evaluation of the first (and in some cases, the second) derivative of the CF with respect to the design variables or the linearisation of the constraints and the objective functions (Schmit and Farshi3). The advantage of this two-step approach is that the design variables during the first iterative optimisation process are always two for any number of layers; this makes the repeated FEM calculations faster and easier. Moreover, the objective functions in the [Ell, E22] domain, where the optimisation operates, are generally smoother than in the layer orientation domain and
310
the optimisation process becomes more reliable and reaches the optimum after only a few iterations. Within the two-step optimisation, many algorithms could be used; Powell's algorithm 7, Baldur's algorithm 4, etc. Up to now, however, the simplest to implement, but somewhat inefficient, algorithm was used: the steepest descent algorithm. This work does not focus on the development of an efficient optimisation software but on the usefulness of the multilevel two-step approach, The particular CF chosen is very simple to handle with a FEM program (its derivatives evaluation does not require the stiffness matrix inversion) and balances the low efficiency of the optimisation algorithm.
OPTIMISATION STRATEGY In order to use the optimisation method, some limitations must be observed: (a) Given the geometry of the structure and its weight, the optimisation process searches for the best combination of the thicknesses of the various layers oriented within a limited set of orientations. In the present implementation, eight orientations are considered: 0 °, __+ 30 °, __+ 45 °, ___ 60 °, 90 ° (but a different set could be chosen). (b) Only MBAS laminates subjected to membrane loads are considered, this means that the lay-up sequence is not taken into account in that the in-plane characteristics do not depend on the layer sequence but only on their thickness balance. (c) The optimisation defines the stiffest structure. In the present work the parameter used to maximise the stiffness of a structure is the strain energy absorbed by the structure; it is a compliance index in components with only external loads applied (and not fixed deformations). Another CF could be chosen for different loading conditions.
Conceptual framework and problem formulation The overall thickness, T, is assigned on the basis of the structure weight, then the single layer thickness, t(Oi), is considered as a continuous variable with: 0 <-- t(Oi) <-- T 8
Z t(Oi)= T i:1
Any resulting laminate can be considered as an equivalent orthotropic body with elastic characteristics depending on the i variables t(Oi). The feasible domain, which these elastic characteristics, E x, Ey, i,xy, Gxy must belong to, can be established as shown in the next paragraph. A first optimisation loop, based on the steepest-descent method, is started within the feasible domain leading to optimum equivalent orthotropic laminate. As shown later, the optimum equivalent orthotropic laminate can be obtained with various combination of at most three angles (and their corresponding thickness). This property makes the second optimisation loop very simple: all the allowable
Layer thickness optimisation: P. Conti et al.
Loads
I]
Geometry
II
MaxWeight.
]
Differentsets orientation i
Finite Element Method ~
- ~
I Singlelayer
Feasible Domain
properties
Yes
New
orientation
balance
I )
k. FIRST OPTIMISATION PROCESS
I I I I
Overall laminate character,
I
Eeng II
I I I I
Eeng 22 veng 12
[
Failure criterion i
I Rankiig Method ~ I
OPTIMISATION PROCESS
Thickness of: 0 °, +30 °, ±45°, ±60 °, ±90 ° layers Figure
1 Two step optimisation strategy
triplets are ranked with respect to the first ply failure safety factor and the designer can select the best one. In general, it may happen that the thickness attributed to an orientation is not an integer multiple of the pre-preg thickness; in that case, the system chooses the nearest layer number in an appropriate discontinuous pattern. In Figure 1 a flow diagram of the process is shown.
elastic behaviour of the laminate can be derived (see Appendix 1). We shall choose A ll and independent and Al2 and A66 as dependent. From equation (1), All and A22, can be obtained as:
A22as
All = 011 (0)X1 -~- 011 (30)X2 + 011 (45)X3 + 011(60)X4 T + 011 (90)X5
Feasible domain
A22 = Q22(0)xI -]- Q22(30)X2 + Q22(45)x3 -b Q22(60)X4 T
As only membrane loads are considered, only the in-plane stress laminate stiffness matrix IIAIImay be considered. All the layers with the same orientation can be lumped together irrespective of their real position within the laminate and an overall thickness t(O) can be attributed to the orientation 0. As symmetric and balanced laminates are considered all the elements of the matrix IIA[I are even functions of 0 and the contribution of each couple +_ 0 to any element of IIAII,can be merged in a single term. The matrix IIAIIis defined as: Aij= ZQij(Oik)tk i , j = l , 2 .... 6 k = l .... 5 k
(1)
T=tl+t2+t3+t4+t5 where 0l = 0°; 02 = 30°; 03 = 45°; 04 = 60°; 05 = 90 ° and tk is the layer thickness of the Ok orientation. Only 4 terms of IIAIIare different from zero: All, A 22, A 12 and A66. In the hypotheses of MBAS laminate, two elastic constants can be considered as independent variables while the other two variables that we need to fully characterise the
+ Q22(90)X5 where Xi = ~ and i Xi = 1 With this formulation, using a linear programming terminology, A11/T and Azz]T can be considered as components of a two-dimensional vector and the five set (Oi), Q22(0i)], i = 1,...5 can be considered as components of five other vectors. The vector [All/T, A22/ T] is then a "linear convex combination" of the other five vectors. On the basis of the properties of "linear convex combinations' ', (Muracchini 19) it can be concluded that the feasible domain for the vector [All/T, Azz/T] is a linear convex polygon having the five vectors [Q11 (Oi), Q22(0i)] vertices and that the domain boundaries are straight segments connecting the vertices. Any vertex can be related to the elastic properties of a particular orientation of the initial set. In Figure 2 the feasible domain for the couple A 11 and A22 for a particular composite material is shown. Moreover, any pair of elastic characteristics can be obtained with three orientations only. In order to understand
[QI1
as
311
Layer thickness optimisation: P. Conti et al. point A corresponds to a
A22 A
B
"
B
-
C
"
.
"
O
"
"
"
E
"
"
90" layer .+60"
where: n = total energy q ----displacement vector F = nodal forces vector IIKII = stiffness matrix
+45" -+30"
d
The stiffness matrix IIKIIis a function of many variables. Among them, we shall only consider its dependence on the two Young's moduli of the equivalent orthotropic laminate. Considering first order differentials only, the relation linking IIKIIto 6~n is (Conti and Celia 1~, Contil2):
H c
6~, = - 1/2(Iqf ll6K[llql) E
2 Feasible domain for a laminate with eight possible orientations and thickness balance choice Figure
this property, let us consider that any point P* ~ [A*II, A22]* of the feasible domain can be obtained as a linear combination of three out of the five vertex vectors defining a triangle enclosing P*. The choice of the triplet of vectors is not generally unique and many combinations of three angles can lead to the same elastic characteristics. For a given triangle, the thickness pertaining to the three orientations can be easily calculated by (see Figure 2): t(01)
HCPE BcEHT;
BHPE t(02)= ~ T ;
__PH7.
t ( 0 3 ) : EH '
where: IISKIIis a small variation of IIKII. This simple formulation of the energy variation due to changes of the stiffness matrix is powerful as it allows one to estimate A~ n without any need to invert ]IKII.The gradient components of the energy with respect to the Young's moduli, Ax~(Ex, Ey) and Ay~(Ex, Ey) can be derived just by evaluating the variations II/~K[IE~and 116KIIEy,due to small variation of Ex and Ey, respectively, and, by equations (4), yields: A x ~ ( E x, Ey) = 6~/6E x = - 1/2(IqlTIIAKIIEx[q[)/6Ex
Ay~(Ex, Ey) : ~y~/~Ey :
--
(5)
1/2(IqlTIIAKIIEyIql)/~Ey
Once the FEM analysis is performed, only the two matrices II~KIIEx and II~KIIEymust be evaluated giving very simply, first 6~,, and then Ax((Ex, Ey) and Ax((Ex, Ey) via equations
(5).
(2)
Once All and A 2z are known, all the elastic characteristics of the laminate can be determined (Jones2°): ng A~IA22 -A22 E~ll ~ A22
ALGORITHM DESCRIPTION
Optimal equivalent orthotropic laminate
E~22g = A~1A22 - a~2
All 6666g = A 6 6 1
During the optimisation process, the feasible domain boundary control is carried out in the [All, A22] domain as the linearity of the boundaries gives very simple analytical expressions.
The optimisation process, based on the steepest descent i-b-,eng, *-'22 ~,eng] method, operates in the t~ll J domain and searches the r K,opt ib-opt] optimum values of the elastic characteristics: t,-~l1 , "-'22 J" At ng ng • each step, the design point, [E~I1 , ~222 ], is converted into the [A l l,A 22] domain in order to control its position with respect to the boundaries. The process is continually repeated obtaining a new internal design point at each step (until a boundary is reached) with the following recursive formulation: ~Ex = - ~ ( A x ~ ( E x , Ey))
Cost function The energy stored in a structure with a given thickness was chosen as the cost function to minimise. In other words, the method can be used to obtain the stiffest layered structure with a given overall weight when the laminate is subject to external forces and stresses only. With the usual notation of FEM 2] the energy can be calculated as (Conti and Cella tl, Conti 12):
}= ~'qlrlF]= llqlrllKl'lql
312
(4)
(3)
6Ey= - e(Ay}(Ex, Ey)) where the step magnitude • is evaluated at each cycle on
the basis of the results of the previous cycle, in order to keep 6}I} ------0 . 1 The process can be stopped because a boundary is reached or one of the following events occurs: -the design point is in the neighbourhood of a minimum but the step magnitude is large and the process "jumps" over the minimum point leading to a new design point with ~(i) - ~(i - 1) 0. In this case the optimisation loop
Layer thickness optimisation: P. Conti et al.
,/
Figure 3
Space partition for boundary control
repeats the cycle with a halved value of e. After this procedure has been repeated for an assigned number of times, the optimisation process is stopped. in two successive cycles the greater value of the increment of the elastic moduli, e.g. Ex, is such that 6E~/Ex <-10 -5 - - a maximum number of total cycles has been reached. This limit can be set by the designer in order to prevent the program from looping. -
-
For design points belonging to the linear segments of the boundaries, only two orientations are needed to obtain the corresponding engineering characteristics and their balance can be easily evaluated just by combining with the mixture rule the two orientations corresponding to the segment vertices. Finally, if the design point is on a vertex, the corresponding orientation is sufficient to give the desired elastic characteristics. For any design point, many laminates with different orientation balances are in general acceptable. As overall elastic characteristics will be the same for all these laminates, the laminate deformation fields (and energy storage) will not be affected by the particular orientation balance considered. To the contrary, the true stress field within each lamina will be completely different with different orientation balances. Every acceptable orientation balance will be affected by its own safety factor that can be easily calculated. A post-processor routine gives the designer a ranked list of all the acceptable lay-ups, i.e. orientation balances with their safety factor. It is then possible to choose among these lay-ups the most convenient not only on the basis of their failure behaviour but also taking into account other aspects such as fatigue behaviour, interlaminar stresses, manufacturing needs, etc.
APPLICATIONS AND CONCLUSIONS
Boundary control
Applications
At each step a routine controls whether the design point (All, A22) is inside the feasible domain or outside. The external space is divided into ten regions as shown in Figure 3. If the suggested design point (A~I,A~z) is outside the feasible domain, a new feasible design point is obtained by an orthogonal projection of (A~IA~2) on the boundary. If the suggested design point belongs to regions 1,2,3,4 and 5, the corresponding vertex is chosen as new design point.
Some applications have been examined in order to verify the performance of the method. A special software has been developed in order to integrate the optimisation program with a commercial FEM code; up to now, the method was tested using the SuperSAP program ALGOR 22. The software developed allows the SuperSAP output to be processed and re-entered as a new input in an iterative way. All the tests were carried out considering a particular unidirectional carbon fibre pre-preg with the following characteristics:
Final orientation balance At the end of the first optimisation level, an optimal laminate is obtained (in which orthotropy axes are aligned with the overall reference system) and its engineering characteristics are known. The next problem is to choose the angle triplets, which can yield these engineering characteristics. All the triangles with their vertices corresponding to the vertices of the boundary of the feasible domain in the [All, A22] space (in Figure 2) are examined*. Each feasible design point (not lying on the domain boundaries) is internal to some triangle and is consistent with the orientations pertaining to the triangle vertices. The orientation balance is obtained using equations (2). If the point is external to some triangle, it means that the corresponding orientation triplet is not consistent with the design point and the triangle is skipped. * The number of triangles is n ! / [ ( n boundary vertices
-
3)! * 3!] where " n " is the number of
El = 131.4MPa E2 = 13.0 MPa Gj2 ----8.27 MPa ul2 = 0.35 thickness = 0.1 mm The simplest example was the optimisation of a laminate subject to pure tensile stress in the " y " direction. Figure 4 shows the optimisation pattern in the IAI~,A22] plane. Point P is the starting design point corresponding to a quasi isotropic laminate with Ej = E2 = 54.95 MPa, Gj2 = 21.26 MPa, and u12 ----.293. Such a laminate is obtained with a quasi isotropic lay-up: [0/90/ + 45]s. The line PQR represents the optimisation path (all the successive design points are on the line), and point R corresponds to the suggested optimal orientation. As one would expect, the optimal (stiffest) orientation is obtained with all the fibres aligned with the external load (E~ = 13 MPa, Ey = 131.4 MPa, Gt2 = 8.27 MPa, vl2 = 0.35).
313
Layer thickness optimisation: P. Conti et al. Ay
Ay R
A
Ax
Figure 4 tension
Optimisation pattern for a laminate in pure shear and in pure
Ay D
D'
Ax
Figure 5 Optimisation pattern for different starting points for a laminate in biaxial tension (o3' = 2ax)
The same laminate with the same starting lay-up was tested with a pure shear load and again the expected result - - all fibres oriented at _+ 45 ° with respect to the reference axes - - was reached (point R in Figure 4). In some applications it may happen that the objective function (energy stored in the structure) does not present a single optimal value but is characterised by a sort of linear fiat gorge. In this case the steepest descent algorithm is not able to choose a real optimum as all the points belonging to the bottom of the gorge are nearly equivalent. The only way to choose the thicknesses is to test many points belonging to the minimum line and then to choose the point corresponding to the higher safety factor. A simple example is shown: the same laminate used in the two previous examples was loaded with a combined load
314
Ax
Figure 6 Optimisation pattern with different starting points for a laminate with a hole under tensile loading
characterised by try = 2ax. The optimisation sequence in the [A ll,A22] plane is shown in Figure 5. If different starting points (A, B, C, D) are used, the optimisation converges to different optimum points (A', B', C', D'). All these points belong to the line A', B', C', D', which is the fiat bottom line of a "valley" of the objective function. The gradient along the "valley" is very low and only the comparison of the energy pertaining to the points A', B', C', D' allows one to find the true optimum, which is located at D' corresponding to a cross-ply laminate with Ex ----4.7 MPa, Ey = 9.86 MPa. The layer balance is 71.5% of the laminate with a 90 ° orientation (some possible laminates approaching this optimum balance could be [18(90°)/7(0°)]s (with 90 ° plies versus total number of plies equal to 72%) or [9(90°)/ 4(0°)]s (69% of 90 ° plies) or [7(900)/3(0°)] s (70% of 90 ° plies) depending on the actual value of oy and ax. Since the true orientation balance must be obtained with an integer numbers of plies, only with a high number of plies the optimal balance can be closely approximated. It is interesting to note that the optimal ratio of 90 ° layers versus 0 ° layers is 2.5, which is close but not equal to the ratio 0"y/O" x a s expected. Two more simple examples are shown. The first one concerns a fiat square plate with a hole (width to hole diameter ratio equal to 3) loaded in one direction. The optimisation paths in the [A ii, A22] plane (Figure 6) shows that independently of the starting point, the optimisation converges to the same optimum point, A, corresponding to a bi-directional laminate with 95% of the fibres oriented in the 90 ° direction and the remaining 5% oriented _+ 60 °. This result, with a strong predominant orientation of 90 °, can be explained: the central hole introduces a perturbation of the stress field, which is confined to the vicinity of the hole and does not significantly affect the overall energy storage. The second example is represented by a cylindrical shell loaded first with radial and axial internal pressure and then with pure torsional load at its ends. Figure 7 a shows the
Layer thickness optimisation: P. Conti et al.
(a)
Deformation Energy
(b) Deformation Energy
Figure 8
Sketch of the F1 aileron
Ey
Ay
Figure 7 (a) Optimisation patterns with different starting points (1-6) for a cylindrical shell with radial and axial pressure. (b) Energy surface (objective function) Ex
[b~eng,~22 b~eng] optimisation path in the t~ll J domain (which is oneto-one related to the [All, A22] domain) for the first load condition while Figure 7(b) shows the energy surface. The optimum is reached for the same orientation balance pertaining to a simple plate loaded in two directions with fix = 2ay; again the same "gorge" was found (and displayed in Figure 7(b)). This result was expected as the hoop stresses are twice as large as the axial stresses. Again, this result does not closely correspond to the first approximation balance: the 0 ° to 90 ° ratio is close but not equal to the ratio between the axial and hoop stresses. The method cannot be used in the present implementation for structures subject to bending or buckling. This limitation seems to be the most severe one; nevertheless, when such loads are present, sandwich structures are often used. In this case, when the core is thick enough, the sandwich can be approximated as two in-plane loaded skins - - one in tention and one under compression - - kept apart one from another by the internal unloaded core. When sandwich structures are modeled using this coarse approximation approach, the
Figure 9 Optimisation pattern for the skins of the sandwich used in the main wing construction
method can be successfully used to optimise sandwich structures too. A design application of the optimisation of a sandwich structure concerns an aileron of a F1 race car shown in Figure 8; the optimisation of the main wing will be only presented here. A detailed description of the analysis can be found in Conti and Luparello 13. The aerodynamic loads on the aileron generate both rotations and "squeezing" of the wing contour. These effects must be minimised in order not to reduce the aerodynamic effectiveness of the wing. For this reason the CF was a weighted combination of the energies pertaining only to these deformation modes and the contribution of the beam type bending of the wing, which is not relevant from an aerodynamic point of view, was neglected. The "squeezing" deformation energies were evaluated by equal virtual distributed forces applied with a different sign on the upper and lower wing surfaces while rotation energies were calculated with opposite virtual
315
Layer thickness optirnisation: P. Conti et al. characteristics orientations.
(a)
D (b)
can
be
obtained
using
only
three
The "geometric" analysis of the convex feasible domain considered as a vectorial combination is a useful approach as it very simply leads to a ranked list of possible lay-ups among which the designer can make a choice. The main advantage of this optimisation strategy is that an optimum - - the stiffest angle combination - - is suggested using only a first-order optimisation method. The procedure described can be extended to various cost functions (CFs) different from stiffness as its main features are independent of the cost function. The particular cost function was chosen because the gradient components can be evaluated without matrix inversion. The present implementation was developed on the basis of the SuperSAP program ALGOR 22, but it can be easily transferred to any finite element program. Presently the MARC implementation is under testing. The method is complementary to the method previously developed by Conti and Cella ~1 and Conti ~2. In fact the optimisation problem presents a duality in that both the thicknesses of given angles and the angles themselves can be optimised.
REFERENCES
Figure 10 Comparison of the deformations of the optimal structure (a) and the structure before the optimisation (b)
vertical concentrated forces acting along the leading and trailing edges. The aerodynamic pressures acting on the aileron were evaluated using a computer code based on a panel method and verified with true scale measure of the resultant aerodynamics loads. The optimisation process led to a laminate with 20% of the layers parallel to the wing axis and 80% with a -+ 30 ° orientations. Figure 9 shows the optimisation path in the equivalent elastic constraints domain. The correctness of the result was tested through the comparison of different configurations of the skins; Figure 10 shows a comparison between the optimal solution and a different lay-up.
1 2
3
4 5
6
7
8
9
CONCLUDING REMARKS The optimisation method developed is suitable for the optimisation of membrane loaded composite structures and carries out the optimisation with the same procedure that a designer would adopt, i.e. the selection of the thickness of single layers chosen among a limited set. The method, which can be implemented in any finite element program is characterised by two features: (a) The optimisation is carried out in two steps, first evaluating the optimal elastic characteristics and then the thicknesses themselves. (b) The optimum laminate always contains at most three different orientations since any set of elastic
316
10
11
12 13 '
14
Baniciuk, N. V. Optimization of elastic proprieties ofanisotropic bodies. Plenum Press, 1983, pp. 181-207. Schmit, L. A. and Farshi, B. Optimum laminate design for strength and stiffness. Int. J. for Numerical Methods in Eng., 1973, 7, 519536. Schmit, L. A. and Farshi, B. Optimum design of laminated fibre composite plates. Int. J. for Numerical Methods in Eng., 1977, 11, 623-640. Baldur, R. Structural optimization by inscribed hyperspheres. Proc. ASCE. Engineering Mechanics Division, 1972, 98(3), 503-518. Stroud, J. and Anderson, M. S. PASCO: Structural Panel Analysis and Sizing Code, Capability and Analytical Foundations. NASA TM-80181, 1980. Hirano, Y. Optimization of laminated composite plates and shells. Proc. IUTAM Syrup., eds Hashin and C. T. Herakovich, 1982, pp. 355-365. Powell, M. J. D. An efficient method for finding the minimum of a function of several variables without calculating derivatives. Computer Journal, 1964, 7, 155-162. Haftka, R. T. and Walsh, J. L. Stacking sequence optimization for buckling of laminated plates by integer programming. AIAA J., 1992, 30(3), 814-819. Nagendra, S., Haftka, R. T. and Gtirdal, Z. Stacking sequence optimization of simply supported laminates with stability and strain constraints. In Proc. 33rd SDM Conf., Dallas, Texas, AIAA-922310-CP., 1992, pp. 2526-2535. Saravanos, D. A. and Chamis, C. C. Multiobjective optimum design of structures with fibre-composite stiffened panel elements. A1AA Journal, 1992, 30(3), 805-813. Conti, P. and Cella, A. An optimal design of multilayered laminates based on finite element stress analysis. Composite Material Technology, eds D. Hui, T. J. Kozik and O. O. Ochoa, ASME, PD-45, 1992, Book No. G00653. Conti, P. Computer aided optimization of composite structures. Numerical Methods for Structural Integrity Evaluation, eds E. Vitale and S. Kossev, ETS Editrice, Pisa, Italy, 1991, pp. 329-364. Conti, P. and Luparello, S. Design optimisation of the fibres pattern for a F1 car stabilizer in composite material. Proc. of "Modena Motori, Tradizioni, Tecnologie e Futuro dell'Auto Sportiva". II Conf.-Modena, Italy, 1995. Miki, M, Material Design of Composite Laminates with Required In-Plane Elastic Properties. Progress in Science and Engineering of
Layer thickness optimisation: P. Conti et al.
15
16
17 18 19 20 21 22
Composites. (eds T. Hayashi, K. Kawata and S. Umekawa), ICCMIV, Tokio, 1982, pp. 1725-1731. Miki, M. Design of laminated fibrous composite plates with required flexural stiffness. Recent advances in Composites in the U.S. and Japan. ASTM STP 864, eds J. R. Vinson and M. Taya, 1983, pp. 387-400. Tsai, S. W. and Pagano, N. J. Invariant proprieties of composite materials. Composite Materials Workshop, Technomic publishing Co., Wesport, Conn., eds S. W. Tsai, L. J. C. Halpin and N. J. Pagano, 1963, pp. 233-253. Tsai, S. W. Composite design. Think Composites. Section 13, 1986. Schmit, L. A. and Mehrinfar, Multilevel optimum design of structures with fibre-composite stiffened-panel components. AIAA J., 1982, 20, 138-147. Muracchini, L. Programmazione matematica. UTET, Torino (in Italian), 1975. Jones, R. M. Mechanics of Composite Materials. McGraw-Hill Kogakusha Ltd., Tokio, 1975. Zienkewicz, O. C. The Finite Element Method. Second Ed. McGraw-Hill, 1971. ALGOR, Finite element Analysis System, Composite Decoder and Element Release Notes, Algor Interactive Systems, INC., 260 Alpha Drive, Pittsburgh, PA, 15238, April 1990.
Al2 = u4VoA + u3V3A A66 = U5VoA --~ 1,13V3A
with
VoA=T ViA = ~
[COS(/+ 1)akl(hk - hk-1)
K=I
where: n = n u m b e r o f layers hk = d i s t a n c e f r o m the m i d surface o f the k th l a y e r interface ui = i th r o t a t i o n i n v a r i a n t
Combination of equations (A1) leads to: A l l + A 2 2 nt- 2A12 = T(u5 - u4)
(A2)
A66 - AI2 = T(u 5 - u4)
APPENDIX Using the notation of Tsai ~7, in a balanced and symmetric laminate, the matrix [IAfrcan be expressed in terms of invariants by: All = u l VOA -I- u2 VIA "~ U3 V3A
Equations (A2) show that two out of the four components of IIAII can be derived from the other two and therefore, only two components are really independent and can fully define the in-plane behaviour of a symmetric and balanced laminate in its orthotropy plane.
(A1)
A22 = u I VOA + u2 VIA -+ u3 V3A
317
PII: S 1 3 5 9 - 8 3 6 8 ( 9 6 ) 0 0 0 4 8 - 0
ELSEVIER
© 1997 Elsevier Science Limited Printed in Great Britain. All rights reserved 1359-8368/97/$17.00
Layer thickness optimisation in a laminated composite
P. Conti*, S. Luparello and A. Pasta Department of Mechanics and Aeronautics, University of Palermo, Viale delle Science, 90128-Palermo, Italy (Received 7 December 1994; accepted 19 February 1996) The paper describes a method to optimise the thickness balance within a composite laminate with layers oriented according to a limited set of angles. The laminate must be symmetric, balanced and loaded in-plane. The optimisation process is particularly suited to be used in conjunction with a finite element program. It provides the designer with the optimal overall engineering characteristics anda list of all the possible orientation combinations ranked with respect to their safety factor. The optimisation method is based on the first order gradient optimum search method and operates iteratively in the engineering elastic characteristics field. The cost function implemented up to now is the structure stiffness but no conceptual limitations exist for different functions. Some simple applications are listed at the end of the paper in order to verify the capabilities of the method. © 1997 Elsevier Science Limited
INTRODUCTION Most composite structures are constructed by joining together many elementary layers (pre-pregs) with different orientations with respect to an arbitrary common reference axis. Although the single layer orientation could be freely chosen between 0 ° and 360 °, in normal practice the orientations are selected among a limited number comprising at most eight angles; the aim of this work is to provide the designer with an automatic mean, based on the finite elements method (FEM), to choose the best orientation balance in a laminated composite structure with discrete values of orientations. Many optimisation methods suitable for composite materials have been developed over the years. An interesting theoretical study on optimisation of anisotropic bodies Baniciuk 1, gives a powerful conceptual framework to optimise the orientation of the fibres at every point of a structure. The method suggested cannot be easily extrapolated to provide the most convenient orientations set within a large laminated component. In earlier works Schmit and Farshi 2'3 have developed an efficient optimisation algorithm adapted to composite laminates based on the hyper sphere method (Baldur4). The method consists of the sequential search of the optimal design assuming a linearisation of both the constraint equations and the cost function (CF). At each step, the derivatives of the CF with respect to all the design variables - - in the present case, all the orientation angles - - must be * P.Conti is presently at the University of Perugia, (Istituto di Energetica, Via S. Lucia Canetola, 06125 - PERUGIA) but the work was developed during his stay at Palermo University
evaluated and the constraints linearised. As long as strain constraints are used, the method needs the evaluation of the strains in each layer of each element at each step. A large number of studies is devoted to buckling analysis and optimisation. Stroud and Anderson 5, have developed the PASCO computer code, which allows the optimisation of stiffened and unstiffened rectangular plates composed of multiple balanced angle plied symmetric laminates (MBAS laminates) with respect to buckling instability with limits on strength and stiffness. The optimisation method suggested is based on Taylor expansions of both the CF and the constraints. The extension of the method to laminates with complex shapes - - analysed by FEM - - should require the evaluation of all the derivatives at each design point and may be difficult to deal with when the CF has steep "peaks and valleys" (as often happens when the angles themselves are assumed as design variables) or when the design variables belong to a discrete set (as in the present work). Hirano 6, proposed the optimisation of fiat plates and cylindrical shells under axial compression. The most interesting aspect, from a methodological point of view, is the recourse to Powell's algorithm 7, which does not require the evaluation of the CF's derivatives. The paper yields the optimal angle, or, for many load conditions for single angleplied laminates only. Haftka and Walsh 8 and Nagendra, Haftka, G i i r d a l 9 have developed a very powerful method to optimise composite plates subjected to buckling with strain constraints. The method is based on integer programming techniques and yields the optimal stacking sequence of layers oriented according to a limited set of angles. The method is particularly suited to buckling and bending problems as it takes into account the distance of each layer from the
309
Layer thickness optimisation: P. Conti et al. mid-surface. In their present implementation, closed-form formulations of the cost function are considered. The extension to structures analysed with FEM seems to require the repeated analysis of the whole structure for all the design points selected by the integer programming algorithm. Saravanos and Chamis 1° have developed a multiobjective shape and material optimisation, which provides the best angle combination and thickness distribution. The method relies upon the choice of some control location within the structure where the thicknesses are optimised and upon the search of the angles of each single layer, which yields a design point as a compromise among all the objective functions. The use of this kind of approach is appropriate when many CF must be taken into account. In a previous work, one of the authors developed a multilevel optimisation method to choose the orientations with a finite element program, without any limit on the orientations themselves but with a limit on the thicknesses (Conti and Celia H and ContilZ). The method, however, yields orientation sets that are completely free and often are not acceptable for manufacturing reasons. The same method is used in this paper but, instead of leaving the orientation set completely free, the layer angles are chosen among a limited number and their balance is optimised. The first particular feature of the method is the design strategy, which relies upon the definition of the allowable values of the engineering elastic constants (Ell, E22, ~'~z, G t2) that can be obtained with any combination of the layer angles for MBAS laminates. An analogous design strategy has been proposed by Miki for in-plane loaded laminates Miki 14 and for flexural loaded plates Miki 15. In Miki's approach, the possible values of the elastic characteristics were deduced from the allowable values of two trigonometric functions, V~ and V2, introduced by Tsai and Pagano 16. In the authors' previous works and in the present one, the allowable characteristics of an equivalent laminate are obtained and displayed directly. This second approach is more convenient for optimisation purpose and closer to the current design practice. The second particularity is the optimisation strategy, which splits the procedure into two steps: (i) choice of the optimum overall engineering constants of the laminate and (ii) choice of an orientation set, which can provide the optimum engineering constants. A similar multilevel approach is applied for stiffened panels (Schmit and Mehrinfarl8). Most of the optimisation methods listed above require the evaluation of the first (and in some cases, the second) derivative of the CF with respect to the design variables or the linearisation of the constraints and the objective functions (Schmit and Farshi3). The advantage of this two-step approach is that the design variables during the first iterative optimisation process are always two for any number of layers; this makes the repeated FEM calculations faster and easier. Moreover, the objective functions in the [Ell, E22] domain, where the optimisation operates, are generally smoother than in the layer orientation domain and
310
the optimisation process becomes more reliable and reaches the optimum after only a few iterations. Within the two-step optimisation, many algorithms could be used; Powell's algorithm 7, Baldur's algorithm 4, etc. Up to now, however, the simplest to implement, but somewhat inefficient, algorithm was used: the steepest descent algorithm. This work does not focus on the development of an efficient optimisation software but on the usefulness of the multilevel two-step approach, The particular CF chosen is very simple to handle with a FEM program (its derivatives evaluation does not require the stiffness matrix inversion) and balances the low efficiency of the optimisation algorithm.
OPTIMISATION STRATEGY In order to use the optimisation method, some limitations must be observed: (a) Given the geometry of the structure and its weight, the optimisation process searches for the best combination of the thicknesses of the various layers oriented within a limited set of orientations. In the present implementation, eight orientations are considered: 0 °, __+ 30 °, __+ 45 °, ___ 60 °, 90 ° (but a different set could be chosen). (b) Only MBAS laminates subjected to membrane loads are considered, this means that the lay-up sequence is not taken into account in that the in-plane characteristics do not depend on the layer sequence but only on their thickness balance. (c) The optimisation defines the stiffest structure. In the present work the parameter used to maximise the stiffness of a structure is the strain energy absorbed by the structure; it is a compliance index in components with only external loads applied (and not fixed deformations). Another CF could be chosen for different loading conditions.
Conceptual framework and problem formulation The overall thickness, T, is assigned on the basis of the structure weight, then the single layer thickness, t(Oi), is considered as a continuous variable with: 0 <-- t(Oi) <-- T 8
Z t(Oi)= T i:1
Any resulting laminate can be considered as an equivalent orthotropic body with elastic characteristics depending on the i variables t(Oi). The feasible domain, which these elastic characteristics, E x, Ey, i,xy, Gxy must belong to, can be established as shown in the next paragraph. A first optimisation loop, based on the steepest-descent method, is started within the feasible domain leading to optimum equivalent orthotropic laminate. As shown later, the optimum equivalent orthotropic laminate can be obtained with various combination of at most three angles (and their corresponding thickness). This property makes the second optimisation loop very simple: all the allowable
Layer thickness optimisation: P. Conti et al.
Loads
I]
Geometry
II
MaxWeight.
]
Differentsets orientation i
Finite Element Method ~
- ~
I Singlelayer
Feasible Domain
properties
Yes
New
orientation
balance
I )
k. FIRST OPTIMISATION PROCESS
I I I I
Overall laminate character,
I
Eeng II
I I I I
Eeng 22 veng 12
[
Failure criterion i
I Rankiig Method ~ I
OPTIMISATION PROCESS
Thickness of: 0 °, +30 °, ±45°, ±60 °, ±90 ° layers Figure
1 Two step optimisation strategy
triplets are ranked with respect to the first ply failure safety factor and the designer can select the best one. In general, it may happen that the thickness attributed to an orientation is not an integer multiple of the pre-preg thickness; in that case, the system chooses the nearest layer number in an appropriate discontinuous pattern. In Figure 1 a flow diagram of the process is shown.
elastic behaviour of the laminate can be derived (see Appendix 1). We shall choose A ll and independent and Al2 and A66 as dependent. From equation (1), All and A22, can be obtained as:
A22as
All = 011 (0)X1 -~- 011 (30)X2 + 011 (45)X3 + 011(60)X4 T + 011 (90)X5
Feasible domain
A22 = Q22(0)xI -]- Q22(30)X2 + Q22(45)x3 -b Q22(60)X4 T
As only membrane loads are considered, only the in-plane stress laminate stiffness matrix IIAIImay be considered. All the layers with the same orientation can be lumped together irrespective of their real position within the laminate and an overall thickness t(O) can be attributed to the orientation 0. As symmetric and balanced laminates are considered all the elements of the matrix IIA[I are even functions of 0 and the contribution of each couple +_ 0 to any element of IIAII,can be merged in a single term. The matrix IIAIIis defined as: Aij= ZQij(Oik)tk i , j = l , 2 .... 6 k = l .... 5 k
(1)
T=tl+t2+t3+t4+t5 where 0l = 0°; 02 = 30°; 03 = 45°; 04 = 60°; 05 = 90 ° and tk is the layer thickness of the Ok orientation. Only 4 terms of IIAIIare different from zero: All, A 22, A 12 and A66. In the hypotheses of MBAS laminate, two elastic constants can be considered as independent variables while the other two variables that we need to fully characterise the
+ Q22(90)X5 where Xi = ~ and i Xi = 1 With this formulation, using a linear programming terminology, A11/T and Azz]T can be considered as components of a two-dimensional vector and the five set (Oi), Q22(0i)], i = 1,...5 can be considered as components of five other vectors. The vector [All/T, A22/ T] is then a "linear convex combination" of the other five vectors. On the basis of the properties of "linear convex combinations' ', (Muracchini 19) it can be concluded that the feasible domain for the vector [All/T, Azz/T] is a linear convex polygon having the five vectors [Q11 (Oi), Q22(0i)] vertices and that the domain boundaries are straight segments connecting the vertices. Any vertex can be related to the elastic properties of a particular orientation of the initial set. In Figure 2 the feasible domain for the couple A 11 and A22 for a particular composite material is shown. Moreover, any pair of elastic characteristics can be obtained with three orientations only. In order to understand
[QI1
as
311
Layer thickness optimisation: P. Conti et al. point A corresponds to a
A22 A
B
"
B
-
C
"
.
"
O
"
"
"
E
"
"
90" layer .+60"
where: n = total energy q ----displacement vector F = nodal forces vector IIKII = stiffness matrix
+45" -+30"
d
The stiffness matrix IIKIIis a function of many variables. Among them, we shall only consider its dependence on the two Young's moduli of the equivalent orthotropic laminate. Considering first order differentials only, the relation linking IIKIIto 6~n is (Conti and Celia 1~, Contil2):
H c
6~, = - 1/2(Iqf ll6K[llql) E
2 Feasible domain for a laminate with eight possible orientations and thickness balance choice Figure
this property, let us consider that any point P* ~ [A*II, A22]* of the feasible domain can be obtained as a linear combination of three out of the five vertex vectors defining a triangle enclosing P*. The choice of the triplet of vectors is not generally unique and many combinations of three angles can lead to the same elastic characteristics. For a given triangle, the thickness pertaining to the three orientations can be easily calculated by (see Figure 2): t(01)
HCPE BcEHT;
BHPE t(02)= ~ T ;
__PH7.
t ( 0 3 ) : EH '
where: IISKIIis a small variation of IIKII. This simple formulation of the energy variation due to changes of the stiffness matrix is powerful as it allows one to estimate A~ n without any need to invert ]IKII.The gradient components of the energy with respect to the Young's moduli, Ax~(Ex, Ey) and Ay~(Ex, Ey) can be derived just by evaluating the variations II/~K[IE~and 116KIIEy,due to small variation of Ex and Ey, respectively, and, by equations (4), yields: A x ~ ( E x, Ey) = 6~/6E x = - 1/2(IqlTIIAKIIEx[q[)/6Ex
Ay~(Ex, Ey) : ~y~/~Ey :
--
(5)
1/2(IqlTIIAKIIEyIql)/~Ey
Once the FEM analysis is performed, only the two matrices II~KIIEx and II~KIIEymust be evaluated giving very simply, first 6~,, and then Ax((Ex, Ey) and Ax((Ex, Ey) via equations
(5).
(2)
Once All and A 2z are known, all the elastic characteristics of the laminate can be determined (Jones2°): ng A~IA22 -A22 E~ll ~ A22
ALGORITHM DESCRIPTION
Optimal equivalent orthotropic laminate
E~22g = A~1A22 - a~2
All 6666g = A 6 6 1
During the optimisation process, the feasible domain boundary control is carried out in the [All, A22] domain as the linearity of the boundaries gives very simple analytical expressions.
The optimisation process, based on the steepest descent i-b-,eng, *-'22 ~,eng] method, operates in the t~ll J domain and searches the r K,opt ib-opt] optimum values of the elastic characteristics: t,-~l1 , "-'22 J" At ng ng • each step, the design point, [E~I1 , ~222 ], is converted into the [A l l,A 22] domain in order to control its position with respect to the boundaries. The process is continually repeated obtaining a new internal design point at each step (until a boundary is reached) with the following recursive formulation: ~Ex = - ~ ( A x ~ ( E x , Ey))
Cost function The energy stored in a structure with a given thickness was chosen as the cost function to minimise. In other words, the method can be used to obtain the stiffest layered structure with a given overall weight when the laminate is subject to external forces and stresses only. With the usual notation of FEM 2] the energy can be calculated as (Conti and Cella tl, Conti 12):
}= ~'qlrlF]= llqlrllKl'lql
312
(4)
(3)
6Ey= - e(Ay}(Ex, Ey)) where the step magnitude • is evaluated at each cycle on
the basis of the results of the previous cycle, in order to keep 6}I} ------0 . 1 The process can be stopped because a boundary is reached or one of the following events occurs: -the design point is in the neighbourhood of a minimum but the step magnitude is large and the process "jumps" over the minimum point leading to a new design point with ~(i) - ~(i - 1) 0. In this case the optimisation loop
Layer thickness optimisation: P. Conti et al.
,/
Figure 3
Space partition for boundary control
repeats the cycle with a halved value of e. After this procedure has been repeated for an assigned number of times, the optimisation process is stopped. in two successive cycles the greater value of the increment of the elastic moduli, e.g. Ex, is such that 6E~/Ex <-10 -5 - - a maximum number of total cycles has been reached. This limit can be set by the designer in order to prevent the program from looping. -
-
For design points belonging to the linear segments of the boundaries, only two orientations are needed to obtain the corresponding engineering characteristics and their balance can be easily evaluated just by combining with the mixture rule the two orientations corresponding to the segment vertices. Finally, if the design point is on a vertex, the corresponding orientation is sufficient to give the desired elastic characteristics. For any design point, many laminates with different orientation balances are in general acceptable. As overall elastic characteristics will be the same for all these laminates, the laminate deformation fields (and energy storage) will not be affected by the particular orientation balance considered. To the contrary, the true stress field within each lamina will be completely different with different orientation balances. Every acceptable orientation balance will be affected by its own safety factor that can be easily calculated. A post-processor routine gives the designer a ranked list of all the acceptable lay-ups, i.e. orientation balances with their safety factor. It is then possible to choose among these lay-ups the most convenient not only on the basis of their failure behaviour but also taking into account other aspects such as fatigue behaviour, interlaminar stresses, manufacturing needs, etc.
APPLICATIONS AND CONCLUSIONS
Boundary control
Applications
At each step a routine controls whether the design point (All, A22) is inside the feasible domain or outside. The external space is divided into ten regions as shown in Figure 3. If the suggested design point (A~I,A~z) is outside the feasible domain, a new feasible design point is obtained by an orthogonal projection of (A~IA~2) on the boundary. If the suggested design point belongs to regions 1,2,3,4 and 5, the corresponding vertex is chosen as new design point.
Some applications have been examined in order to verify the performance of the method. A special software has been developed in order to integrate the optimisation program with a commercial FEM code; up to now, the method was tested using the SuperSAP program ALGOR 22. The software developed allows the SuperSAP output to be processed and re-entered as a new input in an iterative way. All the tests were carried out considering a particular unidirectional carbon fibre pre-preg with the following characteristics:
Final orientation balance At the end of the first optimisation level, an optimal laminate is obtained (in which orthotropy axes are aligned with the overall reference system) and its engineering characteristics are known. The next problem is to choose the angle triplets, which can yield these engineering characteristics. All the triangles with their vertices corresponding to the vertices of the boundary of the feasible domain in the [All, A22] space (in Figure 2) are examined*. Each feasible design point (not lying on the domain boundaries) is internal to some triangle and is consistent with the orientations pertaining to the triangle vertices. The orientation balance is obtained using equations (2). If the point is external to some triangle, it means that the corresponding orientation triplet is not consistent with the design point and the triangle is skipped. * The number of triangles is n ! / [ ( n boundary vertices
-
3)! * 3!] where " n " is the number of
El = 131.4MPa E2 = 13.0 MPa Gj2 ----8.27 MPa ul2 = 0.35 thickness = 0.1 mm The simplest example was the optimisation of a laminate subject to pure tensile stress in the " y " direction. Figure 4 shows the optimisation pattern in the IAI~,A22] plane. Point P is the starting design point corresponding to a quasi isotropic laminate with Ej = E2 = 54.95 MPa, Gj2 = 21.26 MPa, and u12 ----.293. Such a laminate is obtained with a quasi isotropic lay-up: [0/90/ + 45]s. The line PQR represents the optimisation path (all the successive design points are on the line), and point R corresponds to the suggested optimal orientation. As one would expect, the optimal (stiffest) orientation is obtained with all the fibres aligned with the external load (E~ = 13 MPa, Ey = 131.4 MPa, Gt2 = 8.27 MPa, vl2 = 0.35).
313
Layer thickness optimisation: P. Conti et al. Ay
Ay R
A
Ax
Figure 4 tension
Optimisation pattern for a laminate in pure shear and in pure
Ay D
D'
Ax
Figure 5 Optimisation pattern for different starting points for a laminate in biaxial tension (o3' = 2ax)
The same laminate with the same starting lay-up was tested with a pure shear load and again the expected result - - all fibres oriented at _+ 45 ° with respect to the reference axes - - was reached (point R in Figure 4). In some applications it may happen that the objective function (energy stored in the structure) does not present a single optimal value but is characterised by a sort of linear fiat gorge. In this case the steepest descent algorithm is not able to choose a real optimum as all the points belonging to the bottom of the gorge are nearly equivalent. The only way to choose the thicknesses is to test many points belonging to the minimum line and then to choose the point corresponding to the higher safety factor. A simple example is shown: the same laminate used in the two previous examples was loaded with a combined load
314
Ax
Figure 6 Optimisation pattern with different starting points for a laminate with a hole under tensile loading
characterised by try = 2ax. The optimisation sequence in the [A ll,A22] plane is shown in Figure 5. If different starting points (A, B, C, D) are used, the optimisation converges to different optimum points (A', B', C', D'). All these points belong to the line A', B', C', D', which is the fiat bottom line of a "valley" of the objective function. The gradient along the "valley" is very low and only the comparison of the energy pertaining to the points A', B', C', D' allows one to find the true optimum, which is located at D' corresponding to a cross-ply laminate with Ex ----4.7 MPa, Ey = 9.86 MPa. The layer balance is 71.5% of the laminate with a 90 ° orientation (some possible laminates approaching this optimum balance could be [18(90°)/7(0°)]s (with 90 ° plies versus total number of plies equal to 72%) or [9(90°)/ 4(0°)]s (69% of 90 ° plies) or [7(900)/3(0°)] s (70% of 90 ° plies) depending on the actual value of oy and ax. Since the true orientation balance must be obtained with an integer numbers of plies, only with a high number of plies the optimal balance can be closely approximated. It is interesting to note that the optimal ratio of 90 ° layers versus 0 ° layers is 2.5, which is close but not equal to the ratio 0"y/O" x a s expected. Two more simple examples are shown. The first one concerns a fiat square plate with a hole (width to hole diameter ratio equal to 3) loaded in one direction. The optimisation paths in the [A ii, A22] plane (Figure 6) shows that independently of the starting point, the optimisation converges to the same optimum point, A, corresponding to a bi-directional laminate with 95% of the fibres oriented in the 90 ° direction and the remaining 5% oriented _+ 60 °. This result, with a strong predominant orientation of 90 °, can be explained: the central hole introduces a perturbation of the stress field, which is confined to the vicinity of the hole and does not significantly affect the overall energy storage. The second example is represented by a cylindrical shell loaded first with radial and axial internal pressure and then with pure torsional load at its ends. Figure 7 a shows the
Layer thickness optimisation: P. Conti et al.
(a)
Deformation Energy
(b) Deformation Energy
Figure 8
Sketch of the F1 aileron
Ey
Ay
Figure 7 (a) Optimisation patterns with different starting points (1-6) for a cylindrical shell with radial and axial pressure. (b) Energy surface (objective function) Ex
[b~eng,~22 b~eng] optimisation path in the t~ll J domain (which is oneto-one related to the [All, A22] domain) for the first load condition while Figure 7(b) shows the energy surface. The optimum is reached for the same orientation balance pertaining to a simple plate loaded in two directions with fix = 2ay; again the same "gorge" was found (and displayed in Figure 7(b)). This result was expected as the hoop stresses are twice as large as the axial stresses. Again, this result does not closely correspond to the first approximation balance: the 0 ° to 90 ° ratio is close but not equal to the ratio between the axial and hoop stresses. The method cannot be used in the present implementation for structures subject to bending or buckling. This limitation seems to be the most severe one; nevertheless, when such loads are present, sandwich structures are often used. In this case, when the core is thick enough, the sandwich can be approximated as two in-plane loaded skins - - one in tention and one under compression - - kept apart one from another by the internal unloaded core. When sandwich structures are modeled using this coarse approximation approach, the
Figure 9 Optimisation pattern for the skins of the sandwich used in the main wing construction
method can be successfully used to optimise sandwich structures too. A design application of the optimisation of a sandwich structure concerns an aileron of a F1 race car shown in Figure 8; the optimisation of the main wing will be only presented here. A detailed description of the analysis can be found in Conti and Luparello 13. The aerodynamic loads on the aileron generate both rotations and "squeezing" of the wing contour. These effects must be minimised in order not to reduce the aerodynamic effectiveness of the wing. For this reason the CF was a weighted combination of the energies pertaining only to these deformation modes and the contribution of the beam type bending of the wing, which is not relevant from an aerodynamic point of view, was neglected. The "squeezing" deformation energies were evaluated by equal virtual distributed forces applied with a different sign on the upper and lower wing surfaces while rotation energies were calculated with opposite virtual
315
Layer thickness optirnisation: P. Conti et al. characteristics orientations.
(a)
D (b)
can
be
obtained
using
only
three
The "geometric" analysis of the convex feasible domain considered as a vectorial combination is a useful approach as it very simply leads to a ranked list of possible lay-ups among which the designer can make a choice. The main advantage of this optimisation strategy is that an optimum - - the stiffest angle combination - - is suggested using only a first-order optimisation method. The procedure described can be extended to various cost functions (CFs) different from stiffness as its main features are independent of the cost function. The particular cost function was chosen because the gradient components can be evaluated without matrix inversion. The present implementation was developed on the basis of the SuperSAP program ALGOR 22, but it can be easily transferred to any finite element program. Presently the MARC implementation is under testing. The method is complementary to the method previously developed by Conti and Cella ~1 and Conti ~2. In fact the optimisation problem presents a duality in that both the thicknesses of given angles and the angles themselves can be optimised.
REFERENCES
Figure 10 Comparison of the deformations of the optimal structure (a) and the structure before the optimisation (b)
vertical concentrated forces acting along the leading and trailing edges. The aerodynamic pressures acting on the aileron were evaluated using a computer code based on a panel method and verified with true scale measure of the resultant aerodynamics loads. The optimisation process led to a laminate with 20% of the layers parallel to the wing axis and 80% with a -+ 30 ° orientations. Figure 9 shows the optimisation path in the equivalent elastic constraints domain. The correctness of the result was tested through the comparison of different configurations of the skins; Figure 10 shows a comparison between the optimal solution and a different lay-up.
1 2
3
4 5
6
7
8
9
CONCLUDING REMARKS The optimisation method developed is suitable for the optimisation of membrane loaded composite structures and carries out the optimisation with the same procedure that a designer would adopt, i.e. the selection of the thickness of single layers chosen among a limited set. The method, which can be implemented in any finite element program is characterised by two features: (a) The optimisation is carried out in two steps, first evaluating the optimal elastic characteristics and then the thicknesses themselves. (b) The optimum laminate always contains at most three different orientations since any set of elastic
316
10
11
12 13 '
14
Baniciuk, N. V. Optimization of elastic proprieties ofanisotropic bodies. Plenum Press, 1983, pp. 181-207. Schmit, L. A. and Farshi, B. Optimum laminate design for strength and stiffness. Int. J. for Numerical Methods in Eng., 1973, 7, 519536. Schmit, L. A. and Farshi, B. Optimum design of laminated fibre composite plates. Int. J. for Numerical Methods in Eng., 1977, 11, 623-640. Baldur, R. Structural optimization by inscribed hyperspheres. Proc. ASCE. Engineering Mechanics Division, 1972, 98(3), 503-518. Stroud, J. and Anderson, M. S. PASCO: Structural Panel Analysis and Sizing Code, Capability and Analytical Foundations. NASA TM-80181, 1980. Hirano, Y. Optimization of laminated composite plates and shells. Proc. IUTAM Syrup., eds Hashin and C. T. Herakovich, 1982, pp. 355-365. Powell, M. J. D. An efficient method for finding the minimum of a function of several variables without calculating derivatives. Computer Journal, 1964, 7, 155-162. Haftka, R. T. and Walsh, J. L. Stacking sequence optimization for buckling of laminated plates by integer programming. AIAA J., 1992, 30(3), 814-819. Nagendra, S., Haftka, R. T. and Gtirdal, Z. Stacking sequence optimization of simply supported laminates with stability and strain constraints. In Proc. 33rd SDM Conf., Dallas, Texas, AIAA-922310-CP., 1992, pp. 2526-2535. Saravanos, D. A. and Chamis, C. C. Multiobjective optimum design of structures with fibre-composite stiffened panel elements. A1AA Journal, 1992, 30(3), 805-813. Conti, P. and Cella, A. An optimal design of multilayered laminates based on finite element stress analysis. Composite Material Technology, eds D. Hui, T. J. Kozik and O. O. Ochoa, ASME, PD-45, 1992, Book No. G00653. Conti, P. Computer aided optimization of composite structures. Numerical Methods for Structural Integrity Evaluation, eds E. Vitale and S. Kossev, ETS Editrice, Pisa, Italy, 1991, pp. 329-364. Conti, P. and Luparello, S. Design optimisation of the fibres pattern for a F1 car stabilizer in composite material. Proc. of "Modena Motori, Tradizioni, Tecnologie e Futuro dell'Auto Sportiva". II Conf.-Modena, Italy, 1995. Miki, M, Material Design of Composite Laminates with Required In-Plane Elastic Properties. Progress in Science and Engineering of
Layer thickness optimisation: P. Conti et al.
15
16
17 18 19 20 21 22
Composites. (eds T. Hayashi, K. Kawata and S. Umekawa), ICCMIV, Tokio, 1982, pp. 1725-1731. Miki, M. Design of laminated fibrous composite plates with required flexural stiffness. Recent advances in Composites in the U.S. and Japan. ASTM STP 864, eds J. R. Vinson and M. Taya, 1983, pp. 387-400. Tsai, S. W. and Pagano, N. J. Invariant proprieties of composite materials. Composite Materials Workshop, Technomic publishing Co., Wesport, Conn., eds S. W. Tsai, L. J. C. Halpin and N. J. Pagano, 1963, pp. 233-253. Tsai, S. W. Composite design. Think Composites. Section 13, 1986. Schmit, L. A. and Mehrinfar, Multilevel optimum design of structures with fibre-composite stiffened-panel components. AIAA J., 1982, 20, 138-147. Muracchini, L. Programmazione matematica. UTET, Torino (in Italian), 1975. Jones, R. M. Mechanics of Composite Materials. McGraw-Hill Kogakusha Ltd., Tokio, 1975. Zienkewicz, O. C. The Finite Element Method. Second Ed. McGraw-Hill, 1971. ALGOR, Finite element Analysis System, Composite Decoder and Element Release Notes, Algor Interactive Systems, INC., 260 Alpha Drive, Pittsburgh, PA, 15238, April 1990.
Al2 = u4VoA + u3V3A A66 = U5VoA --~ 1,13V3A
with
VoA=T ViA = ~
[COS(/+ 1)akl(hk - hk-1)
K=I
where: n = n u m b e r o f layers hk = d i s t a n c e f r o m the m i d surface o f the k th l a y e r interface ui = i th r o t a t i o n i n v a r i a n t
Combination of equations (A1) leads to: A l l + A 2 2 nt- 2A12 = T(u5 - u4)
(A2)
A66 - AI2 = T(u 5 - u4)
APPENDIX Using the notation of Tsai ~7, in a balanced and symmetric laminate, the matrix [IAfrcan be expressed in terms of invariants by: All = u l VOA -I- u2 VIA "~ U3 V3A
Equations (A2) show that two out of the four components of IIAII can be derived from the other two and therefore, only two components are really independent and can fully define the in-plane behaviour of a symmetric and balanced laminate in its orthotropy plane.
(A1)
A22 = u I VOA + u2 VIA -+ u3 V3A
317