# Efficient Method for Aeroelastic Tailoring of Composite Wing to Minimize Gust Response

Authors: Yang Yu; Zhengjie Wang; Shijun Guo
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

Aeroelastic tailoring of laminated composite structure demands relatively high computational time especially for dynamic problem. This paper presents an efficient method for aeroelastic dynamic response analysis with significantly reduced computational time. In this method, a relationship is established between the maximum aeroelastic response and quasi-steady deflection of a wing subject to a dynamic loading. Based on this relationship, the time consuming dynamic response can be approximated by a quasi-steady deflection analysis in a large proportion of the optimization process. This method has been applied to the aeroelastic tailoring of a composite wing of a tailless aircraft for minimum gust response. The results have shown that 20%–36% gust response reduction has been achieved for this case. The computational time of the optimization process has been reduced by 90% at the cost of accuracy reduction of 2~4% comparing with the traditional dynamic response analysis.


---

## Body

## 1. Introduction

In aircraft design, the use of composite materials and increase of wing aspect ratio for higher aerodynamic efficiency becomes a trend [1]. As a result, the wing becomes more flexible and prone to not only static aeroelastic deformation and aileron reversal, but also aeroelastic stability [2]. As one of the approaches to solve the problem, aeroelastic tailoring techniques has been applied to the optimal design of composite wing structures [3, 4]. In addition, gust response is a critical concern and has attracted research attention especially for lightweight aircraft of high aspect ratio wing. In this particular subject, substantial research activities have been carried out to develop various active and passive control technologies for gust alleviation of aircraft since 1990s [5–10]. Therefore aeroelasticity should be considered in the early phase of aircraft design process [11]. For this purpose, an efficient method to predict the gust response of an aircraft with the aeroelastic coupling effect of the wing is required.This current investigation is to develop an efficient method for gust response analysis applicable to optimization of a composite wing structure for minimum gust response. In aircraft design, there are two types of gust model and methods to predict the response according to the airworthiness regulation EASA CS-25 [12]. The discrete gust model is the simple one; while the continuous gust based on Von Kaman’s gust power spectrum density model is usually solved by using random process theory. To focus our study on the efficient dynamic response analysis, the discrete gust model is chosen in the current investigation.When the dynamic aeroelastic coupling effect of the wing is considered, the calculation for maximum gust response performed in either frequency-domain or time-domain is time consuming. For aeroelastic tailoring, this becomes more concerned since such calculation is usually conducted in the gradient calculation of each design variable and each iterative cycle in a traditional optimization process. Although efforts have been made to simplify the structural models, the time consuming analysis of gust response remains as a challenge in this field [13, 14]. In a previous study of aircraft structure optimization subject to gust response, the efficiency was improved in two ways [15], firstly by saving the modal analysis time for reduced structural model at each variation step, and secondly calculating the gust response only in limited iteration steps to check if the optimization was kept in track like the other work [13]. In the FFAST (“Future Fast Aeroelastic Simulation Technologies”) project, modeling and optimization methods were studied to predict the worst-case gust load in fast and efficient manner [16]. Their work also considered changes made to aircraft structure as part of the evolving design process [17].However, when more detailed structure optimization such as aeroelastic tailoring is required, dynamic response result is required in each iteration step or more frequently. Therefore an efficient method for calculating the gust response rather than cutting off the calculation is required. This paper presents a method for establishing a relationship of the maximum gust response with quasi-steady deflection of a wing. Based on this relationship, the time consuming dynamic response analysis can be replaced by a quasi-steady deflection analysis in a large proportion of the iteration cycles during the optimization process.In the optimization of a wing structure subject to gust, a response peak value at a representative point is normally taken as a measurement of the gust response. The center of gravity of a rigid aircraft or the wingtip for a flexible wing is usually taken as the measurement point. This implies that it is unnecessary to perform the time consuming traditional dynamic response analysis in time history or frequency domain. However, the existing approach to evaluate the response peak value is to perform the traditional dynamic response analysis unless a more efficient alternative method is available. This current study presents a quasi-static (Q-S) method to predict the gust response peak value based on a mixture of static and dynamic response analysis in time-domain. This method increases the efficiency of gust response analysis significantly especially when applied to structure optimization where many iterations of the response calculation are performed. The investigation was conducted in three stages. In the first stage, an expression was established for the energy conversion between the work done by external dynamic force and structural deformation at the time when the response reaches its first maximum value. A proportional relationship was created between the elastic energy of a structure and work done by external force and resulted in a Q-S ratio. In the second stage, the relationship was extended to the optimization process in which the structural stiffness has a small variation in the iteration cycles. Based on the relationship, the dynamic response peak value of a structure at an iteration cycle can be predicted based on the Q-S ratio and the response obtained in previous cycle. In the third stage, the time period to reach the first response peak value was determined for the structure having a small stiffness variation during the optimization. By applying this method, the computational time for a structural optimization subject to dynamic response can be reduced significantly. To demonstrate the method, a small aircraft has been taken as example to minimize its gust response by tailoring the composite wing laminate layup. Different cases are studied including gust frequency and wing sweepback angles. The results show that this method is able to predict the gust response accurately with the computational time over the optimization reduced by over 85% at the cost of accuracy reduction by less than 5% comparing with the traditional method of traditional dynamic response analysis.


## 2. Theoretical Study of the Method

### 2.1. The Quasi-Static Model for Dynamic Response Analysis

