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

Authors: Peng Jin; Qihao Zheng; Wei Wei
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

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.


---

## Body

## 1. 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 flutter 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–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–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–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–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.Section2 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.


## 2. 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].Figure 1
X-33 metallic TPS.Figure 2
ARMOR TPS.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].Figure 3
Rectangular panel model.


### 2.2. Structural Model with Lamination Parameters

According to Classical Laminated Plate Theory (CLPT), the load-strain relationship for a laminated panel can be described by(1)NM=ABBDεχ−NTMT,where N and M denote the force and moment resultants, respectively; NT and MT denote the thermal force and moment resultants, respectively; ε and χ denote the strains and curvatures at the middle plane; and A, B, and D are membrane, coupling and bending stiffness matrices, respectively.
(2)A,B,D=∑Q¯k1,z,z2dz,where Q¯k is the transformed elastic constant, defined by
(3)Q¯=UTQU.U and Q are the modulus transformation matrix and material elastic modulus matrix, respectively.
(4)U=cos2θsin2θcosθsinθsin2θcos2θ−cosθsinθ−2cosθsinθ−2cosθsinθcos2θ−sin2θ,Q=Q11Q12Q16Q12Q22Q26Q16Q26Q66.Q11=E11/1−ν12ν21, Q22=E22/1−ν12ν21, Q12=ν12E2/1−ν12ν21, Q66=G12, and ν21=ν12E22/E11.A, B, and D can be expressed by lamination parameters and the material stiffness invariants U¯.
(5)A11A22A12A66A16A26=hU¯1V1AV2AU¯1−V1AV2AU¯40−V2AU¯50−V2A0V3A2V4A0V3A2V4A1U¯2U¯3,(6)B11B22B12B66B16B26=h24V1BV2B−V1BV2B0‐V2B0‐V2BV3B2V4BV3B2‐V4BU¯2U¯3,D11D22D12D66D16D26=h312U¯1V1DV2DU¯1−V1DV2DU¯40−V2DU¯50−V2D0V3D2V4D0V3D2−V4D1U¯2U¯3.The material stiffness invariantsU¯ are calculated as follows:
(7)U¯1U¯2U¯3U¯4U¯5=18323440−401−21−4161−41−214Q11Q12Q22Q66.LPs are as shown below:(8)V1234A=1h∫−h/2h/2cos2θcos4θsin2θsin4θdz,V1234B=4h2∫−h/2h/2cos2θcos4θsin2θsin4θzdz,V1234D=12h3∫−h/2h/2cos2θcos4θsin2θsin4θz2dz,where ViA,ViB,ViDi=1,2,3,4 are membrane, coupling, and bending LPs, respectively.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 (ViB=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 (V4A=V4D=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.
(9)αεαχ=A00D−1Q¯αεQ¯αχ.The following coefficients are used to determine equivalent thermal properties [48].
(10)Q¯αε=∑k=1NQ¯kαkzk−zk−1,Q¯αχ=13∑k=1NQ¯kαkzk3−zk−13,where αk is the thermal expansion coefficient vector depending on the thermal expansion coefficients αL and αT of the lamina,
(11)αk=U−1αLαT0.The membrane and bending coefficients of thermal expansions (CTEs) are determined as follows:(12)αε=A−1Q¯αε=A−1∑k=1NQ¯kαkzk−zk−1,αχ=D−1Q¯αχ=13D−1∑k=1NQ¯kαkzk3−zk−13.Membrane CTEs can be expressed with LPs.(13)αε=A−1Q¯αε=A−1∑k=1NQ¯kαkzk−zk−1=A−1∑k=1NUkTQkUkUk−1αLαT0zk−zk−1=A−1∑k=1NUkTQkzk−zk−1=A−1∑k=1Ncos2θsin2θ‐2cosθsinθsin2θcos2θ2cosθsinθcosθsinθ−cosθsinθcos2θ−sin2θkQ11Q12Q16Q12Q22Q26Q16Q26Q66kαLαT0zk−zk−1=A−1∑k=1N1+cos2θ21−cos2θ2‐sin2θ1−cos2θ21+cos2θ2sin2θ12sin2θ−12sin2θcos2θkQ11αL+Q12αTQ12αL+Q22αT0zk−zk−1=A−1h21+V1Ah21‐V1A‐hV3Ah21‐V1Ah21+V1AhV3Ah2V3Ah2V3AhV1AQ11αL+Q12αTQ12αL+Q22αT0.Bending CTEs can be represented similarly.(14)αχ=D−1Q¯αχ=13D−1∑k=1NQ¯kαkzk3−zk−13=13D−1h381+V1Dh381‐V1D‐h34V3Dh381‐V1Dh381+V1Dh34V3Dh38V3Dh38V3Dh34V1DQ11αL+Q12αTQ12αL+Q22αT0.The structural parameters and material properties of the panel are as follows:a=0.24m, b=0.16m, t=0.002m, ρ=1580kg/m3, E1=181×109N/m2, E2=10.3×109N/m2, μ12=0.28, G12=7.13×109N/m2, αL=0.45×10−6/°C, αT=30.1×10−6/°C, where E1, E2, and G12 are elastic and shear moduli and μ12 is the Poisson’s ratio of the material. αL and αT denote the coefficients of thermal expansion in the longitudinal and transverse direction, respectively. Young’s modulus and Poisson’s ratio of the brackets are 70×109N/m2 and 0.3, respectively. And the brackets are all with lengths of 0.04 m and square sections of 0.012m×0.012m. With MSC.NASTRAN, the panel is modeled with 100 square CQUAD4 elements, and each bracket is modeled with 8 CBAR elements, with DOFs of 4157. The end points of four brackets are clamped, shown in Figure 4(a).Figure 4
Structural model and aerodynamic model.
(a)
Structure model(b)
Aerodynamic model


### 2.3. Unsteady Aerodynamics Modeling

ZONA51 supersonic lifting surface theory is used for unsteady aerodynamics modeling. This panel uses ZONA51 theory, and15×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.


### 2.4. 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).


## 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].Figure 1
X-33 metallic TPS.Figure 2
ARMOR TPS.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].Figure 3
Rectangular panel model.


## 2.2. Structural Model with Lamination Parameters

According to Classical Laminated Plate Theory (CLPT), the load-strain relationship for a laminated panel can be described by(1)NM=ABBDεχ−NTMT,where N and M denote the force and moment resultants, respectively; NT and MT denote the thermal force and moment resultants, respectively; ε and χ denote the strains and curvatures at the middle plane; and A, B, and D are membrane, coupling and bending stiffness matrices, respectively.
(2)A,B,D=∑Q¯k1,z,z2dz,where Q¯k is the transformed elastic constant, defined by
(3)Q¯=UTQU.U and Q are the modulus transformation matrix and material elastic modulus matrix, respectively.
(4)U=cos2θsin2θcosθsinθsin2θcos2θ−cosθsinθ−2cosθsinθ−2cosθsinθcos2θ−sin2θ,Q=Q11Q12Q16Q12Q22Q26Q16Q26Q66.Q11=E11/1−ν12ν21, Q22=E22/1−ν12ν21, Q12=ν12E2/1−ν12ν21, Q66=G12, and ν21=ν12E22/E11.A, B, and D can be expressed by lamination parameters and the material stiffness invariants U¯.
(5)A11A22A12A66A16A26=hU¯1V1AV2AU¯1−V1AV2AU¯40−V2AU¯50−V2A0V3A2V4A0V3A2V4A1U¯2U¯3,(6)B11B22B12B66B16B26=h24V1BV2B−V1BV2B0‐V2B0‐V2BV3B2V4BV3B2‐V4BU¯2U¯3,D11D22D12D66D16D26=h312U¯1V1DV2DU¯1−V1DV2DU¯40−V2DU¯50−V2D0V3D2V4D0V3D2−V4D1U¯2U¯3.The material stiffness invariantsU¯ are calculated as follows:
(7)U¯1U¯2U¯3U¯4U¯5=18323440−401−21−4161−41−214Q11Q12Q22Q66.LPs are as shown below:(8)V1234A=1h∫−h/2h/2cos2θcos4θsin2θsin4θdz,V1234B=4h2∫−h/2h/2cos2θcos4θsin2θsin4θzdz,V1234D=12h3∫−h/2h/2cos2θcos4θsin2θsin4θz2dz,where ViA,ViB,ViDi=1,2,3,4 are membrane, coupling, and bending LPs, respectively.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 (ViB=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 (V4A=V4D=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.
(9)αεαχ=A00D−1Q¯αεQ¯αχ.The following coefficients are used to determine equivalent thermal properties [48].
(10)Q¯αε=∑k=1NQ¯kαkzk−zk−1,Q¯αχ=13∑k=1NQ¯kαkzk3−zk−13,where αk is the thermal expansion coefficient vector depending on the thermal expansion coefficients αL and αT of the lamina,
(11)αk=U−1αLαT0.The membrane and bending coefficients of thermal expansions (CTEs) are determined as follows:(12)αε=A−1Q¯αε=A−1∑k=1NQ¯kαkzk−zk−1,αχ=D−1Q¯αχ=13D−1∑k=1NQ¯kαkzk3−zk−13.Membrane CTEs can be expressed with LPs.(13)αε=A−1Q¯αε=A−1∑k=1NQ¯kαkzk−zk−1=A−1∑k=1NUkTQkUkUk−1αLαT0zk−zk−1=A−1∑k=1NUkTQkzk−zk−1=A−1∑k=1Ncos2θsin2θ‐2cosθsinθsin2θcos2θ2cosθsinθcosθsinθ−cosθsinθcos2θ−sin2θkQ11Q12Q16Q12Q22Q26Q16Q26Q66kαLαT0zk−zk−1=A−1∑k=1N1+cos2θ21−cos2θ2‐sin2θ1−cos2θ21+cos2θ2sin2θ12sin2θ−12sin2θcos2θkQ11αL+Q12αTQ12αL+Q22αT0zk−zk−1=A−1h21+V1Ah21‐V1A‐hV3Ah21‐V1Ah21+V1AhV3Ah2V3Ah2V3AhV1AQ11αL+Q12αTQ12αL+Q22αT0.Bending CTEs can be represented similarly.(14)αχ=D−1Q¯αχ=13D−1∑k=1NQ¯kαkzk3−zk−13=13D−1h381+V1Dh381‐V1D‐h34V3Dh381‐V1Dh381+V1Dh34V3Dh38V3Dh38V3Dh34V1DQ11αL+Q12αTQ12αL+Q22αT0.The structural parameters and material properties of the panel are as follows:a=0.24m, b=0.16m, t=0.002m, ρ=1580kg/m3, E1=181×109N/m2, E2=10.3×109N/m2, μ12=0.28, G12=7.13×109N/m2, αL=0.45×10−6/°C, αT=30.1×10−6/°C, where E1, E2, and G12 are elastic and shear moduli and μ12 is the Poisson’s ratio of the material. αL and αT denote the coefficients of thermal expansion in the longitudinal and transverse direction, respectively. Young’s modulus and Poisson’s ratio of the brackets are 70×109N/m2 and 0.3, respectively. And the brackets are all with lengths of 0.04 m and square sections of 0.012m×0.012m. With MSC.NASTRAN, the panel is modeled with 100 square CQUAD4 elements, and each bracket is modeled with 8 CBAR elements, with DOFs of 4157. The end points of four brackets are clamped, shown in Figure 4(a).Figure 4
Structural model and aerodynamic model.
(a)
Structure model(b)
Aerodynamic model


## 2.3. Unsteady Aerodynamics Modeling

ZONA51 supersonic lifting surface theory is used for unsteady aerodynamics modeling. This panel uses ZONA51 theory, and15×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.


## 2.4. 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).


## 3. 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. First-Stage Optimization

#### 3.1.1. Feasible Region of Lamination Parameters

Explicit expressions for solving feasible regions of LPs have been advanced progressively by many authors [25, 50–52]. LPs are restricted by the constraints as Equation (17) for symmetric laminate with layer angles limited to 0°, 90°, +45°, and -45°.
(15)2V1X−V2X−1≤0,2V3X+V2X−1≤0,X=A,D,4ViD+rViA+r≥ViA+r4,i=1,2,3,r=±1,162V1D+rV2D+rs2V1A+rV2A+rs≥2V1A+rV2A+rs4,r=±1,s=1,−3,162V3D+rV2D+rs2V3A+rV2A+rs≥2V3A+rV2A+rs4,r=±1,s=1,−3,4V1D+rV3D+sV1A+rV3A+s≥V1A+rV3A+s4,r=±1,s=±1.


#### 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].
(16)FindX=h,V1A,V2A,V3A,V1D,V2D,V3D,,(17)TominimizeWX,SubjecttoVf≥VreqandEq.18.The constraints are all transformed as(18)dcjX≤0,j=1,2,⋯,nc,where dcjX is the jth design constraint and nc is the number of constraints. The penalty function method is used as follows [38]:
(19)fpenaltyX=1+ε1×vε2×WX;v=∑i=1ncmax0,dcjX,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 Section2.4. This solution is used within the optimization process aiming at minimizing the penalty function as Equation (19). The sensitivities of the penalty function can be obtained by the difference method.


### 3.2. 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]:
(20)vi=xi+λ1×x−bestxi−λ2×xworst−xi,i=1,2,⋯,N,where N is the population size, xi is the ith solution in the population, and xbest and xworst are the best and worst solutions in the population, respectively. Further, λ1 and λ2 are two random numbers in [0, 1]. Better solution of xi and vi will be selected into the population to accelerate convergence.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
(21)M=1N∑i=1Nxi,i=1,2,⋯,N.Pu, a proposed upper local attract point parameter to represent the solution between the current best solution and the current mean solution, is calculated by
(22)Pu=λ3×xbest+1−λ3×M.Pl, a proposed lower local attract point parameter to describe the solution between the current worst solution and the current mean solution, is calculated by
(23)Pl=λ4×xworst+1−λ4×M.The historical populationXoldXold=xold,1,xold,2,⋯,xold,N was first generated by
(24)Xold=X,ifPswitch≤0.5,Xold,otherwise.Then,Xold is processed by Xold=permutingXold, where permuting is a random shuffling function, which is to sort all individuals of Xold in random order [41].The local exploitation strategy of EJAYA can be described as(25)vi=xi+λ5×Pu−xi−λ6×Pl−xi,i=1,2,⋯,N.The global exploration strategy can be described as(26)vi=xi+κ×xold,i−xi,i=1,2,⋯,N,where λ3,λ4,λ5 and λ6 are all random numbers in [0, 1].


#### 3.2.2. Layup Design with Enhanced JAYA Algorithm

The objective functionDis is calculated by the difference between the optimal LPs and the actual LPs of the real laminate configuration.
(27)FindX=x1,x2,⋯,xneg,xi∈Di,TominimizeDis=∑i=13ViA−Vi,optiA2+ViD−Vi,optiD2,Di=45,−45,0,90.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].Pselect between 0 and 1 is used in EJAYA procedure to selected probability for local or global exploration strategy. If Pselect>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 Pselect≤0.5, the global exploration strategy is selected by Equation (24) and Equation (26).


## 3.1. First-Stage Optimization

### 3.1.1. Feasible Region of Lamination Parameters

Explicit expressions for solving feasible regions of LPs have been advanced progressively by many authors [25, 50–52]. LPs are restricted by the constraints as Equation (17) for symmetric laminate with layer angles limited to 0°, 90°, +45°, and -45°.
(15)2V1X−V2X−1≤0,2V3X+V2X−1≤0,X=A,D,4ViD+rViA+r≥ViA+r4,i=1,2,3,r=±1,162V1D+rV2D+rs2V1A+rV2A+rs≥2V1A+rV2A+rs4,r=±1,s=1,−3,162V3D+rV2D+rs2V3A+rV2A+rs≥2V3A+rV2A+rs4,r=±1,s=1,−3,4V1D+rV3D+sV1A+rV3A+s≥V1A+rV3A+s4,r=±1,s=±1.


### 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].
(16)FindX=h,V1A,V2A,V3A,V1D,V2D,V3D,,(17)TominimizeWX,SubjecttoVf≥VreqandEq.18.The constraints are all transformed as(18)dcjX≤0,j=1,2,⋯,nc,where dcjX is the jth design constraint and nc is the number of constraints. The penalty function method is used as follows [38]:
(19)fpenaltyX=1+ε1×vε2×WX;v=∑i=1ncmax0,dcjX,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 Section2.4. This solution is used within the optimization process aiming at minimizing the penalty function as Equation (19). The sensitivities of the penalty function can be obtained by the difference method.


## 3.1.1. Feasible Region of Lamination Parameters

Explicit expressions for solving feasible regions of LPs have been advanced progressively by many authors [25, 50–52]. LPs are restricted by the constraints as Equation (17) for symmetric laminate with layer angles limited to 0°, 90°, +45°, and -45°.
(15)2V1X−V2X−1≤0,2V3X+V2X−1≤0,X=A,D,4ViD+rViA+r≥ViA+r4,i=1,2,3,r=±1,162V1D+rV2D+rs2V1A+rV2A+rs≥2V1A+rV2A+rs4,r=±1,s=1,−3,162V3D+rV2D+rs2V3A+rV2A+rs≥2V3A+rV2A+rs4,r=±1,s=1,−3,4V1D+rV3D+sV1A+rV3A+s≥V1A+rV3A+s4,r=±1,s=±1.


## 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].
(16)FindX=h,V1A,V2A,V3A,V1D,V2D,V3D,,(17)TominimizeWX,SubjecttoVf≥VreqandEq.18.The constraints are all transformed as(18)dcjX≤0,j=1,2,⋯,nc,where dcjX is the jth design constraint and nc is the number of constraints. The penalty function method is used as follows [38]:
(19)fpenaltyX=1+ε1×vε2×WX;v=∑i=1ncmax0,dcjX,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 Section2.4. This solution is used within the optimization process aiming at minimizing the penalty function as Equation (19). The sensitivities of the penalty function can be obtained by the difference method.


## 3.2. 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]:
(20)vi=xi+λ1×x−bestxi−λ2×xworst−xi,i=1,2,⋯,N,where N is the population size, xi is the ith solution in the population, and xbest and xworst are the best and worst solutions in the population, respectively. Further, λ1 and λ2 are two random numbers in [0, 1]. Better solution of xi and vi will be selected into the population to accelerate convergence.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
(21)M=1N∑i=1Nxi,i=1,2,⋯,N.Pu, a proposed upper local attract point parameter to represent the solution between the current best solution and the current mean solution, is calculated by
(22)Pu=λ3×xbest+1−λ3×M.Pl, a proposed lower local attract point parameter to describe the solution between the current worst solution and the current mean solution, is calculated by
(23)Pl=λ4×xworst+1−λ4×M.The historical populationXoldXold=xold,1,xold,2,⋯,xold,N was first generated by
(24)Xold=X,ifPswitch≤0.5,Xold,otherwise.Then,Xold is processed by Xold=permutingXold, where permuting is a random shuffling function, which is to sort all individuals of Xold in random order [41].The local exploitation strategy of EJAYA can be described as(25)vi=xi+λ5×Pu−xi−λ6×Pl−xi,i=1,2,⋯,N.The global exploration strategy can be described as(26)vi=xi+κ×xold,i−xi,i=1,2,⋯,N,where λ3,λ4,λ5 and λ6 are all random numbers in [0, 1].


### 3.2.2. Layup Design with Enhanced JAYA Algorithm

The objective functionDis is calculated by the difference between the optimal LPs and the actual LPs of the real laminate configuration.
(27)FindX=x1,x2,⋯,xneg,xi∈Di,TominimizeDis=∑i=13ViA−Vi,optiA2+ViD−Vi,optiD2,Di=45,−45,0,90.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].Pselect between 0 and 1 is used in EJAYA procedure to selected probability for local or global exploration strategy. If Pselect>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 Pselect≤0.5, the global exploration strategy is selected by Equation (24) and Equation (26).


## 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]:
(20)vi=xi+λ1×x−bestxi−λ2×xworst−xi,i=1,2,⋯,N,where N is the population size, xi is the ith solution in the population, and xbest and xworst are the best and worst solutions in the population, respectively. Further, λ1 and λ2 are two random numbers in [0, 1]. Better solution of xi and vi will be selected into the population to accelerate convergence.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
(21)M=1N∑i=1Nxi,i=1,2,⋯,N.Pu, a proposed upper local attract point parameter to represent the solution between the current best solution and the current mean solution, is calculated by
(22)Pu=λ3×xbest+1−λ3×M.Pl, a proposed lower local attract point parameter to describe the solution between the current worst solution and the current mean solution, is calculated by
(23)Pl=λ4×xworst+1−λ4×M.The historical populationXoldXold=xold,1,xold,2,⋯,xold,N was first generated by
(24)Xold=X,ifPswitch≤0.5,Xold,otherwise.Then,Xold is processed by Xold=permutingXold, where permuting is a random shuffling function, which is to sort all individuals of Xold in random order [41].The local exploitation strategy of EJAYA can be described as(25)vi=xi+λ5×Pu−xi−λ6×Pl−xi,i=1,2,⋯,N.The global exploration strategy can be described as(26)vi=xi+κ×xold,i−xi,i=1,2,⋯,N,where λ3,λ4,λ5 and λ6 are all random numbers in [0, 1].


## 3.2.2. Layup Design with Enhanced JAYA Algorithm

The objective functionDis is calculated by the difference between the optimal LPs and the actual LPs of the real laminate configuration.
(27)FindX=x1,x2,⋯,xneg,xi∈Di,TominimizeDis=∑i=13ViA−Vi,optiA2+ViD−Vi,optiD2,Di=45,−45,0,90.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].Pselect between 0 and 1 is used in EJAYA procedure to selected probability for local or global exploration strategy. If Pselect>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 Pselect≤0.5, the global exploration strategy is selected by Equation (24) and Equation (26).


