Abstract
This paper presents a framework of aeroelastic optimization design for high-aspect-ratio wing with large deformation. A highly flexible wing model for wind tunnel test is optimized subjected to multiple aeroelastic constraints. Static aeroelastic analysis is carried out for the beamlike wing model, using a geometrically nonlinear beam formulation coupled with the nonplanar vortex lattice method. The flutter solutions are obtained using the - method based on the static equilibrium configuration. The corresponding unsteady aerodynamic forces are calculated by nonplanar doublet-lattice method. This paper obtains linear and nonlinear aeroelastic optimum results, respectively, by the ISIGHT optimization platform. In this optimization problem, parameters of beam cross section are chosen as the design variables to satisfy the displacement, flutter, and strength requirements, while minimizing wing weight. The results indicate that it is necessary to consider geometrical nonlinearity in aeroelastic optimization design. In addition, optimization strategies are explored to simplify the complex optimization process and reduce the computing time. Different criterion values are selected and studied for judging the effects of the simplified method on the computing time and the accuracy of results. In this way, the computing time is reduced by more than 30% on the premise of ensuring the accuracy.
1. Introduction
Highly flexible wings, crucial for high-altitude long-endurance (HALE) unmanned aerial vehicles (UAVs), are characterized by light weight with high aspect ratios. Although the slender wings can maximize the lift-to-drag ratio, they may undergo large deformation of about 25% of wing semispan during normal flight load, exhibiting geometrically nonlinear behaviors. Patil and Hodges [1] studied the aeroelasticity of a HALE aircraft and found that large deformation induced by flexibility will change the aerodynamic load distribution. This will result in significant changes in aeroelastic and flight dynamic responses of the aircraft. Thus, the aeroelastic analysis based on small deformation assumption may be improper. Shearer and Cesnik [2] also studied the nonlinear dynamic response of a flexible aircraft. The nonlinear structural response is governed by the 6-DOF rigid-body motions coupled with aeroelastic equations. They compared three solutions: rigid, linearized, and nonlinear models and highlighted the use of the latter formulation to model the flexible wings. All these studies show that the geometrical nonlinearity due to the flexibility should be properly accounted for a nonlinear aeroelastic formulation. Such flexible wings can be treated as a statically nonlinear but dynamic linear system. Then the question of the dynamic stability of the statically nonlinear aeroelastic system may be addressed by a linear dynamic perturbation analysis about this nonlinear static equilibrium [3–5]. Patil and Hodges [6] studied the nonlinear effects and concluded that the dominant geometrically nonlinear effect comes from the nonzero steady-state curvature. Thereby, a linear analysis for curved wing structures can accurately predict the effect of nonlinear behavior on the linearized stability.
Optimization design of aircraft structures subjected to aeroelastic constraints is not new now in preliminary stage of modern aircraft structural design. The aeroelastic optimization process, in essence, is a unified, multidisciplinary design process which can overcome barriers between different discipline groups while reducing the design cycle time. Previous researches have mainly optimized aircraft structures from the perspectives of both static and dynamic constraints such as stress, displacement, modal frequency, and flutter constraints; see, for example, the previous publications by Sikes et al. [7], Patil [8], Battoo and de Visser [9], Maute and Allen [10], Tischler, and Venkayya [11]. However, these researches have not taken into account the geometric nonlinearity induced by large deformation. It may lead to an inaccurate modeling of aircraft in both structural and aerodynamic perspectives. Butler et al. [12, 13] presented aeroelastic optimization of high-aspect-ratio wings subjected to a minimum flutter speed constraint based on a dynamic stiffness matrix method. Yet, the basic assumption in solving the dynamic problems is linear. Previous researches have already shown that geometrical nonlinearities must be properly considered in aeroelastic analysis for flexible wings. In general, few studies have attempted to deal with the geometric nonlinearities in aeroelastic optimization design for flexible wings with large deformation. In this paper, numerical studies will demonstrate the necessity of considering nonlinear effects in aeroelastic optimization. Moreover, nonlinear aeroelastic optimization strategies are explored in order to simplify the complex optimization process while reducing the computing time.
The main purpose of the present work is to establish a complete framework of aeroelastic optimization design for high-aspect-ratio wing with large deformation. The geometrical nonlinearities in present theoretical modeling are considered from three aspects: structural dynamic, aerodynamics, and fluid-structure interface. The incremental finite element method is used to solve the geometrically nonlinear structural problems, and aerodynamic loads on nonplanar wings are calculated by the nonplanar vortex lattice method. The unsteady aerodynamic forces in flutter analysis are calculated by the nonplanar doublet-lattice method. A surface spline interpolation method [14] is used to exchange force and displacement information between the structure model and aerodynamic model. In this way, the aeroelastic formulation is naturally obtained by combining these three aspects together. A wing model for wind tunnel test is optimized subject to displacement, torsion, flutter, and strength constraints. The nonlinear aeroelastic analysis, as a key part in optimization design, is integrated into the whole optimization process to obtain the static equilibrium configuration and flutter speed. To improve efficiency of nonlinear aeroelastic optimization process, some simplifications in static and dynamic aeroelastic analysis are studied and the results indicate that such changes can greatly reduce the computing time.
2. Theoretical Formulations
A complete theoretical framework of nonlinear aeroelastic analysis has been discussed in [3, 6]. An introduction is presented here, followed by the optimization formulation for achieving the optimum wing shape.
2.1. Geometrical Nonlinear Elasticity
Geometric nonlinearities are based on the kinematic description of the body and the strain on the wing should be defined in terms of local displacement of the wing for dynamic motions. These result in the nonlinear geometric equations including the quadric term of the displacement differential, and the nonlinear force equilibrium equation established on the deformed state of the structure. Geometric nonlinear effects are prominent in two different aspects: geometric stiffening due to initial displacements and stresses and follower forces due to a change in loads as a function of displacements. Both of these two factors are considered in this paper and solved with the nonlinear incremental finite element method [15]. The updated Lagrange formulation (ULF) is used in this work and the main equations are presented below.
The relationship between strain and displacement iswhere means the partial derivative of displacement with respect to the coordinate at time . The stress tensor at time satisfieswhere is direction cosine of small aeroelement at time and is the corresponding surface force. The linear elastic constitution can be described as follows:where is the elastic tensor, which has a different form for isotropic or anisotropic material.
The strain can be decomposed into a linear part and a nonlinear part :
The stress is decomposed by increments, where represents the stress at time , and is the incremental stress to be calculated at each time step.
The integral equation is established by linearization in each incremental step:where is a vector of incremental external force, including the aerodynamic force, gravity, and engine thrust at time , and is the element volume. Using the following shape functions, the relationship between strain and displacements is performed as follows:
Substituting them into (6) yields the element governing equation for static problems:where is the equivalent inner force of the structure. The stiffness matrix can be decomposed into a linear part and nonlinear part . The nonlinear part is related to the deformed configuration, load condition, and strain, which should be updated in each computation step.
The corresponding dynamic equation can be expressed as follows: where is the instant mass matrix at time . and are the structural displacement vector and acceleration vector, respectively.
According to previous researches [1, 3], an assumption of small-amplitude vibration around nonlinear static equilibrium state is suitable for aeroelastic stability problem such as flutter analysis. Ignoring the damping effect, there is where is the large static deformation from (8) and is a small vibration deformation. According to (9) and the static equilibrium condition, the vibration equation of the dynamic system reduces towhere is the inertial matrix of the structure at the static equilibrium. is the corresponding stiffness matrix. Both of them are functions of static deformation and vary with different static equilibrium states. From (11), the mode shapes and frequencies are deduced. Combined with nonlinear aerodynamic computation, the flutter boundary can be determined by - method. However, the flutter boundary is only predicted under a certain static equilibrium state. Different static equilibrium states have different flutter boundaries and the exact boundary should be searched iteratively to make the flutter speed be consistent with the static equilibrium state.
2.2. Nonplanar Aerodynamics
Two methods are used to obtain the aerodynamic loads: nonplanar vortex lattice method (VLM) [16] for solving steady aerodynamic loads in static aeroelastic analysis and doublet-lattice method (DLM) [17, 18] for solving unsteady aerodynamic loads in dynamic stability analysis. The lifting surface mesh in both of them can fit the actual wing surface and be consistent with the structural nonlinearities caused by large deformation.
2.2.1. Nonplanar Vortex Lattice Method
In static aeroelastic analysis, the steady aerodynamic load is calculated by nonplanar vortex lattice method, which is based on full potential equations without any linearization. The exact boundary condition is satisfied on the actual wing surface. A Cartesian coordinate system is established for aeroelastic analysis. The -axis points from the nose to the tail along the free stream, the -axis points to the right side, and -axis is determined by the right-hand rule. In order to solve the steady aerodynamic load, the lifting surface is ideally discretized into trapezoidal panels, with two side edges parallel to the undistributed flow as shown in Figure 1. This spatial surface mesh can reflect the actual wing surface that has camber and various platform shapes well. Some typical panel elements are shown in Figure 2. Each vortex ring consists of four segments of a vortex line, and the leading segment is placed on the chord line. The aerodynamics of the panel act on the midpoint of this segment (represented by “” in Figure 2). The collocation point (represented by “” in Figure 2) is located at the center of the three-quarters chord line, and the actual boundary condition will be implemented at this point.
The velocity induced by a typical vortex ring at an arbitrary point can be calculated by applying the Biot-Savart law to the ring’s four segments, except for the rings located at the trailing edge of the wing. Different from normal vortex rings, two semi-infinite vortex lines are shed into the flow along the -axis at each trailing edge panel. In this way, the induced velocity at all of the collocation points could be expressed as follows:where , , and are the components of the induced velocities along the -, -, and -axes. , , and are the corresponding influence coefficients matrices. is the vortex strength vector of all the vortex rings.
According to Neumann boundary condition [19], for the collocation point of the th vortex ring element, there iswhere is the velocity of the free stream. When the angle of sideslip and the angle of attack are small, there is , in which is the velocity magnitude of the free stream. is the induced velocity at the collocation point of the th vortex ring element and . is the normal vector of the th panel element as shown in Figure 1 and . Since the boundary condition is satisfied on the actual wing surface, the normal vector of each panel element depends on the structural deformation.
Expanding (13) yields
Let , and the boundary condition equations can be expressed in a matrix form.where is the coefficient matrix of the normal induced velocity. is the coefficients vector related to the initial angle of attack of each panel element. and are the coefficients vectors related to the angle of sideslip and the angle of attack, respectively. The solution of the above equation results in the vector of vortex strength .
After obtaining the strength of each vortex ring, the aerodynamic force that acts on the ith vortex ring element can be calculated using Kutta-Joukowski theorem:where is the air density. is the total vortex strength at quarter-chord line of the ith panel element and , in which is the vector that describes the magnitude and direction of the th panel element’s chord line. When the panel is located at the leading edge of the wing, , and otherwise, . It should be noted that is not the value of lift per unit span; it is the force acting on the panel element.
Combined with (15), the aerodynamic forces of all the panel elements can be expressed as a column vector:where is the flight dynamic pressure. is the coefficient matrices along the three axes of the aerodynamic coordinate. is the transfer matrix for the vortex strength.
2.2.2. Unsteady Aerodynamic Modeling
The subsonic doublet-lattice method in frequency domain is used to calculate unsteady aerodynamics in flutter analysis. The lifting surface mesh is the same as that in VLM. The doublet line stands along the quarter-chord line of the panel. The control points are fixed to the center of the three-quarters chord line. Different from the VLM, a doublet lattice is placed on the midpoint (represented by in Figure 1) of the chord line of the aeroelement. The collocation point (represented by in Figure 1) is at the center of the three-quarters chord line. At the control point of the th panel element, there is where is the amplitude of normal wash at the control point of the th panel element. The relationship between the pressure coefficient of the th element and the pressure satisfies the following: . and are the length of centerline chord and chord line of the th element, respectively. is the sweep angle of the th element (i.e., the sweep angle of as shown in Figure 1). is the number of panel elements. is the acceleration potential kernel function for oscillatory subsonic flow.
The aerodynamic pressure on panel elements can be written in a matrix form:where is the air density and is the velocity of the free stream. D is the aerodynamic influence coefficient matrix. is the vector of the normal washes at collocation points. By modal analysis, the structural deformation can be expressed by a certain combination of structural modes, which are calculated based on the nonlinear static analysis. To simulate the real situation of large deformation, the boundary conditions are satisfied on the geometrically exact aerodynamic configuration. Therefore, the generalized boundary condition based on the quasimodes around the static equilibrium configuration can be expressed as follows:where is the mode shape projected in the local normal direction and varies with different static equilibrium configurations. Due to the geometrical nonlinearity discussed above, modal shapes and structural stiffness should be updated with the variation of structural deformations when computing unsteady aerodynamic forces. is the reduced frequency, in which is the reference chord length and is the radius frequency.
2.3. Nonlinear Static Aeroelastic and Flutter Analysis
Both of the nonlinear static aeroelastic and flutter analysis are proceeded iteratively, as shown in Figure 3. The static aeroelastic analysis consists of three aspects: nonlinear static analysis, calculation of the steady aerodynamic force, and coupling iterative computations between the structure and aerodynamic surface. In this paper, a loosely coupling strategy is used to solve the static aeroelastic problem via the nonplanar vortex lattice method under steady cases coupled with nonlinear static structural analysis. Surface spline interpolation is responsible for the information exchanges between the structural deformation and aerodynamics. Thus, the lifting surface can be updated according to the actual structural deformation automatically in each iteration. Meanwhile, the lifting surface mesh and the normal vector in (13) are also updated to satisfy the exact boundary condition. In this way, the geometric nonlinearity is fully considered and the final static equilibrium configuration is obtained by the coupling iterative computation. It is true that this procedure can result in a nonconservation work due to nonconservative aerodynamic loads. Meanwhile, the stability limits can be significantly reduced according to previous publications [20–22]. In Wall’s paper [20], it is stated that every loose coupling scheme inherently suffers from accuracy and stability problems. The statement is based on the character of the loosely coupled schemes in which the aerodynamic loads are solely computed on the predicted deflection state while the structural solver converges to the actual deflection state of the new time step. Therefore, Wall concluded that a sufficiently stable coupling can be obtained by introducing subiterations, in which the aerodynamic forces and the structural deformations exchange not only once but several times per time step. However, it is not the focus of this work to quantify the loss of accuracy. Besides, we found that results obtained by loosely coupling show good agreement with the experimental results [5, 23] referring to flutter and dynamic responses. Considering the purpose of this paper and the computation cost, the loose coupling is used in static aeroelastic analysis to obtain the equilibrium state.
Here is an example of nonlinear static aeroelastic analysis, as shown in Figure 4. The presented wing model will be also used in the following optimization analysis. The airflow speed is 37 m/s. This figure presents two parts: force interpolation and displacement interpolation. For each iterative step, updated aerodynamic loads are applied to the previous structure, and the displacements are then solved by ULF. It is easy to find that the variation of the structural deformation turns to be smaller and smaller after each iteration step, and finally the iteration stops when the static equilibrium is obtained.
On the basis of the static equilibrium configuration, the vibration modes and frequencies of dynamics system equations can be calculated using quasimodal method, which assumes that the structure vibrates slightly around the nonlinear equilibrium state. Therefore, the stability problem at this state can be solved by traditional flutter analysis method. By selecting several main structural modes, the general aeroelastic equation in the frequency domain can be written asin which is the general inertial matrix and is the general stiffness matrix, which is the tangent stiffness based on the static equilibrium state. is the general aerodynamic coefficient matrix, which is a complex function of the reduced frequency. The aerodynamic force here is calculated using the unsteady DLM, considering the nonplanar effect. Equation (21) can be solved by - method, the same way as in the traditional flutter analysis. in which is the complex eigenvalue of the system at a given speed and and are the imaginary and real parts of . The structural damping coefficient is , and the frequency is . The critical flutter speed can be determined by judging the - curve (where goes through negative to positive values).
However, the flutter speed obtained in this way is not the exact solution because the flutter analysis is related to the static equilibrium, which is determined by the given flight conditions. It means that the damping value and the frequency in - curve and - curve are accurate only when the given speed equals the flight speed. Thus, in order to get the exact flutter speed, the iterative processes should be repeated on a series of flight speeds. Then the accurate - and - curves can be obtained by plotting the damping values and the frequencies with respect to the corresponding flight speed. The exact flutter speed can be identified through the - curve, and the flutter frequency is obtained from - curve. As shown in Figure 3, to solve the exact flutter speed, two iteration cycles exist in the whole process. In this study, the flutter speed will be one of the constraints in optimization problem. However, the key point of nonlinear optimization is not to obtain the exact flutter speed but to evaluate the dynamic stability at the flutter constraint. In order to save computing time in optimization process, the flutter constraint is set as the input of the nonlinear aeroelastic analysis and then the calculated structural damping value is judged whether it is greater than zero.
2.4. Optimization Problem
A wing model designed for wind tunnel test is optimized in the present work. Keeping the material properties and the wing’s span unchanged, design variables of the optimization problem are where are parameters of the th cross section of the wing beam as shown in Figure 6. The wing masses are determined by these parameters. The analysis can be expressed as a constrained optimization problem:
The optimum solution has to satisfy several constraints. In order to study the nonlinear aeroelastic characters, the wing model is required to produce large deformation. Thus, the designed flutter speed of the wing model cannot be too low to meet the referred requirement or too high in consideration of the difficulties of protecting the wing model in a wind tunnel. The constraint of flutter speed can be described as
This makes a big difference from the optimization design of a wing model for the flight test, in which the flutter speed is always constrained to be higher than the flight speed to ensure a safe flight.
In the same way, the constraint of the wingtip displacement in the trimmed wing configuration is considered:
Furthermore, some other variables should also be constrained within certain limits, includingwhere is the torsion angle at the wing tip and is the maximum stress in the main beam, which is required to be no more than the structural strength limit.
The optimum solutions are obtained by the ISIGHT (version 5.7) platform, which can be taken as a toolbox containing various optimization algorithms. In present work, the “Direct Search” method proposed by Hooke and Jeeves [24] is used for obtaining the optimum wing shape. This technique uses a combination of objective and constraints penalty as the objective function and it is well-suited for linear and nonlinear design spaces. Besides the optimization algorithm, an external program for nonlinear aeroelastic analysis has to been integrated into the optimization process to supply the inputs and outputs. The optimization process is shown in Figure 5. An input file containing initial values of design variables is given firstly, followed by setting the iterative parameters and initializing the system environment. Subsequently, the maximum stress, the wingtip displacement, and torsion angle are calculated by nonlinear static aeroelastic analysis. The flutter speed is determined by nonlinear flutter analysis. According to the above results, the design variables are updated automatically through the optimization algorithm. Finally, all constraints are satisfied and the optimum solution is obtained.
3. Description of the Wing Model
In this section, a wing model is designed for wind tunnel test. Unlike a real wing, a simplified FEM model without control surfaces was designed as shown in Figure 6. Its detailed parameters are listed in Table 1. The stiffness of the wing is supplied by the aluminum beam along 40% chord. The length of beam is 1542 mm. The crisscross section of beam is shown in Figure 6. This kind of section design has a good designability in both vertical and horizontal stiffness. The width and thickness of the beam section vary along the spanwise direction, and the initial design variables of all 11 segments are listed in Table 2. The total mass of the wing consists of two parts: the beam mass and lumped mass. The initial value of total mass is 1.983 Kg. The mass properties except for the beam mass were given in a lumped form. 22 lumped elements are set up along the leading and trailing edge and the mass values are listed in Table 3.
4. Numerical Results
The wing model is firstly optimized without considering the geometrical nonlinearity, followed by the nonlinear aeroelastic reanalysis of the optimized model. Subsequently, geometric nonlinearities are fully considered in optimization process, and some optimization strategies are studied to improve the computation efficiency.
4.1. Linear Optimization
In this section, the wing model is firstly optimized without considering the structural nonlinearity caused by large deformation. Therefore, traditional method for aeroelastic analysis is applicable for solving this optimization problem. Modal analysis results can be directly applied to the flutter analysis instead of the complex iterative processes for obtaining the equilibrium configuration. The trim analysis is conducted in a given inflow speed and an attack angle. Considering the real conditions of a wind tunnel test, the initial values of the inflow velocity and attack angle are set as 36 m/s and 4.0°, respectively. The whole optimization process was conducted on a computer with Inter Xeon CPU E3-1226 and 8.0 GB RAM.
Figure 7 shows the linear optimization process. The wing weight, wingtip displacement, wingtip torsion angle, and flutter speed are monitored during the iterative processes. It can be found that the Direct Search makes two types of moves—the procedure of going from one given point to the following point. The first is called exploratory move, which is designed to acquire knowledge concerning the searching direction. Most of the moves in Figure 7 are exploratory moves. The second type is pattern move, which changes suddenly, followed by a sequence of exploratory moves. After 983 iterations, the optimum solution was obtained, and all constraints are satisfied.
(a) Wing weight
(b) Wingtip displacement
(c) Flutter speed
(d) Torsion angle
In Table 4, the wing weight drops by 12.1% after optimization. Meanwhile, the bending and torsional stiffness are both reduced in accordance with the increase of the wingtip displacement and torsion angle. From the modal results listed in Table 5, one can find that the frequencies of the first 6 main modes are dropped by more than 30% after linear optimization. In other words, the wing becomes lighter and more flexible. The flutter speed after optimization becomes 38.4 m/s, dropping by 34%.
4.2. Nonlinear Reanalysis
In order to study the influence of geometric nonlinearity on aeroelastic characteristics, nonlinear aeroelastic analysis is conducted on the optimized wing model after linear optimization. The results are listed in Table 6. For the same wing model, there is a drastic change in structural and aeroelastic characteristics. It can be seen that wingtip displacement is beyond the upper boundary of the constraint, and the flutter speed is lower than the minimum of the constraint condition. To better understand this variation between linear and nonlinear analysis, the modal results are listed and compared in Table 7. It can be seen that the frequencies of the listed six modes all decreased in nonlinear reanalysis. The first vertical bending and first chordwise bending are influenced mostly by the geometrical nonlinearity, and the frequencies of these two modes are dropped by 25.6% and 27.4%, respectively. However, this kind of reduction is not the same as that in Table 5, which is mainly caused by the decrease of the total mass of the wing model. The aerodynamic loads and the associated deflection lead to a drastic change in the structural dynamic characteristics of the wing and consequently a catastrophic change in aeroelastic stability. According to the results in Table 7, the frequency differences tend to be smaller and smaller with the increase of the modal frequency. Combined with the results in Table 6, the decrease in flutter speed can be attributed to the decrease in the low-order mode frequencies. There is about a 10% decrease in flutter speed which significantly compromises the flight boundary of the test wing model. The results of nonlinear reanalysis show that the optimum solutions of linear case may be inaccurate. Therefore, it is necessary to consider the structural nonlinearity in aeroelastic optimization design for suchlike flexible wing in present work.
4.3. Nonlinear Optimization
In this section, geometrical nonlinearities are considered. Initial values of the airflow speed and the attack angle are the same as those in the linear optimization. The trim analysis is not a direct one-step process, but complex step-by-step iterative processes. Meanwhile, the flutter analysis also needs lots of repeated iterative calculations. However, in this study, some optimization strategies are adopted in nonlinear aeroelastic analysis to simplify the iterative processes, which will be described detailed in next section.
Figure 8 shows the searching process. The optimum solutions are obtained after 531 iterations. From Figures 8(a) and 8(c), one can see that the changing trend of the flutter speed is the same as that of the beam weight. However, both of the wingtip displacement and torsion angle show a contrary tendency with the beam weight. This is because the reduction of the beam weight leads to the structural stiffness decreasing, which makes the structural deformation increase under the same aerodynamic loads. What is more, the mode frequencies also decrease with the reduction of the beam weight, meaning that the bending mode is more likely to be coupled with the torsion mode and causes the dynamic instability. Thus, the flutter speed decreases with the weight decreasing. It is easy to find that there are two pattern moves in the nonlinear optimization process. The first one is around 70 iterations. The beam weight decreases suddenly to about 1.86 Kg. Meanwhile, the flutter speed, the wingtip displacement, and torsion angle show sudden changes and all exceed the corresponding constraint conditions. The second pattern move occurred after about 110 iterations to make the constraints satisfied. Finally, the optimum solutions are obtained after 531 iterations.
(a) Wing weight
(b) Wingtip displacement
(c) Flutter speed
(d) Torsion angle
Nonlinear optimization results are listed in Table 8. After optimization design, the wing weight drops by 4% and the flutter speed becomes 36.17 m/s, which is very close to the lower limit of the flutter speed constraint condition. Similar to the result in the linear optimization, the displacement of the wingtip reaches the upper limit of the constraint condition. These mean that the wing turns to be more flexible after optimization. The results obtained by the nonlinear optimization and linear optimization are both satisfied the constraints. Compared with the linear results, it seems to be unexpected that the optimum wing weight obtained by nonlinear optimization is heavier, since the analysis is closer to the actual situation. From the previous discussion, one can find that the optimum solutions by the linear analysis are inappropriately viewed from the nonlinear reanalysis. When nonlinear effects are considered, the optimum solutions are safer and more accurate. And technically to say, the “Direct Search” method is a local optimization algorithm. Due to the complexity of the nonlinear aeroelastic optimization problem, the searching process is easy to fall into local optimal solutions. However, the aim is not to investigate the optimization algorithm for nonlinear aeroelastic optimization in this paper.
Table 9 shows the change of the frequencies of the first six main vibration modes. It can be seen that frequencies of these modes all decreased. The nonlinear flutter analysis shows that the most important modes contributing to the flutter are the first vertical bending mode, the first chordwise bending mode, and the first torsion mode, whose frequencies drop by 23.6%, 22.6%, and 9.6%, respectively. To further study the differences of the flutter results before optimization from those after optimization, - and - curves are plotted, respectively, in Figure 9. The results show that the first vertical bending mode goes though zero to positive values, and the coupling modes that result in dynamic instability keep unchanged after optimization. The large deformed wing mainly induces the stiffness coupling of vertical bending and torsion, and the frequencies of the two modes are changed obviously.
(a) - curve before optimization
(b) - curve before optimization
(c) - curve after optimization
(d) - curve after optimization
4.4. Optimization Strategy
Nonlinear aeroelastic analysis is quite complex. It is undoubtedly that the process would become more time-consuming and complicated combined with optimization analysis. In order to improve the efficiency of nonlinear aeroelastic optimization, it is necessary to do more research on optimization strategies.
Apparently, the most time-consuming part of the optimization is nonlinear aeroelastic analysis in that too many iterations are needed to obtain the solution. More specifically, the iterative processes start over again every time the variables are varied. For instance, it needs about 2,500 times of iteration assuming that the optimum results can be acquired after 500 steps and 5 iterations are required in each step. This undoubtedly needs lots of computing time. However, in static aeroelastic analysis, the deformation and the distribution of the aerodynamic force on the wing surface in the equilibrium state show little change when the variation of design variables is not big and the wind speed is unchanged. For this kind of situation, the previous equilibrium state can be directly introduced into the next calculation. This means that load calculated in last step can be directly used in the updated wing model without any iteration. This is an approximation method to obtain the nonlinear static aeroelastic results, which can greatly improve computation efficiency. The specific calculation processes can be described as follows:(1)Apply the load calculated in ()th step to the updated wing model and obtain new static deformation denoted as .(2)Denote as the deformation of last iteration step and as the criterion value indicating the difference between these two-step deformation results. If the following inequality holds, the aerodynamic force of the previous equilibrium state can be directly used for the next step of nonlinear static aeroelastic analysis and flutter analysis:Otherwise, restart the iterative calculation to obtain the equilibrium configuration and then calculate the flutter speed.
In this way, the optimization process can be modified as Figure 10 shows. Compared to the optimization process given in Section 2, there is an additional judgment before the static aeroelastic analysis. The flutter analysis will not start until the second condition proposed above is satisfied. In order to study the influence of the criterion value on nonlinear aeroelastic analysis, different values of are selected for optimization analysis, respectively, including the situation of which shows the results without simplified method. The comparison of computing time with different criterion value is shown in Figure 11. It is easy to find that the computing time is greatly reduced using the simplified method, and it takes minimum time when the criterion value is 3%. The maximum computing time occurs in the value of 1%. This is mainly because the criterion is so strict that the iterative processes cannot be omitted in most occasions. However, what is seemingly a paradox is that the calculation time becomes longer when the criterion value is higher than 3%. Actually, it is not difficult to understand that the results become inaccurate with the increase of the criterion value, which simultaneously brings bad influence on the optimization analysis and makes the computing time increase.
Besides the comparisons of computing time, the errors of beam weight, the flutter speed, and the maximum displacement are also considered to make comparisons for determining a proper criterion value. Taking the results obtained when is 0% as the standard criterion, the relative errors of the beam weight, the flutter speed, the maximum displacement, and stress of the beam root are shown in Figure 12. The results show that the error of mass and flutter speed reaches the minimum in the criterion value of 3%, and the relative errors are both no more than 5% in other two aspects. It can be seen that the error of root stress reaches up to 25% or even higher when the criterion value is greater than 3%. This is mainly because the results become unreliable and inaccurate with the increase of the criterion value. Thus, the criterion value cannot be too high in an optimization design.
(a) Error of the mass
(b) Error of the flutter speed
(c) Error of the maximum displacement
(d) Error of the maximum stress
Considering the computing efficiency and the accuracy of optimization results, 3% is taken as the criterion value in nonlinear static aeroelastic iterative processes in this article. In fact, the optimum results in last section are obtained by setting a criterion value as 3%.
Besides the static aeroelastic analysis, it also needs lots of computing time to obtain the accurate flutter speed. As discussed before, there are two methods of nonlinear flutter analysis. The first method is to calculate the damping and frequency according to a series of given inflow velocities and the second one obtains the flutter speed by several iterations. Both of the two methods inevitably need repeatedly nonlinear static aeroelastic iterations, which is a time-consuming process. However, in optimization analysis, the only focus is whether the flutter speed satisfies the constraint condition so that the precise result becomes unnecessary. Therefore, a simple approach is to set the boundary of the flutter constraint as the input of the nonlinear aeroelastic analysis and then judge whether the calculated structural damping value is greater than 0. A damping value greater than 0 means the constraint of flutter speed is satisfied. In this way, the computing time can be greatly reduced.
5. Conclusions
This paper focused on a metal beamlike high-aspect-ratio wing model and optimized the beam section subjected to multiple aeroelastic constraints including the wingtip vertical displacement, the wingtip torsion angle, and the flutter speed. Synthesizing the analysis in this paper, the following conclusions can be drawn.
The study established a theoretical framework of aeroelastic optimization design for high-aspect-ratio wing considering structural nonlinear effects. This study can be applicable for not only aeroelastic problems but also for other similar optimization problems especially for flexible structures. The theoretical formulation for aeroelastic analysis can also be applied to solving other fluid-structure coupling problems, such as coupling vibrations and related dynamic stabilities.
The results of nonlinear reanalysis show that the optimum solutions of linear case might be inaccurate when the wing produces large deformation, and it is necessary to consider the geometric nonlinearity in optimization design. The nonlinear optimum results highlight a number of important features of optimum wings design with constraints on flutter speed and wingtip deformation. Briefly, the optimization decreased the section parameters and made the beam more flexible. The wingtip displacement firstly reached the upper limit in the optimizing process, which means that it is the most important constraint for this problem.
To simplify the complex optimization process and reduce the computing time, 3% is selected as the criterion value in simplifying the nonlinear static aeroelastic analysis. In this way, the computing time is reduced by more than 30% on the premise of ensuring the accuracy.
Of course, this paper focuses on establishing a theoretical framework for optimization design of very flexible wings. There remain some works to be improved. One is to adopt a new global search method to solve the optimization problem instead of the Direct Search method to avoid locally optimal solutions. In addition, more research work should be done to verify the framework for complex full-aircraft model.
Conflicts of Interest
The authors declare that they have no conflicts of interest.