The work done on an elastic body by external forces is absorbed and transformed to kinetic and elastic potential energy in a form of movement and deformation. For a linear elastic body subjected to a dynamic load, the total energyE includes both elastic potential and kinetic energy EK written as [18](1)E=12XTKX+EK,where X is the elastic deformation vector and K is the stiffness matrix of a structure.An estimation of the stiffness of a structure,KE, can be obtained by the force applied against the deflection measured at representative points. In this article, an estimation of stiffness KE of the wing structure is calculated from wingtip displacement and the total bending moment and torque on the wing obtained from static aeroelastic analysis. The estimated equivalent KE is used in the following Q-Static method for gust response analysis.AssumingTP is time when the dynamic response reaches the first peak value as shown in Figure 1. The dynamic response such as the wingtip displacement reaches a peak value XE at TP. At this moment, the velocity and kinetic energy EK of the structure becomes zero [19]. The system is in Q-Static status and total energy E can be simplified to an expression of elastic potential energy EEP.(2)E=EEP=12KEXE2.Figure 1
An illustration ofTP, XE, and deviations between two optimization circles.The work done by the external dynamic forceQ(t) through the structural deformation over a time period until a specified time such as the peak value of dynamic response can be written as(3)WP=∫0TPQtHtdt,where H(t) is the corresponding velocity of the structural response to Q(t).It is noted that theWP gives an approximate total work done based on the estimated EEP. For a small variation of structural KE at each optimization cycle, it is a relative value against the previous one and will be updated using a “calibration process” as described later in Section 2.5 during the optimization.In the current optimization study, the response variation between two contiguous iteration cycles is caused by a small stiffness variation by tailoring the laminate layup of a wing structure. The equivalent structure stiffness at iteration cycleN can be represented by the stiffness at previous cycle N-1:(4)KE∣N=ηKE∣N-1,in whichη is the ratio of the stiffness variation between two contiguous iteration cycles. For linear elastic structural deformation, the response peaks at the two iteration steps are also in the same proportion as the stiffness:(5)XE∣N=1ηXE∣N-1.This leads to the relationship of theEEP between two iteration cycles:(6)EEP∣N=12KE∣NXE∣N2=1ηEEP∣N-1.As structural stiffness is varied, the time to reach the response peak value is also changed from the previous iteration cycle:(7)TP∣N=TP∣N-1+ΔTη,in whichΔTη represents the time difference of peak gust response between cycle N and N-1. The work done by the external force between two iteration cycles is related by(8)WEP∣N=∫0TP∣NQtHE∣Ntdt=AN∫0TP∣N-1QtH-E∣Ntdt+∫TP∣N-1TP∣N-1+ΔTQtH-E∣Ntdt=AN∫0TP∣N-1QtH-E∣N-1t+H-E∣Nt-H-E∣N-1tdt+∫TP∣N-1TP∣N-1+ΔTQtH-E∣Ntdt=ANAN-1AN-1∫0TP∣N-1QtH-E∣N-1tdt+AN-1∫0TP∣N-1QtH-E∣Nt-H-E∣N-1tdt+∫TP∣N-1TP∣N-1+ΔTQtH-E∣Ntdt=ANAN-1WEP∣N-1+ΔW,in whichHt=AH-t, and AN/AN-1 is the Q-S ratio of response amplitude in two contiguous iteration cycles given by(9)ANAN-1=XE∣NXE∣N-1=1η.From (6) and (8), the following relation between EEP and WEP is obtained in two contiguous iteration cycles:(10)WEP∣NEEP∣N=η-1WEP∣N-1+ΔWη-1EEP∣N-1=WEP∣N-1EEP∣N-1+ΔWEEP∣N-1,where ΔW/EEP∣N-1 represents a deviation of the ratio WEP/EEP between two contiguous iteration cycles.In (10), ΔW is defined in (8) as(11)ΔW=AN-1∫0TP∣N-1QtH-E∣Nt-H-E∣N-1tdt+∫TP∣N-1TP∣N-1+ΔTQtH-E∣Ntdt.TheΔW is expressed in two terms: the first term is due to the variation of deformation; the second term is due to the variation of peak response time all caused by tailoring the structure such as the laminate layup. For a small ΔW/EEP∣N-1 below a value, it can be negligible as discussed in Section 2.5. In this case, a simplified ratio can be obtained from (10):(12)EEP∣NEEP∣N-1=WEP∣NWEP∣N-1.When a ratio is obtained at iteration stepN-1, an approximate XE at iteration cycle N can be obtained based on (12) instead of performing the traditional dynamic response analysis again. The more iteration cycles it performs, the more computational time to save.In this current study, the (1-cos) discrete gust model as defined in airworthiness regulation CS-25 is considered:(13)Ug=Vds21-cos⁡ωgt,where Vds is a specified design gust velocity and ωg=πs/tH is the gust frequency, in which s=Vt is the penetrating distance of an aircraft at an equivalent flight speed V into a discrete gust at time t. H is gust gradient distance.The symmetric aerodynamic force due to the gust can be calculated by(14)Qt=12ρV2Sa1FUgV=qVds21-cos⁡ωgt,where ρ is the air density, S is the reference wing area, and a1 is the lift curve slope of the wing. F is the gust alleviating factor expressed by(15)F=0.88μ5.3+μ,μ=2mρSc-a1.It is noted that the aboveF is defined for a typical gust gradient distance H=12.5c where c is the mean chord. However the Q-static method is not related to the factor. It is only the accuracy of the gust response that will be affected when the factor is adapted for other H cases. Corresponding to the (1-cos) gust, the responding wingtip deflection can also be expressed as a (1-cos) function. Hence the response velocity HE(t) before approaching the peak time TP can be expressed in a sinusoidal function by differentiating the wingtip deflection:(16)HEt=AπTP-1sin⁡πTP-1t,where A represents the wingtip deflection amplitude at the response timeTP. Based on (3), (14), and (16), the WEP can be rewritten as(17)WEP=∫0TPQt×HEtdt=∫0TPq×Vde2×1-cos⁡ωgt×AπTP-1sin⁡πTP-1tdt=A∫0TPQt×H-Etdt.According to (2) and (17), (12) can be rewritten as(18)KE∣NXE∣N2KE∣N-1XE∣N-12=∫0TP∣Nq×Vde/2V1-cos⁡ωgt×ANπTP∣N-1sin⁡πTP∣N-1tdt∫0TP∣N-1q×Vde/2V1-cos⁡ωgt×AN-1πTP∣N-1-1sin⁡πTP∣N-1-1tdt.By extracting the parametersq, π, TP, AN, and AN-1 out of the time integration and substituting (9) into (18), a simple recursive relation of XE∣N at TP∣N can be obtained from the XE∣N-1 as shown:(19)XE∣N=TP∣N-1∫0TP∣N1-cos⁡ωgt×sin⁡πTP∣N-1tdtTP∣N-1-1∫0TP∣N-11-cos⁡ωgt×sin⁡πTP∣N-1-1tdt×KE∣N-1KE∣NXE∣N-1.In (19), since the TP∣N in the iteration cycle N is different from previous iteration cycle N-1, it will be predicted in each iteration cycle. The details of how to predict the TP∣N are presented in the following section.Although the deviationΔW/EEP|N-1 as shown in (10) may be negligible, it will accumulate when (19) is used to predict the XE in subsequent cycles beyond cycleN. After number of cycles, it is necessary to update the approximated ratio as shown in (12) by performing a full dynamic response analysis without approximation. This so-called “calibration process” is described in Section 2.5.


### 2.2. Prediction of the Peak Response Time

