OptimumDesign of Fibre Orientation in Composite Laminate Plates for Out-Plane Stresses

Previous studies have shown that composite fibre orientations can be optimised for specific load cases such as longitudinal or inplane loading. However, the methodologies utilised in these studies cannot be used for general analysis of such problems. In this research, an extra term is added to the optimisation penalty function in order to consider the transverse shear effect. This modified penalty function leads to a new methodology whereby the thickness of laminated composite plate is minimized by optimising the fibre orientations for different load cases. Therefore, the effect of transverse shear forces is considered in this study. Simulated annealing (SA) is used to search for the optimal design. This optimisation algorithm has been shown to be reliable as it is not based on the starting point, and it can escape from the local optimum points. In this research, the Tsai-Wu failure and maximum stress criteria for composite laminate are chosen. By applying two failure criteria at the same time the results are more reliable. Experimentally generated results show a very good agreement with the numerical results, validating the simulated model used. Finally, to validate the methodology the numerical results are compared to the results of previous research with specific loading.


Introduction
The demand for high strength, high modulus and low density industrial materials has generated an increased number of applications for fibre laminated composite structures in many different fields such as in submarines, sport equipment, medical instruments, civil engineering, enabling technologies, primary and secondary marine and aerospace structures, astronavigation and many more industries [1].Composite constructions are usually multilayer produced structures, mostly made of flat and curved panels, built up from several layers or laminae, which are bonded together [2].
In the last half century, the use of composite materials has grown rapidly.These materials are ideal for structural applications that require high strength and low weight.They have good fatigue characteristics and are resistant to corrosion.They provide some flexibility in design through the variation of the fibre orientation or stacking sequence of fibre and matrix materials [3,4].Another advantage of fibre laminate composites is the capability to design the physical structure and mechanical properties prior to manufacture.The mechanical behaviour of laminates strongly depends on the orientation of fibres and thickness of lamina.Accordingly, the lamina should be designed to satisfy the specific requirements of each particular application in order to obtain the maximum advantages from the directional properties of materials.Accurate and efficient structural analysis, design sensitivity analysis, and optimisation procedures for size and shape and the orientation of fibres within the material are also required.This provides a good opportunity to tailor the material properties to the specific application [5,6].However, it increases the complexity of the design problem.This complexity exists, not only because of numerous design variables, but also because of having a multimodal and variable-dimensional optimisation problem with unattainable or costly derivatives [7].
Optimum strength designs of continuous fibre-reinforced composite laminates have been used since the early days of these materials.The first research to investigate the fibre orientation of a unidirectional lamina yielding maximum strength under in-plane stress conditions has 2 Advances in Materials Science and Engineering State of stress at a point of a continuum [36].
been carried out by Sandhu [8].Brandmaier found that the strength of a unidirectional lamina under in-plane stresses could be maximised analytically with respect to the fibre orientation [9].The results based upon Tsai-Hill failure criterion indicated that the optimum fibre orientation depended upon the stress state and the relative value of the transverse and in-plane shear strengths of the lamina.When the strength of a multidirectional composite laminate is to be maximised, more complicated and explicit optimisation techniques are needed [10].Chao et al. were the pioneers that sought the optimum strength design of multidirectional laminates using a search technique [11].Latterly, many studies have been devoted to the optimum strength design of multidirectional laminate.Among these are the works by Park [12], Fukunaga and Chou [13], Miravete [14], Fukunaga and Vanderplaats [15], Gurdal et al. [16,17], Spallino et al. [18], Weaver [19], Chattopadhyay et al. [20], Luersen and Riche.[21], and Ghiasi et al. [22].
In this paper, the thickness of laminated composite plates is minimised by optimising the fibre orientation angles for different load cases.The novelty of the research presented in this paper is that the effect of transverse shear forces, and therefore, the induced twist angle are considered.2.1.Analysis of a Laminate Composite Plate.The state of stress at a point in a general continuum can be represented by nine stress components σ i j (i, j = 1, 2, 3) acting on the sides of an elemental cube with sides parallel to the axes of a reference coordinate system (Figure 1).

Governing Equation
In the most general case the stress and strain components are related by the generalised Hook's law as follows: or where C i jkl is the stiffness components [30].By considering the symmetry of the stress and strain tensors and the energy relations, it is proven that the stiffness matrices are symmetric.Thus, the state of stress (strain) at a point can be described by six components of stress (strain), and the stress-strain equations are expressed in terms of 21 independent stiffness constants [36].

