Thermal Aeroelastic Tailoring for Laminated Panel with Lamination Parameters and JAYA Algorithm

This work investigates thermal aeroelastic tailoring of a laminated composite panel with a lamination parameter-based method. Equivalent membrane and bending coefficients of thermal expansions for symmetric laminated panel are derived and represented with lamination parameters using Classical Laminated Plate Theory. The relationship between thermal flutter behavior and lamination parameters is examined. The optimization process is split into two stages. In the first stage, lamination parameters and laminate thicknesses are as design variables to minimize the structure mass, subject to thermal flutter behavior and feasible region constraints of lamination parameters. In the second-stage, instead of using conventional genetic algorithm, the enhanced JAYA method is extended to search the laminate configuration to target the optimal lamination parameters. The effectiveness of the presented method is demonstrated through a thermal aeroelastic tailoring problem.


Introduction
Panel thermal flutter is a self-excited oscillation phenomenon resulting from the interactions of the inertial force, elastic force, the aerodynamic pressure, and the induced temperature load due to supersonic airflow [1]. Usually, when the Mach number of flight vehicle is greater than 3, the aerodynamic heating effect on panel structure will become obvious. The thermal load generated by boundary constraint during thermal expansion can reduce the structural stiffness and produce the so-called stiffness "softening" phenomenon. Mechanical stiffness softening can directly change the dynamic characteristics and usually reduce the thermal flutter boundary of the panel.
Fiber reinforced laminated composite materials have been increasingly used in aerospace for skin panels of aircraft fuselage and wing to reduce the structural weight. Increasing attentions have been paid to thermal flutter investigation for laminated composite structures. Ganapathi and Touratier [2] presented an investigation of the effects of various temperature loading, laminate thickness, and boundary conditions on the thermal flutter behavior. Supersonic flut-ter analysis of stiffened laminated plates has been presented by Lee et al. [3]; results showed that a proper stiffening method can improve the flutter behavior for stiffened laminated panels. Recently, Zhao and Cao [4] investigated the flutter characteristics of a laminated panel with aerodynamic, thermal, and acoustic loads. Chai et al. [5] investigated the flutter and thermal buckling behavior of laminated shells with various boundary conditions. Huang et al. [6] detailed a parametric study of the relationship between the number of plies and the flutter characteristics of laminated plates. With the development of advanced manufacturing technologies such as automated fiber placement (AFP), fiber angles can be placed continuously and accurately in any direction theoretically [7]. Also, thermal flutter characteristics of curvilinear fiber path variable stiffness composite has been paid attention recently [8][9][10].
Both active control method and passive control method can be used to improve thermal flutter behavior [11]. Smart materials [12,13], advanced control techniques [14], and viscoelastic materials [15] are used for flutter control. Aeroelastic tailoring design is also a passive form for thermal flutter avoidance. Aeroelastic tailoring is a structural design strategy in which the arrangement of mechanical properties leads to more control on the static and dynamic deformations, thereby impacting on the aeroelastic performance of the aircraft [16]. There has been many researches on the application of composite structure for aeroelastic tailoring design [17][18][19][20][21][22]. Not too many research has been encountered for the specific task of aeroelastic tailoring for panel thermal flutter avoidance in the literature. Marques et al. [23] used layer angles as design variables to maximize the thermal flutter dynamic pressure for various stiffened laminated panels, using conventional genetic algorithm (GA). Results showed that the flutter behavior can be improved by altering the flutter mechanism to passively avoid flutter instabilities. However, excessive number of structural and aeroelastic analysis are required using GA as optimization algorithm.
Lamination parameters (LPs), first presented by Tsai and Pagano [24], can be used as design variables instead of conventional layer angles for laminate optimization. The stiffness matrix (A, B, D) of laminated panel can be expressed as a linear function of the LPs [7]. LPs seem to be efficient design variables for laminate configuration defined by a convex feasible region, which makes them suitable as design variables in gradient-based optimization algorithms for composite structure design problem [25]. The polar parameter (PP) modeling method proposed by Vannucci and Verchery firstly can also be used for laminated panel modeling [26,27], and a multiscale two-level optimization strategy (MS2LOS) based on PPs has been developed by Montemurro et al. [28].
The LP-based optimization procedure is usually split into two stages. In the first stage, gradient-based optimization algorithm is adopted with LPs and laminate thickness as design variables. In the second stage, logic-based methods [29] and GA [30][31][32][33] are usually used to find the real laminate configuration to match the optimal LPs obtained in the first stage.
Lots of metaheuristic algorithms with various inspiration mechanism have been developed in recent years. In this paper, we have focused on the JAYA algorithm proposed by Rao [34], to be used in the second optimization stage for layup match. JAYA algorithm is inspired by a simple concept that the solution should move towards the best solution and away from the worst solution. JAYA and its variants have been used for many engineering problems in various fields [35][36][37][38][39][40]. Recently, Zhang et al. [41] proposed an enhanced JAYA algorithm (EJAYA) for global optimization. Both results of numerical optimization and engineering design problems demonstrate the efficiency of the improved strategies.
This paper presents a study on thermal aeroelastic tailoring of a symmetric laminated composite panel in supersonic airflow. Laminated panel finite element model considering thermal load constructed with lamination parameters and a ZONA51 supersonic lifting surface theory are used to predict the flutter boundary. With parametric study, thermal flutter behavior relationships with LPs, a two-stage layup optimization strategy is employed for thermal aeroelastic tailoring. In the second stage, the enhanced JAYA algorithm is extended to search the real layups for the optimum LPs obtained in the first stage. Section 2 of this paper outlines the mathematical model of panel thermal flutter. Section 3 presents the proposed method for thermal aeroelastic tailoring; Section 4 describes the relationships between thermal flutter behaviors and lamination parameters; thermal aeroelastic tailoring results are also analyzed; brief conclusions are given in Section 5.