As an earlier description for (2), the wing structure was simplified to a second-order system with an equivalent stiffness. The system response time TP to reach its peak value as illustrated in Figure 1 can be expressed in (19), in which ωn is the structural natural frequency and ζ the damping ratio:(20)TP=πωn1-ζ2.An expression of the angular value at peak response timeTP in general form can be written as(21)TPωn=fωg,ωn,ζ,where ωg is the gust frequency. Since the constant damping ζ has little influence to the frequency of the system having a small variation of stiffness, the response frequency in the iteration cycle N has a small difference from that in cycle N-1 due to a small stiffness variation:(22)ωN=ωN-1+Δω.Similarly the angular value of response functionTP∣NωN has a small difference from TP∣N-1ωN-1:(23)TP∣NωN=TP∣N-1ωN-1+eΔω.In this study, the first bending mode of the wing is taken as an approximateωN since it is in the same frequency range and hence dominates the response to the gust. This approximation may produce a small deviation eΔω for the TP∣N prediction. Equation (23) can be further simplified as below when this deviation is initially ignored for small Δω:(24)TP∣NωN=TP∣N-1ωN-1.During optimization process, the angleTP∣N-1ωN-1 in cycle N-1 is obtained from the traditional dynamic response analysis. For the cycle N and following circles, the time TP at the peak gust response can be calculated based on (24). If TP∣NωNR represents the correct value, the deviation of the above predicted result TP∣NωNP due to the eΔω can be evaluated by(25)eTP∣NωNP=1-TP∣NωNPTP∣NωNR.Since theeΔω may increase with the variation of structure stiffness during the optimization, it is also necessary to obtain the correct gust response by performing a complete gust response analysis. The correct value will be used for calibrating the TP∣N+nωN+n to ensure the optimization converge. This is part of the “calibration process” similar to the TP∣N+n prediction process in the Q-Static method.IfXE∣P and XE∣R represent the predicted and the correct gust response values, respectively, the deviation of the predicted result due to all the approximations can be evaluated by(26)eP=1-XE∣PXE∣R.


### 2.3. Gust Response Analysis

The equations used in gust response analysis can be expressed in a simple form as [20](27)MX¨+DX˙+KX=FsteadyX+FunsteadyX,X˙,X¨X0=XN-1X˙0=X˙N-1,in whichM is mass matrix and D is structural damping matrix. Unsteady aerodynamic forces are also included in simulation, which is calculated based on Theodorsen function. The structural analysis is based on thin-walled wing box beam model with the stiffness of each section depending upon the fiber orientations of the skin and spar web laminate. A program “AeroBeamSaw” for static aeroelastic analysis and model analysis is linked with the gust response analysis performed in MATLAB and an optimization program performed in FORTRAN environment [12, 21].


### 2.4. Optimization Method and Objective Function

A gradient-based deterministic optimization method was employed to tailor a composite wing structure for minimum gust loads with the objective function expressed below:(28)MinimizeFdα=LSαLSα0η×XEαXEα0ξ2,in whichLS represents total shear force at wing root, XE represents the peak value of wingtip bending deflection, and αo and α represent a set of initial and tailored design variables respectively, that is, the fiber orientations in a range of 0~±90° of each layer of the laminated composite wing structure. In (28), η and ξ are weight coefficients, η is set at 5, and ξ is set at 0.15.


### 2.5. Calibration Process in the Optimization

In optimization procedure, there are two major parts, which are shown in Figure2. The first major part is Q-S method, which is used to predict wingtip deflection in this article. The second major part is calibration process, which is used to eliminate the error accumulation due to the simplification of the Q-S method.Figure 2
The flowchart of the Q-S method in structure optimization procedure.


#### 2.5.1. Q-S Method

Q-S method developed in the current investigation is based on a relationship between the dynamic response of a structure and the work done by an external force. For structural optimization, the relationship is extended to the responding deflections when the structure stiffness has a small variation between two iteration cycles as expressed in (12) and (24).At the starting point of the optimization when the iteration cycleN-1(N-1=0), a dynamic response analysis has to be performed for the initial structure design. Based on the results obtained in the cycle N-1, the deflection of the structure of a small stiffness variation in the iteration cycle N and subsequent cycles N+n(n=1,2,3,…) can be predicted.In this process, instead of performing a dynamic response analysis like the starting point, Q-S method should be used in each iteration cycle repeatedly by recalculating the parametersKE∣N, TP∣N and the Q-S deflection XE∣N-1. It means that the results from one dynamic response analysis can be used for number of iteration cycles with the parameters in (12) and (24) kept updating.


#### 2.5.2. Calibration Process

The calibration process in current study is designed on the purpose of improving prediction accuracy. There will be error accumulated in iteration cycles. The main cause of error comes from the variation of structure stiffness, which will affect (12) and (24), as error terms that have connection with small changes made to stiffness are omitted in order to get a simplified recursive relation. When calibration is launched, a traditional dynamic response analysis will be performed to update parameters used in recursive equations for the subsequent optimization process.The launch condition of calibration process is vital for the performance of Q-S method in optimization. There are two major problems.The first problem is that, during optimization procedure, with only a few iteration cycles having “real” data, there is no reliable direct observation method for the assessment of prediction accuracy. Indirect measurement will be taken to set the launch condition. In the current study, equivalent stiffnessKE∣N obtained from static aeroelastic analysis is used as the indirect measurement to trigger calibration process.The second problem is how much indirect measurement varies that can make it necessary to do the calibration. If the conditions are set too strict, it may lower down the efficiency of Q-S method. Otherwise, it may deteriorate Q-S method’s accuracy.In the current study, 4% stiffness variation is set as the criteria for performing the calibration process during optimization.


## 2.1. The Quasi-Static Model for Dynamic Response Analysis