In-Plane Stress.
The simplest equivalent single-layer (ESL) laminated plate theory, based on the displacement field, is the classical laminated plate theory (CLPT) [37][38][39][40][41][42][43].The two-dimensional classical theory of plates was initiated by Kirchhoff [44] in the 19th century, and then was continued by Love [45] and Timoshenko [46] during the early 20th century.The principal assumption in CLPT is that normal lines to the midplane before deformation remain straight and normal to the plane after deformation.Although this assumption leads to simple constitutive equations, it is the main deficiency of the theory.The effect of the transverse shear strains on the deformation of the elastic two-dimensional structure are ignored and some of the deformation mode constraints by reducing the model to a single degree of freedom (DoF) results are neglected.This is a consequence of the basic assumptions made.It is also worth mentioning that neglecting shear stresses leads to a reduction or removal of the three natural boundary conditions that should be satisfied along the free edges.These boundary conditions being the normal force, bending moment, and twisting couple [47].
For solving in-plane stress normally the classical laminate theory is used.It is assumed that plane stress components are taken as zero.With respect to the coordinate system shown in (Figure 2) the in-plane stress components are related to the strain components as where k is the lamina number, Q i j are the off-axis stiffness components, which can be explained in terms of principal stiffness components, Q i j , using the tensor transformation rules [48] as × sin 2 θcos 2 θ + Q 22 sin 4 θ, × sin 2 θcos 2 θ + Q 22 cos 4 θ, The principal stiffness terms, Q i j , are related to elastic properties of the material along the principal directions, E 1 , E 2 , G 12 , ϑ 12 , and ϑ 21 [48].The effect of transverse shear stress is not considered by previous work [35] because in their work the laminate is only subject to in-plane loads.Therefore, strain components defined with respect to x-y-z coordinate system are the same for each ply regardless of the fibre orientation.For the same reason, the mechanical response of the laminate is independent of the stacking sequence.Stress resultants, or forces per unit length of the cross section, are obtained as ( Here, m is the number of distinct laminae, n k is the number of plies in the kth lamina.Here, lamina is meant to be a group of plies with the same orientation angle.Substituting the stress-strain relation given by ( 3) into ( 5) we have where A i j , components of extensional stiffness matrix, are given by Principal stress components can be obtained using the following transformation [35]: Despite its limitations, CLPT is still a common approach utilised to determine quick and simple predictions especially for the behaviour of thinplated structures.The main simplification is that three-dimensional thick structural plates or shells are treated as two-dimensional plates or shells located through midthickness, which results in a significant reduction in the total number of variables and equations, consequently saving a lot of computational time and effort.
The governing equations are easier to solve and present in closed-form solutions, which normally provides more physical or practical interpretation.This approach remains popular as it is well-known and has become the foundation for further composite plate analysis methods.This method works relatively well for structures that are made out of a symmetric and balanced laminate, experiencing pure bending or pure tension.The error induced/introduced by neglecting the effect of transverse shear stresses becomes trivial on or close to the edges and corners of thick-sectioned configurations.The induced error increases for thick plates made of composite layers, for which the ratio of longitudinal to transverse shear elastic modulii is relatively large compared to isotropic materials [49].It neglects transverse shear strains, underpredicts deflections and overestimates natural frequencies and buckling loads.

Out-Plane Stress.
As discussed in Section 2.2, in the classical lamination theory, it was assumed that the laminate is thin compared to its lateral dimensions and that straight lines normal to the middle surface remain straight and normal to that surface after deformation.Therefore, the transverse shear stress and strain are neglected.These assumptions are not valid in the case of thicker laminates and laminates with low stiffness central plies undergoing significant transverse shear deformations.In order to overcome these limitations several theories have been proposed to analyse thicker laminated composite plates in order to consider the transfer shear effect.Most of these theories are extensions of the conventional theories developed by Reissner [47] and Mindlin [50], which are known as the shear deformation plate theories.These theories are based on the assumption that the displacement w is constant through the thickness while the displacements u and v vary linearly through the thickness of each layer (constant cross-sectional rotations w x and w y ).Generally these theories are known as first-order shear deformation theories (FSDT) [3].According to this theory, transverse straight lines before deformation will still be straight after deformation but they are not normal to the midplane after deformation.This theory assumes constant transverse shear stress.
In the following, referred to as first-order shear deformation laminate plate theory, the assumption of normality of straight lines is removed compared to CLPT.On the other hand straight lines normal to the middle surface remain straight but not normal to that surface after deformation [43].
where the components of this section stiffness matrix are given by 2.4.Failure Criteria.As it is shown in previous research [35] using one failure method is not reliable enough for evaluating the results.Also in each fibre orientation the type of the failure is switched between fibre, shear, and transverse shear failure [43].Therefore, it is logical that the results would be checked by at least two failure criteria.In this research maximum stress and Tsai-Wu criteria are chosen.

Maximum Stress Criterion.
Maximum stress criterion is one of the simplest failure methods to apply.According Advances in Materials Science and Engineering 5 to this criterion, failure is predicted whenever one of the principal stress components exceeds its corresponding strength.It is expressed in the form of the following subcriteria: F 1t and F 1c are the longitudinal tensile and compressive strengths, F 2t and F 2c are the transverse longitudinal tensile and compressive strengths, and F 6 is the in-plane shear strength.Four additional lamina strength parameters, which are relevant in three-dimensional analysis, are the out-plane or interlaminar tensile, compressive, and shear strengths, F 3t , F 3c , F 4 , and F 5 .

Tsai-Wu Criterion.
The Tsai-Wu failure criterion is one of the most reliable static failure criteria as it provides a simple analytical expression taking components.Reddy [43] proposed a modified tensor polynomial theory by assuming the existence of a failure in the stress space.In contracted notation it takes the form: where f i and f i j are second-and fourth-order strength tensors, and i, j = 1, 2, . . ., 6 [36].By applying assumptions, some of f i and f i j are identified.Finally it is reduced to a failure envelope for constant values of shear stress τ 6 = kF 6 or equivalently (16)

Optimisation
An optimised composite laminate requires finding the minimum number of layers, and the best fibre orientation and thickness for each layer.Several optimisation methods have been introduced to solve this challenging problem, which is often nonlinear, nonconvex, multimodal, and multidimensional.Nowadays usually stochastic non-linear optmisation methods are utilised for this problem as they can avoid the local minimums.One of the best algorithms in this category is simulated annealing (SA) method which is used in similar problems [7].
3.1.Simulated Annealing Algorithm.Kirkpatrick et al. [51] proposed simulated annealing as a powerful stochastic search technique in 1983.The method gets its name from the physical process whereby the temperature of a solid is raised to a melting point, where the atoms can move freely, and then slowly cooled.The method attempts to model the behaviour of the atoms in forming arrangements in solid material during annealing.Although the atoms move randomly, as their natural behaviour they tend to form lower-energy configurations [52].However, this is a time-driven process.When a material is crystallised from its liquid phase, it must be cooled slowly if it is to assume its highly ordered, lowestenergy, perfect crystalline state.At each temperature level during this annealing process, the material should reach equilibrium.As the temperature decreases, the arrangement of the atoms gets closer and closer to the lower energy state.This process continues until the temperature finally reaches freezing point [52].The temperature is initially assigned a higher value, which corresponds to more probability of accepting a bad move and is gradually reduced by a user-defined cooling schedule.Retaining the best solution is recommended in order to preserve the good solution [52].
At each iteration of the simulated annealing algorithm, a new point is randomly generated.The distance of the new point from the current point, or the extent of the search, is based on a probability distribution with a scale proportional to the temperature.The algorithm accepts all new points that lower the objective, but also, with a certain probability, points that raise the objective.The algorithm avoids being trapped in local minima, by accepting points that raise the objective, and is able to explore globally for more possible solutions.An annealing schedule is selected to systematically decrease the temperature as the algorithm proceeds.As the temperature decreases, the algorithm reduces the extent of its search to converge to a minimum.
If a set of configurations is considered, in each iteration the speed convergency would be increased.In this paper the SA proposed by Erdal and Sonmez [52] is applied.The number of these configurations depends on the dimension of the problem: where n is the dimension of the problem, that is, the number of design variables [52].
Advances in Materials Science and Engineering

Penalty Function and Optimisation
Procedure.In this step a penalty function is expressed, and then this function has to be optimized: n k is the number of plies in the kth lamina, in which the orientation angle is θ k ; m is the total number of distinct lamina; the second and third terms represent the penalty values introduced to increase the value of the objective function for designs for which failure is predicted, and thus to restrict the search to the feasible design space; P MS and P TW are penalty values calculated based on the maximum stress criterion and the Tsai-Wu criterion, respectively.SF MS and SF TW are equal to the safety factors according to the maximum stress, and Tsai-Wu criteria, respectively, if they are greater than 1.0, otherwise these terms are equal to zero; w i are suitable coefficients [35].This penalty function is the same as the one defined by Akbulut and Sonmez [35] in 2008 except the last term of (18).In their work, the ply angles are optimised for inplane plate and the effect of shear stress and as a result induced twist angle of the plate was neglected.However, in this research the induced twist angles are considered.Therefore the maximum acceptable twist angle is defined, and by assuming this maximum twist angle, the appropriate coefficient for w 3i is obtained.
Implying several tests by finite element method (FEM) software shows that maximum twist for each material happens at the specific angles (θ max ) i i = 1, 2, 3, . ..,i is the number of possible θ max .The range of fibre orientations [−90, 90] is divided into several areas, each θ max is the centre of the area and in each iteration the program find the θ k to the range that it belongs to and then the program works with the appropriate θ max and corresponding w 3i .Except for the related w 3i , the other w 3i are equal to zero.Δθ is defined as So, when the ply angle in each layer (θ k ) would be close to θ max , the amount of cos(Δθ) is bigger, and thus the last term of penalty function and as a result penalty function would be larger.
For a proper w 3 , the amount of induced twist (α) always will be less than maximum acceptable twist (α max ).For example, in a composite product, the final induced twist is desired to be less than a certain amount of twist (which here it is called α max ).When fibre orientation for layers approach the θ max (where maximum twist happening), the term Δθ or (θ max − θ k ) in ( 18) tends to zero.Then the term cos(Δθ) and as a result penalty function in (18) are increased.Increasing penalty function means that this fibre orientation is not going to be selected in optimisation process and therefore the induced twist in final product stays less than α max .This technique is what distinguished this model from others as the induced twist is ignored in them.
The reason that the objective is reduced for safe designs is that there may be many feasible designs with the same minimum thickness.Of these designs, the optimum is defined as the one with the largest failure load.Accordingly, the objective function is linearly reduced in proportion to the failure margin [53].Similarly in another study [54], the margins to initial failure were maximized with the minimum feasible number of laminae.The safety factor of the laminate according to the maximum stress criterion, SF MS , is calculated as follows [35]: The safety factor for the kth lamina, SF k TW , according to the Tsai-Wu criterion is defined as the multiplier of the stress components at lamina k, σ k i j , that makes the right hand side of ( 16) equal to 1 then it turns into [35] where The root of the above equation gives the safety factor.Because a negative safety factor is not physically meaningful, the absolute value of the first root is considered as the actual safety factor Then, the minimum of SF k TW is chosen as the safety factor of the laminate: In (24), the |(−b + √ b 2 + 4a)/2a| can be considered, as b is always positive, and the aim is to find the minimum of SF k TW .The penalty value due to the violation of the maximum stress and Tsai-Wu criteria are calculated in ( 26) and ( 27), respectively: The total penalty value for the laminate due to the violation of the maximum stress and Tsai-Wu criteria are then calculated by summing up the penalty values calculated for each lamina (28)

Experimental Results.
In order to validate the FEM model some experimental tests have been performed.For each case six similar laminated plates are manufactured.One of the samples is shown in Figure 3. Carbon fibre is used for all laminated experimental tests, and the size of plate is 500 * 500 (mm).Square laminated composite panels (500×500 mm) were produced by vacuum bagging.In order to minimise the error, the average results of these six plates were calculated to be compared with FEM results obtained from the design tool.As it is shown in Figure 4, the plates are fixed on one side (a) and displacement in z direction is measured on the opposite free side of the plate (b) due to the effect of gravity only.
In Figure 5, t 1 to t 6 are the experimental test results for six similar plates under the constant load.Six plates in each case have the same layup and geometry; t a is the average of t 1 to t 6 and t F is the FEM results.
The experimental tests have been carried out for 15 different cases with different loads and layups.In Figure 6 the percentage difference of deformation between tests and FEM results and in Figure 7 the percentage difference of α are shown.As is shown in Figures 6 and 7 there is a very good agreement between FEM and experimental results.In Figures 6 and 7 the axis which shows error is zoomed in to 10% to distinguish the differences between each case.These Experimental tests have been done to validate the FE model which is used in this work.The good agreement, as it is shown in Figures 6 and 7, indicates the validity of this FEM.In the following case studies, FE results are compared with In this case the maximum acceptable twist angle is α max = .01.In Tables 1 and 2 the results are compared with previous work [35].As it is shown in Table 1, the number of layup increase for some loads and it shows that in these cases the amount of twist angle is more than α max .Clearly by increasing the lay-ups, safety factor for both Tsai-Wu and Maximum stress will be increased.In Table 2, the number of lay-ups and thickness is constant but in some cases the optimum orientation angles is different from the previous work.Although the safety factor is reduced, it avoids passing the acceptable twist angle and the safety factor is still more than 1, therefore, it is still acceptable.When the maximum twist angle is less than the α max , the results are comparable with the work has been done by Akbulut and Sonmez [35].
In the second case study a highly anisotropic material is considered (material II).The elements of stiffness matrix are In Figure 8 the amount of twist angle under a constant load for different ply angles is shown.So the amount of (θ max ) i in ( 20) are (±30, ±60) which are the local maximums.This test is for the first layer in order to find the (θ max ) i and the areas which were explained in Section 3.2.FEM tests show that these (θ max ) i are the same by adding the next layers.If the stiffness matrix is symmetric the curves will be symmetric about y axis.A general stiffness matrix for material II is considered, so there is no mirror about y axis.Optimum lamina orientations under different loads in this case are shown in the Table 3.In this case the pure bending load N zz is also applied to the plate.

Conclusion
In this study, an optimisation methodology of composite plates was presented.A method was proposed in order to overcome the difficulties and shortcomings faced by the previous research.In previous work the effect of transverse shear was neglected, and therefore the induced twist angle is ignored.In some applications the twist angle, which is the direct effect of transverse shear, is undesirable.Therefore, in this research, after optimising the fibre orientations, by considering the induced twist angle as well as safety factor, the induced twist angle always stays less than the acceptable twist angle.One of the other weakness in previous work was that the plate was optimised under specific loads, such as longitudinal or in-plane loading.By the proposed method in this research the out-plane stress optimisation can be solved as well as the in-plane stresses.In order to have a reliable optimisation, simulated annealing, which is one of the stochastic optimisation methods and can escape the local minima is applied and the penalty function for this optimisation method is modified.This modified penalty function forces the induced twist to stay under a predefined induced twist.In addition, two Tsai-Wu and maximum stress failure criteria are used in the algorithm individually to avoid false optimal design.

Figure 1 :
Figure1: State of stress at a point of a continuum[36].

Figure 3 :
Figure 3: Process of making laminated plate.

Figure 4 :
Figure 4: Boundary condition and loading (the plate is clamped on side a. Displacements are measured on side b).

Figure 8 :
Figure 8: Twist angle for a first layer of material (II).
12 A 13 B 11 B 12 B 13 0 0 A 12 A 22 A 23 B 12 B 22 B 23 0 0 A 13 A 23 A 33 B 13 B 23 B 33 0 0 B 11 B 12 B 13 D 11 D 12 D 13 0 0 B 12 B 22 B 23 D 12 D 22 D 23 0 0 B 13 B 23 B 33 D 13 D 23 D 33 0 0 Figure5: An example for one case; t 1 to t 6 are the experimental test results for six similar plates under constant load.t a is the average of t 1 to t 6 and t F is the FEM outputs.

Table 1 :
[35]mum lamina orientations for material T300/5308 under different loads.In this Section, two case studies are considered.The first case study compares the obtained results with those which were found by Akbulut and Sonmez[35]to validate the proposed optimisation method.The graphite/epoxy materials T300/5308 with the properties of E 11 = 40.91GPa, E 22 = 9.88 GPa, G 12 = 2.84 GPa, ϑ 12 = 0.292, F 1t = 779 MPa, F 1c = −1134 MPa, F 2t = 19 MPa, F 2c = −131 MPa, and S = 75 MPa is considered for the first case study.

Table 3 :
Optimum lamina orientations for second case study material under different loads for constant thickness.Loading: N xx , N yy , N xy , N zz (KPa m)