Full Text

Turn on search term navigation

Changchuan Xie 1 and Yi Liu 1 and Chao Yang 1 and J. E. Cooper 2

Academic Editor:Samuel da Silva

1, School of Aeronautics Science and Engineering, Beihang University, Beijing, China
2, Department of Aerospace Engineering, University of Bristol, Bristol, UK

Received 21 December 2015; Revised 22 March 2016; Accepted 3 April 2016

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Large-aspect-ratio wings produce a high lift-drag ratio and good flight performance for high-altitude long-endurance unmanned aerial vehicles (HALE UAVs). Advanced composite materials make wing structures lightweight but also introduce large elastic deformations during aerodynamic loads in flight, which induces significant geometric nonlinearity in aeroelasticity and flight dynamics. In particular, after the NASA Helios mishap, the shortcomings of traditional linear analysis were deemed unsuitable for VFAs with large deflections, and the following recommendation was made by the investigators of the accident [1]:

[There is a need to] develop more advanced, multidisciplinary (structures, aeroelastic, aerodynamics, atmospheric, materials, propulsion, controls, etc.) time domain analysis methods appropriate to highly flexible, morphing vehicles.

Structural geometric nonlinearity, a key feature of HALE UAVs, yield a nonlinear relationship between displacement and strain due to large elastic deformations. Because the structural stiffness and equilibrium equations all depend on instantaneous structural deflections and load condition, the aerodynamic calculations and dynamic equations should be established for a large deformed state [2]. Thus, geometrically nonlinear aeroelasticity could be defined as a subdiscipline of aeroelasticity that considers nonlinear large structural deformations and the aerodynamics on curved aerosurfaces simultaneously.

In traditional linear analysis, aircraft were not particularly flexible, and the geometric nonlinearity were not significant; thus, linear structural finite elements based on the small displacement assumption combined with the planar doublet-lattice method [3] were widely used in engineering analyses and even imbedded in commercial software [4]. As aircraft have become more flexible, researchers have found linear aeroelastic analysis methods to be inaccurate and many new analysis methods for VFAs have been brought up. Hodges developed a nonlinear geometric exact beam element [5], reducing the order of structural nonlinearity, to analyse the nonlinear aeroelastic response of VFAs. He reported that structural geometric nonlinearity had a significant effect on the structural dynamics and dynamic aeroelastic characteristics of a high aspect ratio wing. The geometrically exact calculation of the angle of attack and aerodynamically consistent application of the air loads was also important for accurate aeroelastic characterisation [6]. Combined with Peter's 2D inflow theory [7], an aeroelastic analysis toolbox named NATASHA was developed to analyse the nonlinear behaviour of VFAs. The nonlinear beam model and reduced order model (ROM) combined with 2D or quasi-3D aerodynamic theory help researchers to predict VFA performance [8-10] but prevent wide application in industry because a real structure is overly complex (i.e., it cannot be simply represented with a beam model); thus, it is necessary to find a 3D aerodynamic code to manage VFA aerodynamic computations. To obtain accurate aerodynamic loads, some researchers used CFD/CSD [11] methods to analyse the geometric nonlinearity of VFAs. Smith et al. [12] loosely coupled the Euler solver with geometric exact beam structural analysis and investigated the effects of adding aerodynamic nonlinearity to the elastic behaviour of high aspect ratio wings. For VFA stability analysis, considering the high computing costs of time domain aerodynamic computations and the importance of highlighting basic aeroelastic principles for unconventional wings with high aspect ratios [13], dynamic flexible motion of the system can be assumed with small amplitudes around a nonlinear static equilibrium state [14]; thus the linearized method and frequency domain solution are still a valuable approach in preliminary designs and even in the detailed design stage. Panel aerodynamic methods have been well understood and widely used in engineering design; extending the panel aerodynamic code into 3-Dimentional application can make geometrically nonlinear aeroelastic analysis easy to accept in industrial applications.