The work done on an elastic body by external forces is absorbed and transformed to kinetic and elastic potential energy in a form of movement and deformation. For a linear elastic body subjected to a dynamic load, the total energyE includes both elastic potential and kinetic energy EK written as [18](1)E=12XTKX+EK,where X is the elastic deformation vector and K is the stiffness matrix of a structure.An estimation of the stiffness of a structure,KE, can be obtained by the force applied against the deflection measured at representative points. In this article, an estimation of stiffness KE of the wing structure is calculated from wingtip displacement and the total bending moment and torque on the wing obtained from static aeroelastic analysis. The estimated equivalent KE is used in the following Q-Static method for gust response analysis.AssumingTP is time when the dynamic response reaches the first peak value as shown in Figure 1. The dynamic response such as the wingtip displacement reaches a peak value XE at TP. At this moment, the velocity and kinetic energy EK of the structure becomes zero [19]. The system is in Q-Static status and total energy E can be simplified to an expression of elastic potential energy EEP.(2)E=EEP=12KEXE2.Figure 1
An illustration ofTP, XE, and deviations between two optimization circles.The work done by the external dynamic forceQ(t) through the structural deformation over a time period until a specified time such as the peak value of dynamic response can be written as(3)WP=∫0TPQtHtdt,where H(t) is the corresponding velocity of the structural response to Q(t).It is noted that theWP gives an approximate total work done based on the estimated EEP. For a small variation of structural KE at each optimization cycle, it is a relative value against the previous one and will be updated using a “calibration process” as described later in Section 2.5 during the optimization.In the current optimization study, the response variation between two contiguous iteration cycles is caused by a small stiffness variation by tailoring the laminate layup of a wing structure. The equivalent structure stiffness at iteration cycleN can be represented by the stiffness at previous cycle N-1:(4)KE∣N=ηKE∣N-1,in whichη is the ratio of the stiffness variation between two contiguous iteration cycles. For linear elastic structural deformation, the response peaks at the two iteration steps are also in the same proportion as the stiffness:(5)XE∣N=1ηXE∣N-1.This leads to the relationship of theEEP between two iteration cycles:(6)EEP∣N=12KE∣NXE∣N2=1ηEEP∣N-1.As structural stiffness is varied, the time to reach the response peak value is also changed from the previous iteration cycle:(7)TP∣N=TP∣N-1+ΔTη,in whichΔTη represents the time difference of peak gust response between cycle N and N-1. The work done by the external force between two iteration cycles is related by(8)WEP∣N=∫0TP∣NQtHE∣Ntdt=AN∫0TP∣N-1QtH-E∣Ntdt+∫TP∣N-1TP∣N-1+ΔTQtH-E∣Ntdt=AN∫0TP∣N-1QtH-E∣N-1t+H-E∣Nt-H-E∣N-1tdt+∫TP∣N-1TP∣N-1+ΔTQtH-E∣Ntdt=ANAN-1AN-1∫0TP∣N-1QtH-E∣N-1tdt+AN-1∫0TP∣N-1QtH-E∣Nt-H-E∣N-1tdt+∫TP∣N-1TP∣N-1+ΔTQtH-E∣Ntdt=ANAN-1WEP∣N-1+ΔW,in whichHt=AH-t, and AN/AN-1 is the Q-S ratio of response amplitude in two contiguous iteration cycles given by(9)ANAN-1=XE∣NXE∣N-1=1η.From (6) and (8), the following relation between EEP and WEP is obtained in two contiguous iteration cycles:(10)WEP∣NEEP∣N=η-1WEP∣N-1+ΔWη-1EEP∣N-1=WEP∣N-1EEP∣N-1+ΔWEEP∣N-1,where ΔW/EEP∣N-1 represents a deviation of the ratio WEP/EEP between two contiguous iteration cycles.In (10), ΔW is defined in (8) as(11)ΔW=AN-1∫0TP∣N-1QtH-E∣Nt-H-E∣N-1tdt+∫TP∣N-1TP∣N-1+ΔTQtH-E∣Ntdt.TheΔW is expressed in two terms: the first term is due to the variation of deformation; the second term is due to the variation of peak response time all caused by tailoring the structure such as the laminate layup. For a small ΔW/EEP∣N-1 below a value, it can be negligible as discussed in Section 2.5. In this case, a simplified ratio can be obtained from (10):(12)EEP∣NEEP∣N-1=WEP∣NWEP∣N-1.When a ratio is obtained at iteration stepN-1, an approximate XE at iteration cycle N can be obtained based on (12) instead of performing the traditional dynamic response analysis again. The more iteration cycles it performs, the more computational time to save.In this current study, the (1-cos) discrete gust model as defined in airworthiness regulation CS-25 is considered:(13)Ug=Vds21-cos⁡ωgt,where Vds is a specified design gust velocity and ωg=πs/tH is the gust frequency, in which s=Vt is the penetrating distance of an aircraft at an equivalent flight speed V into a discrete gust at time t. H is gust gradient distance.The symmetric aerodynamic force due to the gust can be calculated by(14)Qt=12ρV2Sa1FUgV=qVds21-cos⁡ωgt,where ρ is the air density, S is the reference wing area, and a1 is the lift curve slope of the wing. F is the gust alleviating factor expressed by(15)F=0.88μ5.3+μ,μ=2mρSc-a1.It is noted that the aboveF is defined for a typical gust gradient distance H=12.5c where c is the mean chord. However the Q-static method is not related to the factor. It is only the accuracy of the gust response that will be affected when the factor is adapted for other H cases. Corresponding to the (1-cos) gust, the responding wingtip deflection can also be expressed as a (1-cos) function. Hence the response velocity HE(t) before approaching the peak time TP can be expressed in a sinusoidal function by differentiating the wingtip deflection:(16)HEt=AπTP-1sin⁡πTP-1t,where A represents the wingtip deflection amplitude at the response timeTP. Based on (3), (14), and (16), the WEP can be rewritten as(17)WEP=∫0TPQt×HEtdt=∫0TPq×Vde2×1-cos⁡ωgt×AπTP-1sin⁡πTP-1tdt=A∫0TPQt×H-Etdt.According to (2) and (17), (12) can be rewritten as(18)KE∣NXE∣N2KE∣N-1XE∣N-12=∫0TP∣Nq×Vde/2V1-cos⁡ωgt×ANπTP∣N-1sin⁡πTP∣N-1tdt∫0TP∣N-1q×Vde/2V1-cos⁡ωgt×AN-1πTP∣N-1-1sin⁡πTP∣N-1-1tdt.By extracting the parametersq, π, TP, AN, and AN-1 out of the time integration and substituting (9) into (18), a simple recursive relation of XE∣N at TP∣N can be obtained from the XE∣N-1 as shown:(19)XE∣N=TP∣N-1∫0TP∣N1-cos⁡ωgt×sin⁡πTP∣N-1tdtTP∣N-1-1∫0TP∣N-11-cos⁡ωgt×sin⁡πTP∣N-1-1tdt×KE∣N-1KE∣NXE∣N-1.In (19), since the TP∣N in the iteration cycle N is different from previous iteration cycle N-1, it will be predicted in each iteration cycle. The details of how to predict the TP∣N are presented in the following section.Although the deviationΔW/EEP|N-1 as shown in (10) may be negligible, it will accumulate when (19) is used to predict the XE in subsequent cycles beyond cycleN. After number of cycles, it is necessary to update the approximated ratio as shown in (12) by performing a full dynamic response analysis without approximation. This so-called “calibration process” is described in Section 2.5.


## 2.2. Prediction of the Peak Response Time

