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.
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 energy includes both elastic potential and kinetic energy written as [18]where is the elastic deformation vector and is the stiffness matrix of a structure.
An estimation of the stiffness of a structure, , can be obtained by the force applied against the deflection measured at representative points. In this article, an estimation of stiffness 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 is used in the following Q-Static method for gust response analysis.
Assuming 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 at . At this moment, the velocity and kinetic energy of the structure becomes zero [19]. The system is in Q-Static status and total energy can be simplified to an expression of elastic potential energy .
The work done by the external dynamic force through the structural deformation over a time period until a specified time such as the peak value of dynamic response can be written aswhere is the corresponding velocity of the structural response to .
It is noted that the gives an approximate total work done based on the estimated . For a small variation of structural 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 cycle can be represented by the stiffness at previous cycle :
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:
This leads to the relationship of the between two iteration cycles:
As structural stiffness is varied, the time to reach the response peak value is also changed from the previous iteration cycle:
in which represents the time difference of peak gust response between cycle and . The work done by the external force between two iteration cycles is related by
in which , and is the Q-S ratio of response amplitude in two contiguous iteration cycles given by
From (6) and (8), the following relation between and is obtained in two contiguous iteration cycles:where represents a deviation of the ratio between two contiguous iteration cycles.
The 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 below a value, it can be negligible as discussed in Section 2.5. In this case, a simplified ratio can be obtained from (10):
When a ratio is obtained at iteration step , an approximate at iteration cycle 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 () discrete gust model as defined in airworthiness regulation CS-25 is considered:where is a specified design gust velocity and is the gust frequency, in which is the penetrating distance of an aircraft at an equivalent flight speed into a discrete gust at time . is gust gradient distance.
The symmetric aerodynamic force due to the gust can be calculated bywhere is the air density, is the reference wing area, and is the lift curve slope of the wing. is the gust alleviating factor expressed by
It is noted that the above is defined for a typical gust gradient distance where 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 cases. Corresponding to the () gust, the responding wingtip deflection can also be expressed as a () function. Hence the response velocity before approaching the peak time can be expressed in a sinusoidal function by differentiating the wingtip deflection:where represents the wingtip deflection amplitude at the response time TP. Based on (3), (14), and (16), the can be rewritten as
According to (2) and (17), (12) can be rewritten as
By extracting the parameters , , , , and out of the time integration and substituting (9) into (18), a simple recursive relation of at can be obtained from the as shown:
In (19), since the in the iteration cycle is different from previous iteration cycle , it will be predicted in each iteration cycle. The details of how to predict the are presented in the following section.
Although the deviation as shown in (10) may be negligible, it will accumulate when (19) is used to predict the in subsequent cycles beyond cycle N. 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 to reach its peak value as illustrated in Figure 1 can be expressed in (19), in which is the structural natural frequency and the damping ratio:
An expression of the angular value at peak response time in general form can be written aswhere 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 has a small difference from that in cycle due to a small stiffness variation:
Similarly the angular value of response function has a small difference from :
In this study, the first bending mode of the wing is taken as an approximate since it is in the same frequency range and hence dominates the response to the gust. This approximation may produce a small deviation for the prediction. Equation (23) can be further simplified as below when this deviation is initially ignored for small :
During optimization process, the angle in cycle is obtained from the traditional dynamic response analysis. For the cycle and following circles, the time at the peak gust response can be calculated based on (24). If represents the correct value, the deviation of the above predicted result due to the can be evaluated by
Since the 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 to ensure the optimization converge. This is part of the “calibration process” similar to the prediction process in the Q-Static method.
If and 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
2.3. Gust Response Analysis
The equations used in gust response analysis can be expressed in a simple form as [20]
in which is mass matrix and 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:
in which represents total shear force at wing root, represents the peak value of wingtip bending deflection, and 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 Figure 2. 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.
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 cycle , a dynamic response analysis has to be performed for the initial structure design. Based on the results obtained in the cycle , the deflection of the structure of a small stiffness variation in the iteration cycle and subsequent cycles 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 parameters , and the Q-S deflection . 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 stiffness 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 Figure 3 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.
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 Figure 5.
Each of the panels is made of 8-layer carbon/epoxy laminate of symmetric layup. The composite material properties in the usual notation are GPa, GPa, GPa, and density 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 Figure 4 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 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.
The optimization process history can be divided into three stages as shown in Figure 6. 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.
3.2. Reduced Computation Time in Optimization
Table 2 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.
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 of prediction or the amount of as shown in (24). The second one is the amount of 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 predicted . 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 of 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 results based on the predicted 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 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.
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 Figure 8. 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.
As shown in Figure 9, 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 and 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.
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 Figure 3 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 () 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.
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 () 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.
Nomenclature
| : | Total energy of a structure’s deformation movement |
| : | Kinetic energy of a structure’s deformation movement |
| : | Elastic potential energy of a structure’s deformation movement |
| : | Stiffness matrix of a structure |
| : | Estimation of a structure’s stiffness |
| : | Elastic deformation vector |
| : | First peak time of a structure’s dynamic response |
| : | A type of dynamic response of a structure at |
| : | Time function of external dynamic force |
| : | Time function of corresponding velocity of structural response to |
| : | Work done by external force through structural deformation over a time period |
| : | Amplitude of structural deflection |
| : | Structural natural frequency |
| : | Total shear force at wing root. |
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors acknowledge Chinese Natural Science Funding Council (NSFC 11302027, NSFC 11672033) for financial support of the project.