## 4. Results

### 4.1. Thermal Flutter Behavior with Lamination Parameter

In this section, thermal flutter behaviors of the laminated panel are discussed in the design space of LPs. Firstly, panel flutter behaviors without thermal effect are considered on LP lane. Secondly, flutter contours on LPs considering temperature variation are further discussed.


#### 4.1.1. Panel Flutter Behavior without Thermal Effect Consideration

When temperature variation is not considered, the flutter speed of panel structure mainly depends on theD matrix, that is, just relates to V1D,V2D,V3D,V4D. If the layer angles are limited to 0°, 90°, +45°, and -45°, then V4D=0. Therefore, the influence of changes of V1D,V2D,V3D on flutter speed is studied firstly. Figure 5 shows the contours of flutter speed on V1D−V2D plane with V3D=0. V3D=0 means that both D16 and D26 are 0, there is no bending-torsion coupling, and the laminate presents an orthotropic characteristic. It can be seen that when V2D is greater than -0.4 (region I), the change of flutter speed mainly depends on the variation of V1D. The maximum critical speed approximately appears at V1D=0. The flutter speed shows a decreasing trend with either increase or decrease of V1D. While when V2D is less than -0.4 (region II), the change of flutter speed mainly depends on variation of V2D. 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 5
V10−V20 plane with V1A,V2A,V3A,V4A,V3D,V4D=0,0,0,0,0,0.Figure6 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.Figure 6
Mode shapes of first six modes.
(a)
Mode 1(b)
Mode 2(c)
Mode 3(d)
Mode 4(e)
Mode 5(f)
Mode 6Figure 7
V-g and v-f curves for point A.
(a)
v-g curve(b)
v-f curveFigure 8
V-g and v-f curves for point B.
(a)
v-g curve(b)
v-f curvePoint AV1D,V2D,V3D=0,0,0.Point BV1D,V2D,V3D=0,−1,0.In order to investigate the influence ofV3D on flutter behavior, flutter speed change curves of point A and point B with V3D variation are given in Figure 9. For point A, the maximum critical speed appears almost at V3D=0; either increase or decrease of V3D will lead to flutter speed reduction, while the minimum critical speed appears almost at V3D=0 for point B.Figure 9
Flutter speed curves withV3D variation of points A and B.
(a)
Point A(b)
Point B