As an earlier description for (2), the wing structure was simplified to a second-order system with an equivalent stiffness. The system response time TP to reach its peak value as illustrated in Figure 1 can be expressed in (19), in which ωn is the structural natural frequency and ζ the damping ratio:(20)TP=πωn1-ζ2.An expression of the angular value at peak response timeTP in general form can be written as(21)TPωn=fωg,ωn,ζ,where ωg is the gust frequency. Since the constant damping ζ has little influence to the frequency of the system having a small variation of stiffness, the response frequency in the iteration cycle N has a small difference from that in cycle N-1 due to a small stiffness variation:(22)ωN=ωN-1+Δω.Similarly the angular value of response functionTP∣NωN has a small difference from TP∣N-1ωN-1:(23)TP∣NωN=TP∣N-1ωN-1+eΔω.In this study, the first bending mode of the wing is taken as an approximateωN since it is in the same frequency range and hence dominates the response to the gust. This approximation may produce a small deviation eΔω for the TP∣N prediction. Equation (23) can be further simplified as below when this deviation is initially ignored for small Δω:(24)TP∣NωN=TP∣N-1ωN-1.During optimization process, the angleTP∣N-1ωN-1 in cycle N-1 is obtained from the traditional dynamic response analysis. For the cycle N and following circles, the time TP at the peak gust response can be calculated based on (24). If TP∣NωNR represents the correct value, the deviation of the above predicted result TP∣NωNP due to the eΔω can be evaluated by(25)eTP∣NωNP=1-TP∣NωNPTP∣NωNR.Since theeΔω may increase with the variation of structure stiffness during the optimization, it is also necessary to obtain the correct gust response by performing a complete gust response analysis. The correct value will be used for calibrating the TP∣N+nωN+n to ensure the optimization converge. This is part of the “calibration process” similar to the TP∣N+n prediction process in the Q-Static method.IfXE∣P and XE∣R represent the predicted and the correct gust response values, respectively, the deviation of the predicted result due to all the approximations can be evaluated by(26)eP=1-XE∣PXE∣R.


## 2.3. Gust Response Analysis

The equations used in gust response analysis can be expressed in a simple form as [20](27)MX¨+DX˙+KX=FsteadyX+FunsteadyX,X˙,X¨X0=XN-1X˙0=X˙N-1,in whichM is mass matrix and D is structural damping matrix. Unsteady aerodynamic forces are also included in simulation, which is calculated based on Theodorsen function. The structural analysis is based on thin-walled wing box beam model with the stiffness of each section depending upon the fiber orientations of the skin and spar web laminate. A program “AeroBeamSaw” for static aeroelastic analysis and model analysis is linked with the gust response analysis performed in MATLAB and an optimization program performed in FORTRAN environment [12, 21].


## 2.4. Optimization Method and Objective Function

A gradient-based deterministic optimization method was employed to tailor a composite wing structure for minimum gust loads with the objective function expressed below:(28)MinimizeFdα=LSαLSα0η×XEαXEα0ξ2,in whichLS represents total shear force at wing root, XE represents the peak value of wingtip bending deflection, and αo and α represent a set of initial and tailored design variables respectively, that is, the fiber orientations in a range of 0~±90° of each layer of the laminated composite wing structure. In (28), η and ξ are weight coefficients, η is set at 5, and ξ is set at 0.15.


## 2.5. Calibration Process in the Optimization

In optimization procedure, there are two major parts, which are shown in Figure2. The first major part is Q-S method, which is used to predict wingtip deflection in this article. The second major part is calibration process, which is used to eliminate the error accumulation due to the simplification of the Q-S method.Figure 2
The flowchart of the Q-S method in structure optimization procedure.


### 2.5.1. Q-S Method

Q-S method developed in the current investigation is based on a relationship between the dynamic response of a structure and the work done by an external force. For structural optimization, the relationship is extended to the responding deflections when the structure stiffness has a small variation between two iteration cycles as expressed in (12) and (24).At the starting point of the optimization when the iteration cycleN-1(N-1=0), a dynamic response analysis has to be performed for the initial structure design. Based on the results obtained in the cycle N-1, the deflection of the structure of a small stiffness variation in the iteration cycle N and subsequent cycles N+n(n=1,2,3,…) can be predicted.In this process, instead of performing a dynamic response analysis like the starting point, Q-S method should be used in each iteration cycle repeatedly by recalculating the parametersKE∣N, TP∣N and the Q-S deflection XE∣N-1. It means that the results from one dynamic response analysis can be used for number of iteration cycles with the parameters in (12) and (24) kept updating.


### 2.5.2. Calibration Process

The calibration process in current study is designed on the purpose of improving prediction accuracy. There will be error accumulated in iteration cycles. The main cause of error comes from the variation of structure stiffness, which will affect (12) and (24), as error terms that have connection with small changes made to stiffness are omitted in order to get a simplified recursive relation. When calibration is launched, a traditional dynamic response analysis will be performed to update parameters used in recursive equations for the subsequent optimization process.The launch condition of calibration process is vital for the performance of Q-S method in optimization. There are two major problems.The first problem is that, during optimization procedure, with only a few iteration cycles having “real” data, there is no reliable direct observation method for the assessment of prediction accuracy. Indirect measurement will be taken to set the launch condition. In the current study, equivalent stiffnessKE∣N obtained from static aeroelastic analysis is used as the indirect measurement to trigger calibration process.The second problem is how much indirect measurement varies that can make it necessary to do the calibration. If the conditions are set too strict, it may lower down the efficiency of Q-S method. Otherwise, it may deteriorate Q-S method’s accuracy.In the current study, 4% stiffness variation is set as the criteria for performing the calibration process during optimization.


## 2.5.1. Q-S Method

Q-S method developed in the current investigation is based on a relationship between the dynamic response of a structure and the work done by an external force. For structural optimization, the relationship is extended to the responding deflections when the structure stiffness has a small variation between two iteration cycles as expressed in (12) and (24).At the starting point of the optimization when the iteration cycleN-1(N-1=0), a dynamic response analysis has to be performed for the initial structure design. Based on the results obtained in the cycle N-1, the deflection of the structure of a small stiffness variation in the iteration cycle N and subsequent cycles N+n(n=1,2,3,…) can be predicted.In this process, instead of performing a dynamic response analysis like the starting point, Q-S method should be used in each iteration cycle repeatedly by recalculating the parametersKE∣N, TP∣N and the Q-S deflection XE∣N-1. It means that the results from one dynamic response analysis can be used for number of iteration cycles with the parameters in (12) and (24) kept updating.


## 2.5.2. Calibration Process

The calibration process in current study is designed on the purpose of improving prediction accuracy. There will be error accumulated in iteration cycles. The main cause of error comes from the variation of structure stiffness, which will affect (12) and (24), as error terms that have connection with small changes made to stiffness are omitted in order to get a simplified recursive relation. When calibration is launched, a traditional dynamic response analysis will be performed to update parameters used in recursive equations for the subsequent optimization process.The launch condition of calibration process is vital for the performance of Q-S method in optimization. There are two major problems.The first problem is that, during optimization procedure, with only a few iteration cycles having “real” data, there is no reliable direct observation method for the assessment of prediction accuracy. Indirect measurement will be taken to set the launch condition. In the current study, equivalent stiffnessKE∣N obtained from static aeroelastic analysis is used as the indirect measurement to trigger calibration process.The second problem is how much indirect measurement varies that can make it necessary to do the calibration. If the conditions are set too strict, it may lower down the efficiency of Q-S method. Otherwise, it may deteriorate Q-S method’s accuracy.In the current study, 4% stiffness variation is set as the criteria for performing the calibration process during optimization.