Mathematical Model
2.1. Boundary Conditions for Thermal Flutter. Most published studies have focused on a clamped or simply supported boundary conditions for flutter and buckling investigation of laminated panel in a supersonic flow, and relaxed boundary condition is usually considered by boundary constraint springs [5]. In this paper, the flexible boundary condition similar as TPSS (Thermal Protection System Support) of X-33 TPS (shown in Figure 1) [42] and ARMOR (Adaptable, Robust, Metallic, Operable, Reusable) TPS (shown in Figure 2) [43] is investigated. Four metal support bracket are used to connect and support the laminated panel, which can accommodate the thermal expansion and prevent buckling [43].
The mathematical model considers a rectangular panel (a × b × t) with four elastic brackets supported under supersonic airflow, shown in Figure 3. The high-speed flow field is supposed at the upper (external) side of the panel, while stagnated air remains in a cavity representing the bottom (internal) side of a vehicle [44].     Usually, symmetric laminates about the midplane are commonly used to avoid bending-extension coupling in aerospace engineering [45]. In the symmetric laminate without extension-shear and bending-twisting coupling (V B i = 0), and if the layer angles are limited to 0°, 90°, +45°, and -45°, the in-plane and out-of-plane elastic properties are governed by 3 membrane and 3 bending LPs (V A 4 = V D 4 = 0). Fukunaga et al. [46] firstly presented the relationships between thermoelastic properties and LPs. According to Refs. [46,47], equivalent thermal properties of symmetric laminates are determined as follows.
The following coefficients are used to determine equivalent thermal properties [48].
where fαg k is the thermal expansion coefficient vector depending on the thermal expansion coefficients α L and α T of the lamina, The membrane and bending coefficients of thermal expansions (CTEs) are determined as follows: Membrane CTEs can be expressed with LPs.
Bending CTEs can be represented similarly.

Unsteady Aerodynamics
Modeling. ZONA51 supersonic lifting surface theory is used for unsteady aerodynamics modeling. This panel uses ZONA51 theory, and 15 × 15 equal boxes are divided for the panel surface (shown in Figure 4(b)). Two regions to either side of the panel are also modeled, so that the center panel behaves as a panel within a surface rather than as a wing [49]. SPLINE4 cards for transferring of loads and displacements between aerodynamics model to the mechanical one.

Solution Method.
Three assumptions are usually adopted during thermal flutter analysis: (1) the temperature field is not affected by deformation of the structure; (2) the time scale of flutter dynamic response is much smaller than that of temperature change, so the temperature field can be regarded as a steady state during thermal flutter analysis; (3) when the temperature rise is not too high, the influence of temperature on material properties is ignored. The solution of thermal flutter analysis in this paper is divided into two processes. Firstly, finite element model of laminate panel with LPs considering thermal load is established based on the method detailed in section B. Through steady-state nonlinear heat conduction analysis (Solution 153 in MSC.NASTRAN), the additional stiffness of the structure caused by thermal load is introduced into the overall stiffness of the structure using RESTART statement in MSC.NASTRAN, and the equivalent stiffness matrix is obtained. Secondly, coupling with ZONA51 supersonic lifting surface theory, PK-method is used to solve the flutter equation of heated laminated panels in supersonic flow (Solution 145 in MSC.NASTRAN).