#### 4.1.2. CTEs with Lamination Parameters

Panels usually have the trend of in-plane expansion when subjected to thermal load. Due to boundary conditions, the panel will produce additional thermally induced in-plane compression load, 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. Influence ofV1A,V2A,V3A 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 V1A,V2A,V3A=0.5,0.4,0, and the minimum value of CTE_y is 8.73e-07 at V1A,V2A,V3A=−0.5,0.4,0. At V1A,V2A,V3A=0,0,0, CTE_x=CTE_y=2.43e−06. The maximum CTE_x and CTE_y of 1.56e-04 both appear at V1A,V2A,V3A=0,‐0.4,‐0.5. Figures 11 and 12 show the contours of CTEs on V1A‐V2A plane for V3A=0. Figure 13 shows the contours of CTEs on V2A‐V3A plane for V1A=0.Figure 10
Discrete points of CTEs.Figure 11
CTE_x contours ofV1AandV2A.Figure 12
CTE_y contours ofV1AandV2A.Figure 13
CTE contours onV2A‐V3A plane.


#### 4.1.3. 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.Point CV1A,V2A,V3A,V4A,V1D,V2D,V3D,V4D=0,0,0,0,0,0,0,0, CTE_x=CTE_y=2.43e−06.Point DV1A,V2A,V3A,V4A,V1D,V2D,V3D,V4D=0,0,0.3,0,0,0,0,0, CTE_x=CTE_y=6.23e−05.For point C, Figure14 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.Figure 14
Frequency change curve with temperature variation for point C.Figure 15
Frequency change curve with temperature variation for point D.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 Figure16(a). The flutter mode is bending-torsion coupling flutter dominated by second-order 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.Figure 16
Flutter change curve with temperature variation.
(a)
Point C(b)
Point D