## 3. Example and Results

### 3.1. An Aircraft Planform and the Wing Structural Model

An aircraft of tailless flying-wing configuration as shown in Figure3 is taken as an example to demonstrate the performance of the Q-S method. The planform and dimension of the aircraft outer wing of 20° swept angle are shown in Figure 4.Figure 3
Planform of a flying wing aircraft model (unit: m).Figure 4
Planform of composite outer wing modeled in eight single-cell wing-box beams.The primary structure of the wing made of laminated carbon/epoxy composite is modeled by assembly of eight segments of thin-walled wing-box beams enclosed between the spars. The leading and trailing edge cells of the wing were considered in mass, inertia, and aerodynamic force calculation. Each of the single-cell sections was further divided into four laminate panels representing the skins and spar webs of the wing-box shown in Figure5.Figure 5
Cross section details of a wing-box beam element.Each of the panels is made of 8-layer carbon/epoxy laminate of symmetric layup. The composite material properties in the usual notation areE1=220.5 GPa, E2=17.7 GPa, G12=7.05 GPa, and density ρ=1900 Kg/m3. Taking the fiber orientation of each ply as a design variable, the total number of independent design variables is 128 in the optimization of the outer wing structure. If the layups of eight span-wise segments were kept uniform, the design variables can be reduced to 16 in the optimization.In this example, the outer wing as shown in Figure4 is clamped at the inboard root section and subjected to the aerodynamic force. The wing box of original laminate layup was optimized to minimize the gust response. The specified gust velocity Vds is 6 m/s and air speed is 80 m/s. The initial design and the optimized laminate layups and response results of the wing box structure are shown in Table 1. The shear force at wing root has been reduced 10.8% from original 142.7 KN to 127.2 KN and torque at wing root has been reduced from −87.7 KN·m to −78.2 KN·m. Meanwhile, the gust response in terms of wingtip deflection of the optimized wing has been reduced by 16.3% from original 0.43 m to 0.36 m as shown in Table 1.Table 1
Initial and optimized laminate layups, stiffness, and deflection.
Status Panels Laminate layup (deg) Stiffness and deformation Initial design Upper skin [0.0/45.0/−45.0/90.0] s EI = 6.80 ⋅105 KN⋅m2 Lower skin [0.0/45.0/−45.0/90.0] s GJ = 3.77 ⋅ 105 KN⋅m2 Upper/lower skin thickness 2 mm for each layer Wingtip deflection: 0.32 m Front spar web [0.0/45.0/−45.0/90.0] s Wingtip twist: −0.94 deg Rear spar web [0.0/45.0/−45.0/90.0] s Shear force at wing root: 142.7 KN Front/rear spar thickness 6 mm for each layer Torque at wing root: −87.7 KN⋅m  Optimized design Upper skin [0.0/76.2/−84.1/90.0] s EI = 5.49 ⋅ 105 KN⋅m2 Lower skin [0.0/76.2/−84.1/90.0] s GJ = 1.09 ⋅ 105 KN⋅m2 Upper/lower skin thickness 2 mm for each layer Wingtip deflection: 0.36 m Front spar web [0.0/89.9/−89.9/90.0] s Wingtip twist: −1.70 deg Rear spar web [0.0/89.9/−89.9/90.0] s Shear force at wing root: 127.2 KN Front/rear spar thickness 6 mm for each layer Torque at wing root: −78.2 KN⋅mThe optimization process history can be divided into three stages as shown in Figure6. In stage A, the gradient of objective function to each of the design variables is calculated. The variation of the objective function value was very small due to the small change of design variables and resulting Δω at each step. In stage C, the variation of design variables and resulting Δω was also small when the variation converged to the final optimized results. In stage B however, the design variables and resulting Δω varied in relatively large steps in the direction determined by the gradient. Therefore the necessary dynamic response analysis to update the parameters in (19) is normally required in stage B.Figure 6
The optimization process using gradient-based deterministic method.


### 3.2. Reduced Computation Time in Optimization

Table2 shows the computational time used for calculating the wing response to the gust by using the Q-S method and the traditional dynamic response analysis in the same optimization case. A relatively small percent of time was used by the “other calculation” including the gradient calculation and fiber orientation tailoring.Table 2
CPU time distribution in optimization process.
Time distribution Q-static method Dynamic response Other calculation Time for a cycle 1.922 s 81.03 s 0.79 s Total cycles 100 6 100 Total time used 192.2 s 486.2 s 79.0 s Share of time 25.4% 64.2% 10.4%In this case, the total number of iteration cycles in the optimization process is 100 with total 12.6 minutes of computational time. Although the traditional dynamic response analysis was performed only 6 times, it took 64.2% of the total computation time. If the traditional dynamic response analysis were performed in all 100 iteration cycles in the traditional approach, the total time to complete the optimization would be 135 minutes. The computation time saving is about 91% by employing the Q-S method.


## 3.1. An Aircraft Planform and the Wing Structural Model