Two-Stage Optimization Procedure
Although the use of FEA-based flutter analysis in aeroelastic tailoring is widely used, the optimization process is time consuming. A two-stage optimization procedure using LPs as design variables for laminated panel is employed for thermal aeroelastic tailoring to reduce optimization time.
3.1.2. Optimization for Lamination Parameters. One of the best nonlinear programming methods used for solving the constrained optimization problems is the sequential. Sequential quadratic programming (SQP) is used in the first stage, which is one of the best nonlinear programming methods for the nonlinear problems (NLPs). The NLP is approximated by a quadratic program (QP) problem during optimization process. VMCWD subroutine by Powell [53] is used in the first-stage optimization, which included the watchdog technique proposed by Chamberlain et al. [54].
To minimize W X f g ð Þ, Subject to V f ≥ V req and Eq: 18 The constraints are all transformed as where dc j ðfXgÞ is the jth design constraint and nc is the number of constraints. The penalty function method is used as follows [38]: where ν denotes the total sum of violated constraints; ε 1 and ε 2 are two parameters to make a good balance between the exploration and exploitation rate. ε 1 is set to 1, and ε 2 starts at 1.5 and then increases linearly to 3 at the last step, during the search procedure [38].
The flutter speed is calculated by MSC.NASTRAN through the procedure described in Section 2.4. This solution is used within the optimization process aiming at min-imizing the penalty function as Equation (19). The sensitivities of the penalty function can be obtained by the difference method.

Second-Stage Optimization
3.2.1. JAYA Algorithm. JAYA is a population-based algorithm proposed for optimization problems [34], which is easily adaptable because of not using any algorithm specific parameters. The solution updating procedure used in the basic JAYA algorithm is as follows [41]: where N is the population size, x i is the ith solution in the population, and x best and x worst are the best and worst   The enhanced JAYA algorithm recently proposed by Zhang et al. [41] is used to search the laminate configuration to match the optimal LPs. Lower and upper local attractors are used for local exploitation strategy; historical population is used for global exploration strategy. This scheme is briefly presented below.
Compared with basic JAYA algorithm, additional parameters of EJAYA are as follows: M, the current mean solution, is calculated by P u , a proposed upper local attract point parameter to represent the solution between the current best solution and the current mean solution, is calculated by P l , a proposed lower local attract point parameter to describe the solution between the current worst solution and the current mean solution, is calculated by The historical population X old ðX old = fx old,1 , x old,2 , ⋯, x old,N gÞ was first generated by Then, X old is processed by X old = permutingðX old Þ, where permuting is a random shuffling function, which is to sort all individuals of X old in random order [41]. The local exploitation strategy of EJAYA can be described as The global exploration strategy can be described as where λ 3 , λ 4 , λ 5 and λ 6 are all random numbers in [0, 1].

Layup
Design with Enhanced JAYA Algorithm. The objective function Dis is calculated by the difference between the optimal LPs and the actual LPs of the real laminate configuration.  Suppose the stacking sequence is symmetric with respect to the center of the laminate, hence just half number of the plies need to be optimized. The layer angles are treated as continuous variables. The discrete value of the variable is chosen from the interval ð−90+,90 + Δα, where Δα represents the minimal angle increment between allowable layer angles. For the case where the layer angles are limited to 0, 90, and ±45, Δα is 45 [55]. P select between 0 and 1 is used in EJAYA procedure to selected probability for local or global exploration strategy. If P select > 0:5, local exploitation strategy is selected. The local exploitation strategy will be performed by Equation (21), Equation (22), Equation (23), and Equation (25). If P select ≤ 0:5, the global exploration strategy is selected by Equation (24) and Equation (26).

12
International Journal of Aerospace Engineering on the D matrix, that is, just relates to If the layer angles are limited to 0°, 90°, +45°, and -45°, then V D 4 = 0. Therefore, the influence of changes of V D 1 , V D 2 , V D 3 on flutter speed is studied firstly. Figure 5 shows the contours of flutter speed on V D 1 − V D 2 plane with V D 3 = 0. V D 3 = 0 means that both D 16 and D 26 are 0, there is no bending-torsion coupling, and the laminate presents an orthotropic characteristic. It can be seen that when V D 2 is greater than -0.4 (region I), the change of flutter speed mainly depends on the variation of V D 1 . The maximum critical speed approximately appears at V D 1 = 0. The flutter speed shows a decreasing trend with either increase or decrease of V D 1 . While when V D 2 is less than -0.4 (region II), the change of flutter speed mainly depends on variation of V D 2 . The switch of the flutter modes results in the discontinuity of the contours of flutter speed between the two regions, as mentioned in Ref. [56]. Figure 6 lists the first six modes of the panel. Compared with the cantilever plate, the first four modes can also be named as first-order bending mode, first-order torsion mode, second-order antisymmetric bending mode, and second-order symmetric bending mode, respectively. In region I, the flutter mode is the bending-torsion coupling mode, dominated by the second-order antisymmetric bending (mode 3) and participated with the first-order torsional mode (mode 2). In region II, the flutter mode is dominated by second-order antisymmetric bending (mode 3) and participated with second-order symmetric bending (mode 4). Two points (point A and point B) are selected from the two regions, respectively; v-g and v-f curves are shown in Figures 7 and 8. For point A, the flutter speed is 3121.2 m/ s due to coupling between mode 3 and mode 2. For point B, the flutter speed is 528.9 m/s due to coupling between mode 3 and mode 4. 14 International Journal of Aerospace Engineering which will affect the dynamic characteristics and aeroelasticity of the panel structure. Additional thermal compression load is related to the temperature variation, in-plane stiffness, and in-plane thermal expansion coefficients. Larger CTEs will lead to greater thermally induced load and more obvious influence on the dynamic characteristics of the panel. The in-plane stiffness and in-plane thermal expansion coefficient can be expressed as functions of membrane LPs.
on CTEs is investigated. The CTE in the x direction is recorded as CTE_x, and CTE_y for y direction. In the feasible region of LPs, discrete points of thermal expansion coefficient are shown in Figure 10. The minimum value of CTE_x is 8.73e-07 at ½V A 1 , V A 2 , V A 3 = ½0:5, 0:4, 0, and the minimum value of CTE_y is 8.73e   Figures 11 and 12 show the contours of CTEs on Figure 13 shows the contours of CTEs