### 4.2. 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 ish=0.0019032m, and V1A,V2A,V3A,V4A,V1D,V2D,V3D,V4D=0.005016,0.4444,0,0,−0.2132,−0.1780,0,0. 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).Figure 17
Frequencies and mode shapes of first six modes for optimal solution.
(a)
f1=99.02HZ(b)
f2=253.98HZ(c)
f3=307.45HZ(d)
f4=430.08HZ(e)
f5=584.52HZ(f)
f6=618.79HZFigure 18
V-g and v-f curves for optimal solution.
(a)
v-g curve(b)
v-f curve


### 4.3. 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 theDis 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.Table 1
Optimal results.
LPsV1AV2AV3AV4AV1DV2DV3DV1ADisTarget LPs0.005016460.4444311800-0.21316306-0.1799920200Optimum 10.10.4-0.10-0.209-0.158-0.03700.023Layup 1[45/-45/902/-45/02/90/02]sOptimum 200.400-0.2175-0.188-0.0060.01040.0022Layup 2[-45/45/75/-75/15/90/-15/90/02]SThe optimized laminate configuration is checked for frequency and thermal flutter performance, shown in Table2. 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.Table 2
Comparisons of frequency and flutter speed.
Mode 1Mode 2Mode 3Mode 4Mode 5Mode 6Flutter speedFlutter modesTarget LPs105.37267.85324.4452.77615.74651.23271.7Modes 3 and 2Layup 1114.48272.87334.05459.27624.8660.893222.2Modes 3 and 2Error8.65%1.87%2.97%1.44%1.47%1.49%Layup 2113.4274.37332.72458.13625.41656.313250.8Modes 3 and 2Error7.62%2.43%2.56%1.18%1.57%0.78%