An aircraft of tailless flying-wing configuration as shown in Figure3 is taken as an example to demonstrate the performance of the Q-S method. The planform and dimension of the aircraft outer wing of 20° swept angle are shown in Figure 4.Figure 3
Planform of a flying wing aircraft model (unit: m).Figure 4
Planform of composite outer wing modeled in eight single-cell wing-box beams.The primary structure of the wing made of laminated carbon/epoxy composite is modeled by assembly of eight segments of thin-walled wing-box beams enclosed between the spars. The leading and trailing edge cells of the wing were considered in mass, inertia, and aerodynamic force calculation. Each of the single-cell sections was further divided into four laminate panels representing the skins and spar webs of the wing-box shown in Figure5.Figure 5
Cross section details of a wing-box beam element.Each of the panels is made of 8-layer carbon/epoxy laminate of symmetric layup. The composite material properties in the usual notation areE1=220.5 GPa, E2=17.7 GPa, G12=7.05 GPa, and density ρ=1900 Kg/m3. Taking the fiber orientation of each ply as a design variable, the total number of independent design variables is 128 in the optimization of the outer wing structure. If the layups of eight span-wise segments were kept uniform, the design variables can be reduced to 16 in the optimization.In this example, the outer wing as shown in Figure4 is clamped at the inboard root section and subjected to the aerodynamic force. The wing box of original laminate layup was optimized to minimize the gust response. The specified gust velocity Vds is 6 m/s and air speed is 80 m/s. The initial design and the optimized laminate layups and response results of the wing box structure are shown in Table 1. The shear force at wing root has been reduced 10.8% from original 142.7 KN to 127.2 KN and torque at wing root has been reduced from −87.7 KN·m to −78.2 KN·m. Meanwhile, the gust response in terms of wingtip deflection of the optimized wing has been reduced by 16.3% from original 0.43 m to 0.36 m as shown in Table 1.Table 1
Initial and optimized laminate layups, stiffness, and deflection.
Status Panels Laminate layup (deg) Stiffness and deformation Initial design Upper skin [0.0/45.0/−45.0/90.0] s EI = 6.80 ⋅105 KN⋅m2 Lower skin [0.0/45.0/−45.0/90.0] s GJ = 3.77 ⋅ 105 KN⋅m2 Upper/lower skin thickness 2 mm for each layer Wingtip deflection: 0.32 m Front spar web [0.0/45.0/−45.0/90.0] s Wingtip twist: −0.94 deg Rear spar web [0.0/45.0/−45.0/90.0] s Shear force at wing root: 142.7 KN Front/rear spar thickness 6 mm for each layer Torque at wing root: −87.7 KN⋅m  Optimized design Upper skin [0.0/76.2/−84.1/90.0] s EI = 5.49 ⋅ 105 KN⋅m2 Lower skin [0.0/76.2/−84.1/90.0] s GJ = 1.09 ⋅ 105 KN⋅m2 Upper/lower skin thickness 2 mm for each layer Wingtip deflection: 0.36 m Front spar web [0.0/89.9/−89.9/90.0] s Wingtip twist: −1.70 deg Rear spar web [0.0/89.9/−89.9/90.0] s Shear force at wing root: 127.2 KN Front/rear spar thickness 6 mm for each layer Torque at wing root: −78.2 KN⋅mThe optimization process history can be divided into three stages as shown in Figure6. In stage A, the gradient of objective function to each of the design variables is calculated. The variation of the objective function value was very small due to the small change of design variables and resulting Δω at each step. In stage C, the variation of design variables and resulting Δω was also small when the variation converged to the final optimized results. In stage B however, the design variables and resulting Δω varied in relatively large steps in the direction determined by the gradient. Therefore the necessary dynamic response analysis to update the parameters in (19) is normally required in stage B.Figure 6
The optimization process using gradient-based deterministic method.


## 3.2. Reduced Computation Time in Optimization

Table2 shows the computational time used for calculating the wing response to the gust by using the Q-S method and the traditional dynamic response analysis in the same optimization case. A relatively small percent of time was used by the “other calculation” including the gradient calculation and fiber orientation tailoring.Table 2
CPU time distribution in optimization process.
Time distribution Q-static method Dynamic response Other calculation Time for a cycle 1.922 s 81.03 s 0.79 s Total cycles 100 6 100 Total time used 192.2 s 486.2 s 79.0 s Share of time 25.4% 64.2% 10.4%In this case, the total number of iteration cycles in the optimization process is 100 with total 12.6 minutes of computational time. Although the traditional dynamic response analysis was performed only 6 times, it took 64.2% of the total computation time. If the traditional dynamic response analysis were performed in all 100 iteration cycles in the traditional approach, the total time to complete the optimization would be 135 minutes. The computation time saving is about 91% by employing the Q-S method.


## 4. Further Study of the Q-S Method Accuracy

Theoretically the accuracy and accuracy of the Q-S method largely depend upon two major factors when used in structural optimization. The first one is the accuracy ofTP∣N prediction or the amount of e(Δω) as shown in (24). The second one is the amount of ΔW/EEP∣N-1 given in (10). Since both factors are related to the stiffness variation step size, the Q-S method accuracy also depends on the step size.In practice, the frequency difference between the first mode of the wing structure and gust load also affects the accuracy of the predictedTP∣N. For a higher gust frequency, the wing structure higher order modes have more influence on gust response. The aircraft free rigid-body modes have also influence on the prediction accuracy. As rigid-body mode is much lower than the gust frequency, the influence is negligible.To evaluate the accuracy of the method, two parametric study cases were carried out and presented in the following subsections. The two major factors were also observed in both two groups, as the relative relation between the variation of prediction errors and the change of stiffness is the main object in the study.


### 4.1. Parametric Study of the Q-S Method

The first case is to quantify the influence of different gust frequency on the accuracy ofTP∣N prediction. Since the gust frequency varies with the flight speed, different flight speed from 75 m/s to 150 m/s corresponding to the gust frequency range from 0.75 Hz to 1.50 Hz as specified in CS-25 was considered. The first three modal frequencies of the wing structure 1.17 Hz, 2.29 Hz, and 5.61 Hz, respectively, were taken into account in the analysis. By using the Q-S method, the TP∣NωN(N=1,2,3,…) results based on the predicted TP∣N can be calculated.In order to evaluate the accumulated deviation during the optimization, the calibration process in the Q-S method was locked in this case study. The deviation of predicted angle as defined in (25) from the correct value obtained from complete gust response analysis in the optimization process is shown in Figure 7. Only limited number of cycle steps was taken to calculate the deviation where a complete dynamic analysis was performed to obtain the correct value as reference. The results show that the deviation of predicted TP∣NωN is less than 4% for all four different gust frequencies in the CS25 specified range. For the higher gust frequencies, the prediction deviation increases because of the influence of higher order modes of the wing structure.Figure 7
The deviation of predicted angles without calibration process for different gust frequencies.When the proposed calibration process was performed, the key parameters affecting the accuracy are updated during the optimization. There are four cases investigated here, with different gust frequencies, the results are shown in Figure8. The resulting deviation of the gust response as defined by (26) for the four different gust frequencies is shown in Figure 9. Again the correct gust response was obtained from performing complete gust response analysis at each optimization cycle, and the gust response was measured as wingtip deflection.Figure 8
Predicted and real results under different gust frequencies with calibration process.Figure 9
The deviation of predicted gust response with calibration process for different gust frequencies.As shown in Figure9, the accumulated maximum deviation is nearly 5%. When the calibration process was performed, the deviation was reduced to zero at two iteration cycle steps as circled in Figure 9 during the optimization. The results show the effectiveness of updating the predicted TP∣N and XE∣P by performing the calibration process in the Q-S method. The resulting deviation of the gust response at the end of optimization was reduced to about 1%. Table 3 lists the deviation of the predicted gust response against the computation time saving compared with the correct results obtained from performing the complete gust response analysis at each step of the optimization.Table 3
Accuracy and computation time saving of the Q-static method in different gust frequencies.
Gust frequency 0.75 Hz 1.0 Hz 1.25 Hz 1.50 Hz Accuracy reduction 1.1% 1.5% 2.1% 2.0% Total cycles 100 100 100 100 Time reduced by ~90% ~90% ~90% ~90%