Based on the above discussion, in this paper, a practicable geometrically nonlinear analysis system is established and a wind tunnel test is conducted to validate the analysis results. Nonlinear finite element method (FEM) is utilized so that the method could be easily used in industry. Actually there are some refined modelling methods like intrinsic beam [5], strain-based beam and plate elements [15], nonlinear substructure, ROM, and so forth. These methods are well developed theoretically and can illustrate the flexible structural characteristics to different degrees. Considering the practical problems we have met, the structures are often so complex that these modelling methods are not convenient to apply. Additionally, these simplified methods cannot well reflect the original structure especially in detailed parts. Moreover, FEM is the most often used modelling method in practical analysis and it will be very efficient and convenient if the FEM model can be directly used to analyse the VFAs structural geometric nonlinearity. The purpose of this paper is to establish an analysis framework easy to implement and able to reveal some nonlinear aeroelasticity of flexible wings. Based on the concerns above, the nonlinear FEM is responsible for structural nonlinear analysis. As to the aerodynamic computation, the importance of the nonplanar aerosurface effect and exact aerodynamic modelling consistent with structural deflection comes from our experience and many reference papers. Structural nonlinearity and nonplanar aerodynamic compotation must be considered simultaneously. Considering the easy programming and good inheritance of conventionally used linear method, NVLM (nonplanar vortex lattice method) and NDLM (nonplanar doublet-lattice method) are adopted and they can well present the nonplanar effect of aerodynamic for flexible wings in stability analysis. All of these methods together can describe the geometric nonlinearity of both the aerodynamics and structure of VFAs and can be conveniently used in industrial design. Wind tunnel test under different load cases indicates that different deformations under varying flight states result in different flutter speeds that may alter the aircraft's flight envelope. Reasonable agreement between the analysis results and test results has been obtained, and all the results demonstrate that the static aeroelastic response and flutter characteristics are quite different from the results obtained from linear analysis.

2. Theory

2.1. Structure Geometric Nonlinearity

Large structural excursions induced by very flexible wings when undergoing aerodynamic loads in the air prevent the use of linear methods based on the small displacement assumption and call for geometrically nonlinear structural analysis. Geometric nonlinearity 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 motion. These result in the nonlinear geometric equations including the quadric term of the displacement differential and require the nonlinear force equilibrium equation established on the deformed state of the structure. Structural geometrically nonlinear problems are often solved by the maturely developed nonlinear incremental finite element method [16] and two formulae called the total Lagrange formulation (TLF) and the updated Lagrange formulation (ULF) [17], which are well known. The ULF is used in this study, and the primary equations are presented briefly below. The core method of structural analysis has already imbedded in commercial software, so it is convenient to be used in engineering analysis.

The relationship between the nonlinear Lagrange/Green strain and displacement is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the partial derivative of displacement component [figure omitted; refer to PDF] to the coordinate [figure omitted; refer to PDF] at time [figure omitted; refer to PDF] .