## 4.1. Thermal Flutter Behavior with Lamination Parameter

In this section, thermal flutter behaviors of the laminated panel are discussed in the design space of LPs. Firstly, panel flutter behaviors without thermal effect are considered on LP lane. Secondly, flutter contours on LPs considering temperature variation are further discussed.


### 4.1.1. Panel Flutter Behavior without Thermal Effect Consideration

When temperature variation is not considered, the flutter speed of panel structure mainly depends on theD matrix, that is, just relates to V1D,V2D,V3D,V4D. If the layer angles are limited to 0°, 90°, +45°, and -45°, then V4D=0. Therefore, the influence of changes of V1D,V2D,V3D on flutter speed is studied firstly. Figure 5 shows the contours of flutter speed on V1D−V2D plane with V3D=0. V3D=0 means that both D16 and D26 are 0, there is no bending-torsion coupling, and the laminate presents an orthotropic characteristic. It can be seen that when V2D is greater than -0.4 (region I), the change of flutter speed mainly depends on the variation of V1D. The maximum critical speed approximately appears at V1D=0. The flutter speed shows a decreasing trend with either increase or decrease of V1D. While when V2D is less than -0.4 (region II), the change of flutter speed mainly depends on variation of V2D. 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 5
V10−V20 plane with V1A,V2A,V3A,V4A,V3D,V4D=0,0,0,0,0,0.Figure6 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.Figure 6
Mode shapes of first six modes.
(a)
Mode 1(b)
Mode 2(c)
Mode 3(d)
Mode 4(e)
Mode 5(f)
Mode 6Figure 7
V-g and v-f curves for point A.
(a)
v-g curve(b)
v-f curveFigure 8
V-g and v-f curves for point B.
(a)
v-g curve(b)
v-f curvePoint AV1D,V2D,V3D=0,0,0.Point BV1D,V2D,V3D=0,−1,0.In order to investigate the influence ofV3D on flutter behavior, flutter speed change curves of point A and point B with V3D variation are given in Figure 9. For point A, the maximum critical speed appears almost at V3D=0; either increase or decrease of V3D will lead to flutter speed reduction, while the minimum critical speed appears almost at V3D=0 for point B.Figure 9
Flutter speed curves withV3D variation of points A and B.
(a)
Point A(b)
Point B