### 4.2. Case Studies of Aircraft Response Including Rigid-Body Motion

The second study case is to evaluate the influence of rigid-body mode on the prediction accuracy. In this case, a constant gust frequency corresponding to flight speed 150 m/s was taken. In this case, the gust response of the whole aircraft as shown in Figure3 is taken in the gust response analysis. The aircraft was set in a trimmed level flight condition with free motion in the symmetric plane before entering the (1-cos) gust load uniformly in span-wise direction. The predicted and real deflection results at specific optimization points are shown in Figure 10, compared with results with clamped wing-root. Figure 11 shows the deviation of the predicted gust response in the optimization of the wing structure with and without the aircraft rigid-body mode. In the analysis, the gust response was measured as the wingtip deflection relative to the wing root. The results in Table 4 show that the rigid-body mode of the aircraft has little effect on the accuracy of the prediction using the Q-S method for the wing optimization.Table 4
Prediction method’s performance in trimmed free-flight condition.
Gust model (1-cos) gust Response reduction 36.7% Accuracy reduction 2%~4% Total cycles 80 Time saved 90%Figure 10
Predicted and real result for the cases with and without the aircraft rigid-body mode.Figure 11
Deviation of the prediction gust response for the cases with and without the aircraft rigid-body mode.


## 4.1. Parametric Study of the Q-S Method

The first case is to quantify the influence of different gust frequency on the accuracy ofTP∣N prediction. Since the gust frequency varies with the flight speed, different flight speed from 75 m/s to 150 m/s corresponding to the gust frequency range from 0.75 Hz to 1.50 Hz as specified in CS-25 was considered. The first three modal frequencies of the wing structure 1.17 Hz, 2.29 Hz, and 5.61 Hz, respectively, were taken into account in the analysis. By using the Q-S method, the TP∣NωN(N=1,2,3,…) results based on the predicted TP∣N can be calculated.In order to evaluate the accumulated deviation during the optimization, the calibration process in the Q-S method was locked in this case study. The deviation of predicted angle as defined in (25) from the correct value obtained from complete gust response analysis in the optimization process is shown in Figure 7. Only limited number of cycle steps was taken to calculate the deviation where a complete dynamic analysis was performed to obtain the correct value as reference. The results show that the deviation of predicted TP∣NωN is less than 4% for all four different gust frequencies in the CS25 specified range. For the higher gust frequencies, the prediction deviation increases because of the influence of higher order modes of the wing structure.Figure 7
The deviation of predicted angles without calibration process for different gust frequencies.When the proposed calibration process was performed, the key parameters affecting the accuracy are updated during the optimization. There are four cases investigated here, with different gust frequencies, the results are shown in Figure8. The resulting deviation of the gust response as defined by (26) for the four different gust frequencies is shown in Figure 9. Again the correct gust response was obtained from performing complete gust response analysis at each optimization cycle, and the gust response was measured as wingtip deflection.Figure 8
Predicted and real results under different gust frequencies with calibration process.Figure 9
The deviation of predicted gust response with calibration process for different gust frequencies.As shown in Figure9, the accumulated maximum deviation is nearly 5%. When the calibration process was performed, the deviation was reduced to zero at two iteration cycle steps as circled in Figure 9 during the optimization. The results show the effectiveness of updating the predicted TP∣N and XE∣P by performing the calibration process in the Q-S method. The resulting deviation of the gust response at the end of optimization was reduced to about 1%. Table 3 lists the deviation of the predicted gust response against the computation time saving compared with the correct results obtained from performing the complete gust response analysis at each step of the optimization.Table 3
Accuracy and computation time saving of the Q-static method in different gust frequencies.
Gust frequency 0.75 Hz 1.0 Hz 1.25 Hz 1.50 Hz Accuracy reduction 1.1% 1.5% 2.1% 2.0% Total cycles 100 100 100 100 Time reduced by ~90% ~90% ~90% ~90%


## 4.2. Case Studies of Aircraft Response Including Rigid-Body Motion

The second study case is to evaluate the influence of rigid-body mode on the prediction accuracy. In this case, a constant gust frequency corresponding to flight speed 150 m/s was taken. In this case, the gust response of the whole aircraft as shown in Figure3 is taken in the gust response analysis. The aircraft was set in a trimmed level flight condition with free motion in the symmetric plane before entering the (1-cos) gust load uniformly in span-wise direction. The predicted and real deflection results at specific optimization points are shown in Figure 10, compared with results with clamped wing-root. Figure 11 shows the deviation of the predicted gust response in the optimization of the wing structure with and without the aircraft rigid-body mode. In the analysis, the gust response was measured as the wingtip deflection relative to the wing root. The results in Table 4 show that the rigid-body mode of the aircraft has little effect on the accuracy of the prediction using the Q-S method for the wing optimization.Table 4
Prediction method’s performance in trimmed free-flight condition.
Gust model (1-cos) gust Response reduction 36.7% Accuracy reduction 2%~4% Total cycles 80 Time saved 90%Figure 10
Predicted and real result for the cases with and without the aircraft rigid-body mode.Figure 11
Deviation of the prediction gust response for the cases with and without the aircraft rigid-body mode.


## 5. Conclusions

This paper presents a quasi-static method to predict the peak value of gust response of a wing structure made of composite efficiently. Significant computational time can be saved when applying this method to optimize the structural laminate layup subjected to gust load. An aircraft of flying-wing configuration has been taken as an example to demonstrate the effectiveness of the method when applied to minimizing response to a discrete (1-cos) gust by optimizing the wing box structure. The validation of the example is shown by the results that the gust response has been reduced by 22% after optimizing the laminate layup of the composite wing structure. By employing the Q-static method, the computational time is reduced by over 90% comparing with the results from conducting the usual complete gust response analysis at each iterative cycle of the optimization. Although the prediction has to be calibrated by conducting the usual gust response analysis when the deviation is over a limit, the method is valid in most of the iterative cycles. In the total 80 times of iteration cycles of the example, only four times to perform the usual gust response analysis were required to update the result by the Q-S method. The parametric study has shown that the Q-static method can be applied to the whole gust frequency range as specified in CS-25. The deviation of the predicted results by the method is below 2.1% with the computation time saving nearly 90%. The rigid-body mode of the aircraft has little effect on the accuracy of the Q-S method.



---

*Source: 1592527-2017-01-26.xml*