Despite a large elastic deformation, the material is still within the elastic limitation for a small strain, so the conjugate Kirchhoff stress tensor [figure omitted; refer to PDF] at time t satisfies [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the direction cosine of a small area element [figure omitted; refer to PDF] at time [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is the corresponding surface force in which the follower force effect is considered. The linear elastic constitutive relation is given as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the elastic tensor, which has a different form for isotropic or anisotropic material.

The finite element method (FEM) based on energy principles is an effective approach to solve structural problems. For geometric nonlinear problems, considering the follower force effect, the incremental FEM is used. The strain [figure omitted; refer to PDF] can be decomposed into a linear part [figure omitted; refer to PDF] and a nonlinear part [figure omitted; refer to PDF] of the current displacement: [figure omitted; refer to PDF] The stress can also be decomposed in increments, where [figure omitted; refer to PDF] represents the equilibrium stress at time [figure omitted; refer to PDF] and [figure omitted; refer to PDF] represents the incremental stress to be calculated at each time step: [figure omitted; refer to PDF] The integral equation is established by linearization in each incremental step: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the incremental outer force including the aerodynamic force, engine thrust, and gravity, at the new time step. Considering a number of shape functions, the relationship between strain and deformation is presented as [figure omitted; refer to PDF] Substituting these shape functions into (6) leads to the element governing equation [18]: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the incremental outer force including the aerodynamic force, engine thrust, and gravity at the new time step. The stiffness matrix in (8) could be decomposed into a linear part [figure omitted; refer to PDF] and nonlinear part [figure omitted; refer to PDF] . The linear part is only related to the structure itself, whereas the nonlinear part is related to the deflected configuration and strain quality, which should be updated in each computation step.

The corresponding dynamic equation can be expressed as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the structural displacement acceleration vector at new time step. The assumption of a small amplitude vibration around the static equilibrium state is suitable for many dynamic problems, including the flexible wing structural dynamic stability: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the large deflecting equilibrium deformation from (8) and [figure omitted; refer to PDF] is a small vibration deformation. According to (9) and the static equilibrium condition, the linearized structural quasimode can be obtained by generalized diagonalization, and the vibration equation of the system under steady forces reduces to [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the inertial matrix of the structure at the nonlinear static equilibrium configuration and [figure omitted; refer to PDF] is the corresponding stiffness matrix. Both of these parameters are nonlinear functions of [figure omitted; refer to PDF] and vary under different equilibrium states, which is a key feature of geometric nonlinear structures. The mode shapes and frequencies can be deduced from (11).

Introducing the harmonic oscillating assumption [figure omitted; refer to PDF] , the vibration equation can be written as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the vibration circular frequency and [figure omitted; refer to PDF] is the vibration mode matrix. If (12) has all-nonzero solutions, that demands [figure omitted; refer to PDF] That is a generalized eigenvalue problem about [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Solving (13), the vibration circular frequency [figure omitted; refer to PDF] should be obtained. Substituting [figure omitted; refer to PDF] into (12), structural eigenvector [figure omitted; refer to PDF] (structural mode shape) can be obtained.

2.2. Nonplanar Vortex Lattice Method

The modelling of flexible aircraft with significant structural deformations requires the incorporation of structural dynamics and aerodynamics in a unified framework. The aerodynamic loads in (9) should be computed on a deformed configuration. This section summarizes the primary characteristics of the NVLM and its application to a curved aerosurface. The NVLM is derived from three-dimensional potential flow theory and is suitable for most normal situations that VAFs may encounter. Because the simple programming effort is required, NVLM is easy to combine with structural dynamic computations to obtain the response results for aeroelastic structures. Additionally, the exact boundary condition is satisfied on the real wing surface in the NVLM, which can have camber and various platform shapes. Thus, the NVLM is convenient for use with very flexible wings, whose aerodynamic surfaces are subjected to large structural deformations [19]. The NVLM is implemented using vortex ring quadrilateral elements to discretize the curved lifting surface along with the wing's deformation, as shown in Figure 1.

Figure 1: Nonplanar vortex lattice model for a thin lifting surface.

[figure omitted; refer to PDF]

A Cartesian coordinate system is shown in Figure 1, whose [figure omitted; refer to PDF] -plane represents the undeformed wing plane, and the [figure omitted; refer to PDF] -axis points in the flow direction when the angle of attack is 0. Both the structural displacement and aerodynamic computation will be computed and described in this unified coordinate system. The leading segment of the vortex ring is placed on the panel's quarter chord line; the aeroforce of the panels acts on the midpoint of the segment, which is represented by "[white circle]" in Figure 2, and the collocation point (represented by "×" in Figure 2) is at the centre of the three-quarter chord line, where the nonpenetration boundary condition will be implemented. Aerofoil warp can be represented by its middle camber surface, and some typical panel elements are shown in Figure 2.

Figure 2: Arrangement of the vortex ring elements.

[figure omitted; refer to PDF]

The velocity induced at an arbitrary point by a typical vortex ring can be calculated by applying the Biot-Savart law [20] to the ring's four segments, but attention should be paid to the rings located at the trailing edge where two semi-infinite trailing vortex lines that model the wake flow are shed along the free stream direction, as shown in Figure 3. Then, the induced velocities at all collocation points could be represented as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are the induced velocity components along the [figure omitted; refer to PDF] -axis, [figure omitted; refer to PDF] -axis, and [figure omitted; refer to PDF] -axis, respectively, and [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are their influence coefficient matrices related to the current large deformed configuration.

Figure 3: Typical vortex elements.

(a) Common vortex ring element

[figure omitted; refer to PDF]

(b) Vortex element at trailing edge

[figure omitted; refer to PDF]

The Neumann boundary condition is implemented to obtain the vortex distribution. The nonpenetration boundary condition is applied on the aerodynamic lattice of the current configuration, which implies that the method is also geometrically nonlinear and can account for large wing deformations. For the collocation point of the [figure omitted; refer to PDF] th panel element, [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the velocity of the free stream; [figure omitted; refer to PDF] is the induced velocity at the [figure omitted; refer to PDF] th collocation point by all vortex elements, [figure omitted; refer to PDF] ; and [figure omitted; refer to PDF] is the normal vector of the [figure omitted; refer to PDF] th panel element reflecting the local spatial deformation.

The circulation of the bound vortices is obtained by solving (15). Thus, the pressure distribution over the lifting surface can be obtained through the Kutta-Joukowski lift theorem [21]. The aerodynamic force that acts on the [figure omitted; refer to PDF] th panel element is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the total vortex strength at the [figure omitted; refer to PDF] th panel element's quarter chord line, [figure omitted; refer to PDF] has different values depending on whether the panel is the leading edge panel or not, and [figure omitted; refer to PDF] is the vector describing the magnitude and direction of the [figure omitted; refer to PDF] th panel element's quarter chord line. Consider [figure omitted; refer to PDF] when the panel is located on the leading edge of the wing and [figure omitted; refer to PDF] when it is not.

The curved surface distributed aerodynamic force obtained by the NVLM combined with the nonlinear FEM with the follower force effect can describe the large deformed static equilibrium state. Some extreme conditions, such as stall and hover, require further consideration by some possible modification or more accurate calculation method, such as CFD. Here, only a normal situation that NVLM can manage is considered, which is sufficient for normal stability analysis for flexible wings.

2.3. Nonplanar Doublet-Lattice Method

The calculation of unsteady loads plays an important role across much of the design and development of an aircraft and has significant impacts on structural design, aerodynamic characteristics, weight, flight control system design, control surface design, and performance [22]. In stability analysis, spatial aerodynamic modelling was found to perform well in situations dominated by small amplitude dynamics around large quasistatic wing deflections [23]; thus, the concept of frequency and modal analysis could also be inherited from vibration theory. Then, the frequency domain aerodynamic methods could also be suitable for VFAs in stability analysis, in which the planar DLM is often used in traditional aeroelastic analysis. In this section, the DLM code is extended into nonplanar cases to account for the 3D unsteady loads of large-aspect-ratio wings with large deflections and can be successively applied in engineering practice.

To meet the demand of nonplanar aerodynamic computations, mesh dividing should be determined on the deformed surface and updated along with the structure deflection, as shown in Figure 4. In addition to the spatial lattices, local coordinates should be established to reflect the exact nonplanar configuration of the wing. The nonplanar effect not only is reflected geometrically but also should be contained in the kernel function [figure omitted; refer to PDF] .

Figure 4: Typical nonplanar lattice on a curved lifting surface.

[figure omitted; refer to PDF]

The relationship between the pressure and normal wash distribution in unsteady potential subsonic flow can be written as [24] [figure omitted; refer to PDF]

Rodemich demonstrated that the kernel function can be written in the following form; a detailed derivation and explanation can be found in [25]: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are geometric relations that represent the local deformations. [figure omitted; refer to PDF] and [figure omitted; refer to PDF] could be evaluated by integrals or via other methods [24].

The critical problem of NDLM is the implementation of exact geometric boundary conditions. The local normal wash velocity, shown in Figure 5, is computed via (17); then, the boundary condition is expressed as [figure omitted; refer to PDF]

Figure 5: Typical doublet-lattice line and its normal vector.

[figure omitted; refer to PDF]

Due to large deformations, the wing cannot be treated as oscillating about the [figure omitted; refer to PDF] -plane; thus, the real boundary condition is related to local geometric nonlinearity. The more generalized boundary condition written based on modes is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the nondimensional induced velocity in the normal direction by the spatial lattice, [figure omitted; refer to PDF] is the reduced frequency, and [figure omitted; refer to PDF] is the modal vibration shape function of the wing structure in a local normal direction and varies under different equilibrium states; [figure omitted; refer to PDF] already includes the large static deformation, which makes it different from the traditional DLM despite its similar format. Because of the nonlinear geometric stiffness effects discussed before, modal shape and structure stiffness should be updated according to the different deformed configurations when computing unsteady aerodynamic forces.

Substituting (17) and (20) into (19), the boundary condition can be written as an aerodynamic influence coefficient (AIC) in matrix form, where [figure omitted; refer to PDF] is the instantaneously geometry-related AIC matrix: [figure omitted; refer to PDF]

Although the nonplanar DLM has been developed and used for many years and it could not consider an excessive number of unsteady effects, it is still capable of use in flexible wing stability analysis because the structure can be linearized, and the hypothesis of small amplitude dynamics around large quasistatic wing deflections can be made.

2.4. Surface Spline Interpolation

The integration of the flexible structure motion and the corresponding aerodynamic computation relies on the coupling between aerodynamics and structure. This is a critical problem in aeroelastic analysis, particularly for VFAs, whose large 3D elastic deformations make conventional 1D or 2D interpolations invalid. The generalized surface spline interpolation [26] is used in this study to couple the aerodynamics and structure. This method can clearly describe the interpolations in arbitrary space dimensions, which makes it applicable for VFA aerodynamic/structure coupling, whose structural configuration is typically considered to be embedded in a 3D space. A brief introduction is presented here.

Consider a given vector set [figure omitted; refer to PDF] in [figure omitted; refer to PDF] -dimension space and the corresponding image vectors [figure omitted; refer to PDF] in M -dimension space. For the [figure omitted; refer to PDF] th component of [figure omitted; refer to PDF] , it can be expressed as [figure omitted; refer to PDF] where the symbols in (22) are the same as those in [27] and the coefficient [figure omitted; refer to PDF] could be given different value with [figure omitted; refer to PDF] . This can be used as the interpolation between structural grids and aerodynamic grids. Considering [figure omitted; refer to PDF] given structural grids with coordinates [figure omitted; refer to PDF] and a corresponding deformation vector [figure omitted; refer to PDF] , the relationship between the coordinates and the deformation of the grids could be written in matrix form: [figure omitted; refer to PDF] and according to (22) we know that [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are constant matrices related to the coordinates [figure omitted; refer to PDF] and corresponding deformation [figure omitted; refer to PDF] . [figure omitted; refer to PDF] is the coefficient matrix of the surface spline fitting function. When [figure omitted; refer to PDF] is nonsingular, [figure omitted; refer to PDF] can be solved as [figure omitted; refer to PDF]

Then, the deformation vector [figure omitted; refer to PDF] of [figure omitted; refer to PDF] aerodynamic grids with the coordinates of [figure omitted; refer to PDF] can be interpolated as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the constant matrix related to the given coordinates of aerodynamic grids [figure omitted; refer to PDF] . Equation (25) can be transformed into [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the spline matrix for displacement interpolation between the aerodynamic grids and structural grids. According to this relationship, the aerodynamic grid and meshes can be updated automatically, ensuring that the real deformed aerosurface is considered in every aerodynamic computation.

The other transformation also occurs here between the aerodynamic and structural force systems to obtain an equivalent force on the structure for the structural deflection analysis. When the aerodynamic forces [figure omitted; refer to PDF] and their equivalent structure forces [figure omitted; refer to PDF] perform the same virtual work on their respective virtual deflections, the structural equivalence of the two force systems is guaranteed by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the arbitrary virtual deflections satisfying (26). The nonplanar aerodynamic force computed is [figure omitted; refer to PDF] . Therefore, [figure omitted; refer to PDF]

The bidirectional coupling between aerodynamics and the structure helps to complete the iterative process to determine the equilibrium states, which is fundamental in nonlinear analysis.

2.5. Flutter Analysis

Nonlinear structural stiffness [28] and nonplanar aerodynamics due to large deformations have dramatic effects on the flutter response, and, thus, the conventional linear approach of aeroelastic stability analysis is not applicable. To determine the stability properties of the structure, a modal analysis is conducted at the nonlinear equilibrium state, and the aeroelastic equation can be derived as [figure omitted; refer to PDF] Equation (29) has the same form as that of the linear case, but each item is obtained from the nonlinear analysis introduced previously. [figure omitted; refer to PDF] describes the generalized unsteady aerodynamics vector, and [figure omitted; refer to PDF] is the unsteady nonplanar aerodynamic influence coefficient complex matrix related to the reduced frequency [figure omitted; refer to PDF] and the deformed configuration given by the NDLM. Using the [figure omitted; refer to PDF] method, (29) can be rewritten as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the generalized aerodynamic stiffness matrix and aerodynamic damping matrix, respectively; [figure omitted; refer to PDF] is the complex eigenvalue; [figure omitted; refer to PDF] is the decaying rate; and [figure omitted; refer to PDF] is the structural damping ratio.

In detailed computations with a certain Mach number, air density ρ , and flying velocity [figure omitted; refer to PDF] , solving (30) yields a series of [figure omitted; refer to PDF] / [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Drawing the [figure omitted; refer to PDF] and [figure omitted; refer to PDF] curves of each mode, when [figure omitted; refer to PDF] , the corresponding flutter velocity and frequency can be obtained. However, in the nonlinear analysis, the flutter velocity obtained here is only a predicted value that may change under different static loads and deformation cases. Thus, the nonlinear static aeroelastic analysis and stability analysis should be considered jointly, not separately, as in the linear analysis. The exact flutter velocity can be searched iteratively to use the geometrically nonlinear aeroelastic analysis flowchart below.

The geometric nonlinear aeroelastic analysis flowchart for a very flexible wing, as shown in Figure 6, can be decomposed into two parts: a nonlinear static aeroelastic iterative analysis and a larger iterative stability analysis that searches for an exact flutter speed. The static aeroelastic analysis is implemented using an automatic aerosurface update via surface spline interpolation to acquire the accurate static aeroelastic results. The geometrically nonlinear static structural analysis is executed by the FEM, which could reduce the limitations of the simple beam model and can be easily used in engineering practice. Typically, after four to five times iterative calculation, the nonlinear static aeroelastic analysis can get a convergent result. Flutter analysis is implemented by coupling the NDLM computation and the quasimodes calculation of the wing at its nonlinear equilibrium state obtained by the static aeroelastic analysis. Once those two iterations reach the same flight speed, the aeroelastic flutter analysis is accomplished. Otherwise, the analysis changes the flight status, such as the Mach number or the dynamic pressure, and repeats the calculation. Commonly, dichotomy is used to search for the exact flutter speed. If the predicted flutter speed is higher than the current computed flight speed, then the intermediate speed can be selected as the next computed flight speed. Most often, after four times dichotomy computation, the exact flutter speed can be found.

Figure 6: Geometric nonlinear aeroelastic analysis flowchart.

[figure omitted; refer to PDF]

3. Analysis Model

3.1. Wind Tunnel Test Model

The geometrically nonlinear characteristics of VFAs have attracted people's attention for many years; however, there are still some basic principles that require a more in-depth description; this situation indicates the importance of wind tunnel tests. Thus, a wind tunnel test model was designed to validate the theoretical analysis established in this paper and to provide insight into VFA aeroelastic characteristics. Although the proposed theoretical analysis framework can be widely used for models of different complexities, a simple wing model here was used to perform the wind tunnel test and illustrate the nonlinearity of the aeroelasticity of very flexible wings. The tested semiwing, whose deformation is quite large during testing, is fixed at the root (right side of the [figure omitted; refer to PDF] lattices, as shown in Figure 8). The FEM model as shown in (Figure 7) has been modified by ground vibration tests, in which different weights were attached to test the static deformation and model frequency under horizontal suspended cases. The wing was also fixed upright to minimize the gravity effect and obtain its approximate natural frequency. The natural vibration information obtained by FEM after model calibration is shown in Table 1. Moreover, the GVT test results and nonlinear analysis results were also presented in the table. Although the wing had been fixed vertically, as shown in Figure 9, the nonlinear characteristics could still be observed. Since the gravity acted along spanwise, the stiffening effect can be noticed from the commonly higher frequency in the test results compared with linear natural frequency. The structural linearized frequency results shown in Table 1 were obtained via nonlinear static analysis and quasimode analysis under nonlinear equilibrium as introduced in Section 2.1. The acceptable small errors between test results and linearized analysis results indicate our structural modelling is credible and nonlinear analyses [figure omitted; refer to PDF] were basically accurate and available to be used in the subsequent analysis.

Table 1: Structure mode frequency.

Mode number	Mode name	Natural frequency (Hz)	Test frequency (Hz)	Linearized frequency (Hz)	The error between test and nonlinear analysis (%)
1	Vertical 1st bend	1.78	2.00	1.95	[figure omitted; refer to PDF]
2	Horizontal 1st bend	10.60	10.56	10.63	[figure omitted; refer to PDF]
3	1st twist	13.75	13.93	13.74	[figure omitted; refer to PDF]
4	Vertical 2nd bend	17.36	18.10	17.60	[figure omitted; refer to PDF]
5	Vertical 3rd bend	54.00	55.83	54.26	[figure omitted; refer to PDF]

Figure 7: Structure of the FEM model.

[figure omitted; refer to PDF]

Figure 8: Aerodynamic model.

[figure omitted; refer to PDF]

Figure 9: Ground vibration test (wing was fixed vertically).

[figure omitted; refer to PDF]

4. Aeroelastic Stability Analysis

4.1. Static Aeroelastic Analysis

For comparison purposes, the static aeroelasticity of the wing model is calculated by the traditional linear method (represented by the label "linear" in Figure 10) and the nonlinear method developed in this paper (represented by the label "nonlinear" in Figure 10). In the linear method, the assumption of infinitesimal deformation is adopted, and the aerodynamic loads are computed by the steady planar vortex lattice method, which do not consider the wing's deflection. Although the conventional planar vortex lattice method is not accurate and even erroneous in this case, it has been adopted for more than thirty years and still has been widely used in engineering analysis. However, that is our paper's aim to show the inaccuracy of planar aerodynamic computation and establish an easily acceptable nonplanar aerodynamic computation to deal with the nonplanar lifting surfaces. Figure 10 shows the vertical displacements of the wing tip that are obtained by these two methods under a 0.5° angle of attack at the root; the results of the wind tunnel test are also presented.

Figure 10: Vertical displacement at the wing tip versus airspeed.

[figure omitted; refer to PDF]

The test results in Figure 10 indicate that vertical displacement of the wing tip increases with increasing wind speed, as does the curve slope of the test result until the velocity reaches approximately 33.5 m/s. When the velocity is further increased, the rate of wing tip displacement increase is slowed down and converges to a limit value due to geometric stiffness effect. The results obtained by the nonlinear method are consistent with the test results despite some small differences that may be caused by systematic measurement errors. However, the wing tip deflection obtained by the linear method increases with wind speed and reaches approximately 456.9 mm at 34.0 m/s, which is far beyond reality and is no longer significant.

According to the test results, the vertical displacement at the wing tip is approximately 37.5% of the span under a 0.5° angle of attack when the wind speed is 34.0 m/s. The initial and final aerosurfaces in the nonlinear analysis are presented in Figure 11. It is a typical nonlinear case of the flexible wing with a large deformation, and the nonplanar effect of the lifting surface is quite significant. Detailed displacements in three directions along the span in this case are presented in Figure 13. Figure 12 shows the linear planar case and nonlinear nonplanar case in static analysis. In the traditional linear analysis, the aerosurface remains planar, and the aerodynamic force [figure omitted; refer to PDF] acts only in the z direction; considering the linear relation between displacement and force, deflections in the x- axisand y -axis are small, shown in Figure 13. However, in the nonlinear analysis, the aerosurface is automatically updated with structural deflections, and the follower force effect is included so that the displacements in the x- axis and y -axis can be considered. As shown in Figure 12, the aerodynamic force [figure omitted; refer to PDF] acts vertically to the aerosurface. Thus, besides the lift, a large component side force is significant in nonlinear cases, which may induce large deflections in the [figure omitted; refer to PDF] -axis ( [figure omitted; refer to PDF] in Figure 12) and reduce the effective lift. There is induced drag in [figure omitted; refer to PDF] -axis in both linear case and nonlinear case, but it is only quite small part of force in [figure omitted; refer to PDF] -axis. The deflection in [figure omitted; refer to PDF] -axis under nonlinear case is mainly caused by the follower force effect other than the induced drag force. Once there is a little torsion, there is a component force of [figure omitted; refer to PDF] in [figure omitted; refer to PDF] -axis that may cause deflections in [figure omitted; refer to PDF] -axis. In linear case, the aerosurface keeps planar and the aero force only acts in [figure omitted; refer to PDF] - axis, and the structural deflection in [figure omitted; refer to PDF] -axis is rather small.

Figure 11: Initial and final aerodynamic model of the flexible wing.

[figure omitted; refer to PDF]

Figure 12: Aerodynamic force in linear and nonlinear case.

[figure omitted; refer to PDF]

Figure 13: Displacement of the wing spar spanwise (0.5° AOA, 34 m/s).

(a) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

Additionally, in the linear analysis, the effective aerodynamic load reduction due to the shortening of aerosurface and structural stiffening effect due to elastic deflections cannot be considered; as a result, the displacements along the [figure omitted; refer to PDF] -axis (Figure 13) and the torsion angle (Figure 14) of the nonlinear analysis are smaller and more reliable than the linear results. Because the wings are fixed at their roots, the torsion angle increases from zero along the span. The tip torsion angle is still within the limitation of potential flow, so the NVLM still can be applied.

Figure 14: Torsional angle of the wing spar spanwise (0.5° AOA, 34 m/s).

[figure omitted; refer to PDF]

The aerodynamic loads along spanwise, including the lift, drag, and side force, are shown in Figure 15. These loads are the directional component forces of the total aerodynamic load obtained from the NVLM. Since the aerosurface keeps planar, the effective lifting surface area is much bigger than the nonlinear case. Additionally, the torsion angle is bigger than nonlinear case, so the lift in linear case is bigger than that in nonlinear case. Because of the follower force effect and the updated aerosurface, aerodynamic force may have significant component force in [figure omitted; refer to PDF] -axis, so the drag in nonlinear case is bigger than that in linear case. Since the aerosurface keeps planar, there is no side force in linear case. In nonlinear case, the aerodynamic force acts vertically to the large deformed aerosurface so there is quite large component force acting along the span as the side force. The side force is comparable to the lift and contributes significantly to the root bend moment; thus, it cannot be ignored in the nonlinear analysis. The elastic torsional angle increases the effective angle of attack; thus, the maximum aerodynamic loads along the span move to the outer wing.

Figure 15: Drag, side force, and lift along spanwise.

(a) Lift

[figure omitted; refer to PDF]

(b) Drag

[figure omitted; refer to PDF]

4.2. Flutter Analysis

The structural dynamics are different under various deflections and force cases due to geometrically nonlinear effects. Figure 16 shows the structural dynamic frequency under different static load cases. The reduction of the horizontal bend frequency and the increase in the twist frequency are distinct. An iterative process is used to obtain the exact flutter boundary under certain flight conditions, and the effect of gravity is also considered. The results of different cases are compared with the experiment results to illustrate the nonlinear flutter characteristics and validate the theoretical analysis.

Figure 16: Structural frequency under nonlinear static deflections.

[figure omitted; refer to PDF]

Figure 17 shows the [figure omitted; refer to PDF] and [figure omitted; refer to PDF] curves of the traditional linear results. A classical bend/twist coupling flutter pattern results in a flutter speed at approximately 33.5 m/s, and the horizontal bend mode does not contribute to the flutter; its frequency remains unchanged. The 1st vertical mode frequency reduces to zero when the speed is over 34 m/s but the double calculation and the wind tunnel test indicate that there is no divergence occurring around 34 m/s.

Figure 17: [figure omitted; refer to PDF] and [figure omitted; refer to PDF] curves of the linear flutter analysis.

[figure omitted; refer to PDF]

The nonlinear stability analysis results are presented below. All [figure omitted; refer to PDF] curves are found to be significantly different from the linear results with an unstable mode switching to a horizontal bend in the nonlinear analysis. This is caused by the geometrically nonlinear effect, which makes the horizontal bending mode contain large torsional components and becomes unstable when the structure experiences a large deformation, as shown in Figure 18. It is a remarkable characteristic of nonlinear stability analysis of VFAs that the horizontal bend mode plays an important role in flutter characteristics.

Figure 18: Horizontal bending mode in a nonlinear equilibrium state.

[figure omitted; refer to PDF]

Figure 19 shows the stability analysis results under the equilibrium states of 20 and 25 m/s. The corresponding predicted flutter speeds are approximately 31 and 33 m/s. Typically, the flutter speed decreases with increasing airspeed in the nonlinear analysis, but the test wing induces a negative deformation when the airspeed is below 28 m/s. Thus, the predicted flutter speed of 25 m/s is larger than that of 20 m/s. Because of the disagreement between the static equilibrium flight speed and predicted flutter speed, only the data near the static state marked red in Figure 19 are reliable. An iterative process is performed to narrow the disagreement and search for the exact flutter speed. Figure 20 shows the [figure omitted; refer to PDF] and [figure omitted; refer to PDF] curves from the dynamic aeroelastic analysis under the static equilibrium flight speed of 31.5 m/s. The predicted critical speed is 31.5 m/s, which is the theoretical exact nonlinear flutter speed of the analysis model under 0.5° angle of attack. Again, the data in Figure 20 is only accurate near 31.5 m/s, and if completely accurate [figure omitted; refer to PDF] and [figure omitted; refer to PDF] curves were desired, the divisional accurate data should be pieced together, as is in Figure 21, which requires more data to make the curve smoother and more reliable.

Figure 19: Predicted flutter speed under two flight cases.

[figure omitted; refer to PDF]

Figure 20: [figure omitted; refer to PDF] and [figure omitted; refer to PDF] curves under 31.5 m/s.

[figure omitted; refer to PDF]

Figure 21: Accurate data pieced together roughly.

[figure omitted; refer to PDF]

The calculated cases above are under a 0.5° angle of attack. Due to geometric nonlinearity, cases under different angles of attack lead to different flutter speeds. The iterative processes searching for exact flutter boundary under different angles of attack have been implemented and the final results are shown in Figure 22. The flutter results are summarized in Table 2. The aerodynamic loads and structural deflections increase as the angle of attack increases. Due to geometric nonlinearity, large deflections and load conditions easily cause the horizontal bending to become unstable. Thus, both the test and nonlinear analysis results indicate that the flutter speed decreases with an increasing angle of attack, whereas no changes are noted in the linear analysis under different cases. The first three cases present nearly identical deviations between the nonlinear analysis results and test results (in Figure 23). This effect may be caused by the conservative aerodynamic computation that can be validated in the deflection curve in Figure 10, where the analysed displacement is larger than that found in the test when the displacement is positive. Thus, it is understandable that the analysed flutter speed is smaller than the test results. The deviation narrows in the last case because the test result changes to a stall flutter depending on the flutter phenomena observed in the test; the stall model is not considered in theoretical aerodynamic computation, but the decline in flutter speed can still be reflected. In our test cases, the linear analysis results sometimes are less conservative than nonlinear results. This is not a referable conclusion. For most cases, the nonlinear analysis results are more accurate and reliable than linear results, especially when the nonlinear flutter speed is lower than the linear results.

Table 2: Flutter speeds under different angles of attack.

Angle of attack	0°	0.5°	1.0°	1.5°	Linear
Computation	Experiment	Computation	Experiment	Computation	Experiment	Computation	Experiment	Computation
Flutter speed (m/s)	33.0	35.7	31.5	34.4	28.0	32.2	26.5	27.0	33.5
Flutter frequency (Hz)	7.60	6.94	7.82	6.88	7.34	7.00	7.53	13.38	7.29
Unstable mode	First horizontal bend in plane	First twist

Figure 22: Nonlinear analysed flutter speed under different angles of attack.

[figure omitted; refer to PDF]

Figure 23: Flutter limit speed at different angles of attack.

[figure omitted; refer to PDF]

5. Conclusion

A theoretical analysis framework has been introduced for very flexible wing structures that present notable geometric nonlinearity. The NVLM and NDLM are combined, and the proposed code has been developed to obtain the steady and unsteady aerodynamic loads for aeroelastic stability analysis. Wind tunnel testing is used to demonstrate the geometrically nonlinear characteristics of a large-aspect-ratio wing and validate the theoretical analysis results. Reasonable agreement between the computational results and test results has been obtained; all the results indicated that the static aeroelastic responses and flutter characteristics of very flexible wings are significantly different compared with the traditional linear analysis results. The geometrically nonlinear aeroelastic stability is related to certain flight conditions, and large elastic deformations make the horizontal bend mode unstable and may decrease the flutter speed, which have a significant influence on the flight envelope. The NVLM/NDLM and nonlinear FEM with the quasimodal dynamic approach make the nonlinear aeroelastic stability analysis feasible and more accurate. The theoretical analysis process established in this paper follows conventional linear analysis ways with some significant modifications considering VFAs' geometrical nonlinearity both theoretically and practically and can thus be easily used to estimate the occurrence of critical aeroelastic behaviours of VFAs in engineering analysis. Efforts to apply the methodology to a flexible aircraft including control surfaces and coupled flight dynamics are ongoing; some aerodynamic derivatives may also be discussed in future work.

Nomenclature