### 4.1.2. CTEs with Lamination Parameters

Panels usually have the trend of in-plane expansion when subjected to thermal load. Due to boundary conditions, the panel will produce additional thermally induced in-plane compression load, 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. Influence ofV1A,V2A,V3A 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 V1A,V2A,V3A=0.5,0.4,0, and the minimum value of CTE_y is 8.73e-07 at V1A,V2A,V3A=−0.5,0.4,0. At V1A,V2A,V3A=0,0,0, CTE_x=CTE_y=2.43e−06. The maximum CTE_x and CTE_y of 1.56e-04 both appear at V1A,V2A,V3A=0,‐0.4,‐0.5. Figures 11 and 12 show the contours of CTEs on V1A‐V2A plane for V3A=0. Figure 13 shows the contours of CTEs on V2A‐V3A plane for V1A=0.Figure 10
Discrete points of CTEs.Figure 11
CTE_x contours ofV1AandV2A.Figure 12
CTE_y contours ofV1AandV2A.Figure 13
CTE contours onV2A‐V3A plane.


### 4.1.3. 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.Point CV1A,V2A,V3A,V4A,V1D,V2D,V3D,V4D=0,0,0,0,0,0,0,0, CTE_x=CTE_y=2.43e−06.Point DV1A,V2A,V3A,V4A,V1D,V2D,V3D,V4D=0,0,0.3,0,0,0,0,0, CTE_x=CTE_y=6.23e−05.For point C, Figure14 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.Figure 14
Frequency change curve with temperature variation for point C.Figure 15
Frequency change curve with temperature variation for point D.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 Figure16(a). The flutter mode is bending-torsion coupling flutter dominated by second-order 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.Figure 16
Flutter change curve with temperature variation.
(a)
Point C(b)
Point D