Panel Flutter Behavior with Thermal Effect Consideration.
As mentioned above, membrane CTEs influence the dynamic characteristics of a panel when subjected to thermal load. Two points with different CTEs and identical bending stiffness are selected for panel flutter behavior investigation.
For point C, Figure 14 describes the frequency change curves of the first six models with temperature increase. It can be seen that the mode shape of each order is consistent with nonthermal modes, and frequency decreases slightly with temperature increase. Figure 15 gives the frequency change curves for point D. It can be seen that when the temperature rises by 20 degrees, the first two modes exchange compared with nonthermal modes. The frequency keeps increasing with temperature increase. The larger coefficient of thermal expansion, the greater influence of temperature on frequency change for the panel.
The relative change of each order frequency will also affect the flutter speed. For point C, the flutter speed change is very small, reaches the maximum at 20 degrees, and then, decreases gradually, shown in Figure 16(a). The flutter mode is bending-torsion coupling flutter dominated by secondorder antisymmetric bending (mode 3) and participated by first-order torsional mode (mode 2). For point D, the flutter speed change is shown in Figure 16(b). The discontinuity of flutter speed change is also due to the switch of flutter mode.

First-Stage Optimization
Results. The laminate thickness and LPs are as design variables for weight reduction using SQP, with flutter speed limit (≥3000 m/s) and feasible regions of LPs as constraints. As mentioned above, when the discontinuity of flutter speed due to the switch of the flutter modes occurs in design space, optimization procedures may break down at the point where the flutter speed changes discontinuously [56]. Different initial values should be attempted for this optimization procedure.
The optimal result is h = 0:0019032 m, and  Figure 17 details the frequencies and mode shapes of first six modes of optimal solution. V-g and v-f curves are shown in Figure 18. The flutter mode is the bending-torsion coupling flutter mode, dominated by the second-order antisymmetric bending (mode 3) and participated with the first-order torsional mode (mode 2).

Second-Stage Optimization Results.
In the second optimization stage, the optimal LPs obtained in the first stage are used as target LPs. The ply thickness is 0.1 mm in this problem, and the ply number is rounding up to the nearest larger integer number 20. Although this will not result in the lightest weight solution, it will always ensure a safe solution [51].
Two different kinds of discrete design variables are used to target the optimum LPs. Optimal 1 describes the layup when the layer angles are restricted as 45°, -45°, 0°, and 90°. Optimal 2 gives the optimal laminate configuration when the layer angles are extended to 15-degree intervals between -75°and 90°. The LPs related to the optimized laminate configurations together with the Dis are listed in Table 1. As can be seen, the Dis are very small, which means that the LPs of the obtained laminate configurations are good matches for the optimal LPs. And larger design space will lead to smaller Dis.
The optimized laminate configuration is checked for frequency and thermal flutter performance, shown in Table 2.  Results show comparison of these frequencies and thermal flutter speeds obtained using the real laminate configurations with those obtained using the optimal LPs directly, there is only a small increase or decrease of frequency and flutter speed caused by the mismatch between the target LPs and the real laminate configuration.

Conclusions
Based on parametric study of the relationships between thermal flutter behavior and LPs, this paper presents a two-stage thermal aeroelastic tailoring methodology using LPs. In the first stage, LPs and laminate thicknesses for symmetric laminates are as design variables using SQP for optimal design. The optimal LPs obtained are used as targets in the second stage in which the enhanced JAYA algorithm is used to search the real layups for the optimum LPs. Also, the switch of flutter modes can result in the discontinuity of the contours of flutter speed, and the membrane CTEs have an obvious effect on the frequency and flutter speed of the laminated panel when subjected to thermal load.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.