## 4.1.1. Panel Flutter Behavior without Thermal Effect Consideration

When temperature variation is not considered, the flutter speed of panel structure mainly depends on theD matrix, that is, just relates to V1D,V2D,V3D,V4D. If the layer angles are limited to 0°, 90°, +45°, and -45°, then V4D=0. Therefore, the influence of changes of V1D,V2D,V3D on flutter speed is studied firstly. Figure 5 shows the contours of flutter speed on V1D−V2D plane with V3D=0. V3D=0 means that both D16 and D26 are 0, there is no bending-torsion coupling, and the laminate presents an orthotropic characteristic. It can be seen that when V2D is greater than -0.4 (region I), the change of flutter speed mainly depends on the variation of V1D. The maximum critical speed approximately appears at V1D=0. The flutter speed shows a decreasing trend with either increase or decrease of V1D. While when V2D is less than -0.4 (region II), the change of flutter speed mainly depends on variation of V2D. 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 5
V10−V20 plane with V1A,V2A,V3A,V4A,V3D,V4D=0,0,0,0,0,0.Figure6 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.Figure 6
Mode shapes of first six modes.
(a)
Mode 1(b)
Mode 2(c)
Mode 3(d)
Mode 4(e)
Mode 5(f)
Mode 6Figure 7
V-g and v-f curves for point A.
(a)
v-g curve(b)
v-f curveFigure 8
V-g and v-f curves for point B.
(a)
v-g curve(b)
v-f curvePoint AV1D,V2D,V3D=0,0,0.Point BV1D,V2D,V3D=0,−1,0.In order to investigate the influence ofV3D on flutter behavior, flutter speed change curves of point A and point B with V3D variation are given in Figure 9. For point A, the maximum critical speed appears almost at V3D=0; either increase or decrease of V3D will lead to flutter speed reduction, while the minimum critical speed appears almost at V3D=0 for point B.Figure 9
Flutter speed curves withV3D variation of points A and B.
(a)
Point A(b)
Point B


## 4.1.2. CTEs with Lamination Parameters

Panels usually have the trend of in-plane expansion when subjected to thermal load. Due to boundary conditions, the panel will produce additional thermally induced in-plane compression load, 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. Influence ofV1A,V2A,V3A 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 V1A,V2A,V3A=0.5,0.4,0, and the minimum value of CTE_y is 8.73e-07 at V1A,V2A,V3A=−0.5,0.4,0. At V1A,V2A,V3A=0,0,0, CTE_x=CTE_y=2.43e−06. The maximum CTE_x and CTE_y of 1.56e-04 both appear at V1A,V2A,V3A=0,‐0.4,‐0.5. Figures 11 and 12 show the contours of CTEs on V1A‐V2A plane for V3A=0. Figure 13 shows the contours of CTEs on V2A‐V3A plane for V1A=0.Figure 10
Discrete points of CTEs.Figure 11
CTE_x contours ofV1AandV2A.Figure 12
CTE_y contours ofV1AandV2A.Figure 13
CTE contours onV2A‐V3A plane.


## 4.1.3. 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.Point CV1A,V2A,V3A,V4A,V1D,V2D,V3D,V4D=0,0,0,0,0,0,0,0, CTE_x=CTE_y=2.43e−06.Point DV1A,V2A,V3A,V4A,V1D,V2D,V3D,V4D=0,0,0.3,0,0,0,0,0, CTE_x=CTE_y=6.23e−05.For point C, Figure14 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.Figure 14
Frequency change curve with temperature variation for point C.Figure 15
Frequency change curve with temperature variation for point D.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 Figure16(a). The flutter mode is bending-torsion coupling flutter dominated by second-order 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.Figure 16
Flutter change curve with temperature variation.
(a)
Point C(b)
Point D


## 4.2. 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 ish=0.0019032m, and V1A,V2A,V3A,V4A,V1D,V2D,V3D,V4D=0.005016,0.4444,0,0,−0.2132,−0.1780,0,0. 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).Figure 17
Frequencies and mode shapes of first six modes for optimal solution.
(a)
f1=99.02HZ(b)
f2=253.98HZ(c)
f3=307.45HZ(d)
f4=430.08HZ(e)
f5=584.52HZ(f)
f6=618.79HZFigure 18
V-g and v-f curves for optimal solution.
(a)
v-g curve(b)
v-f curve


## 4.3. 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 theDis 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.Table 1
Optimal results.
LPsV1AV2AV3AV4AV1DV2DV3DV1ADisTarget LPs0.005016460.4444311800-0.21316306-0.1799920200Optimum 10.10.4-0.10-0.209-0.158-0.03700.023Layup 1[45/-45/902/-45/02/90/02]sOptimum 200.400-0.2175-0.188-0.0060.01040.0022Layup 2[-45/45/75/-75/15/90/-15/90/02]SThe optimized laminate configuration is checked for frequency and thermal flutter performance, shown in Table2. 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.Table 2
Comparisons of frequency and flutter speed.
Mode 1Mode 2Mode 3Mode 4Mode 5Mode 6Flutter speedFlutter modesTarget LPs105.37267.85324.4452.77615.74651.23271.7Modes 3 and 2Layup 1114.48272.87334.05459.27624.8660.893222.2Modes 3 and 2Error8.65%1.87%2.97%1.44%1.47%1.49%Layup 2113.4274.37332.72458.13625.41656.313250.8Modes 3 and 2Error7.62%2.43%2.56%1.18%1.57%0.78%


## 5. 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.



---

*Source: 5201514-2022-03-12.xml*