Nolinear Static Aeroelastic Analysis and Optimization for High-Altitude Solar-Powered UAV With Large Aspect Ratio

According to the static aeroelastic characteristics of Solar-Powered UAV (SP-UAV) with large aspect ratio, this paper develops a nolinear aeroelastic analysis method which takes into account both the calculation accuracy and efficiency. Firstly, we establish the finite element model of the geometric nonlinear structure, which is verified by the static test of the whole UAV to reduce the calculation error caused by the simplified method of structure model. Then, the aerodynamic model of SP-UAV is established by using Curved Vortex Lattice (CVL), and the structure-aerodynamic coupling interpolation is realized by the Thin Plate Spline-interpolation (TPS) method, so as to complete the full static aeroelastic analysis. The analysis indicate that the calculated results of the longitudinal trim are well match the first flight condition. In particular, considering the elevator efficiency loss occurred in longitudinal trim analysis of the UAV, we design an optimization scheme to improve the efficiency of the elevator by about 15%. The simulation results show that the loss of elevator efficiency of flexible solar UAV with high aspect ratio is non-negligible, with the maximum value higher than 40% under high dynamic pressure. In actual flight, the speed of UAV must be reduced as much as possible to minimize the loss of rudder efficiency, while the elastic control method should be considered in flight control.


I. INTRODUCTION
The area located 20-100km high on the earth surface is the near space, which has extremely important application and development value and has been widely concerned by scientific research. In near space, the Solar-Powered UAV (SP-UAV) platforms can use solar energy as power to achieve high altitude flight, long time duration and wide-area operations. Additionally, the SP-UAV is highly survivable due to its low exposure to weather and up-down convection, making the platform particularly suitable for emergency communications, space-based early warning, communications relays, reconnaissance and surveillance, forest fire prevention, and The associate editor coordinating the review of this manuscript and approving it for publication was Ming Xu . disaster monitoring [1], [2]. Compared to the conventionally powered UAVs, SP-UAVs do not need any conventional expendable fuel, and are cheaper to manufacture, operate and maintain, thus having a higher mission-to-cost ratio. Compared with the orbital satellites or aerosols, SP-UAVs have lower launch costs and technical risks. Besides performing diversified tasks, it can be quickly deployed, has flexible maneuverability, stronger controllability and the ability to execute tasks at fixed points, thus it can respond more quickly to the mission requirements in the target area [3].
The near-space SP-UAVs generally adopt large-aspect ratio wing to obtain good lift-drag characteristics, and a large number of composite materials are used to reduce the weight of the structure to meet the mission requirements of long endurance [4]. Due to its light weight and high flexibility, this kind of aircraft structure will produce large elastic deformation under the action of aerodynamic loads, and the deformation of wing-tip will even reach up to more than 15% of the wing span. However, the local deformation of the wing structure still satisfies the linear displacement and linear strain relationship [5], which is a typical geometric nonlinear problem with large displacement and small deformation. Such large deformation usually causes changes in the distribution and action direction of aerodynamic loads, leading to significant differences in the balance position of the structure before and after loading, and its aeroelastic problems become very prominent [6].
The static aeroelastic analysis of high-aspect ratio SP-UAV has been extensively studied by scholars. Patil and Hodges [5], [6], [7] firstly used the eigenbeam model to introduce geometric nonlinearity into the aeroelastic analysis of flexible wings with high-aspect ratio. Wang et al. in [8] conducted nonlinear aeroelastic analysis of high-aspect ratio wings by using nonlinear beam theory and optimized unsteady length vortex aerodynamic model. However, due to the complexity of the real engineering structure, it is difficult to simplify the beam model. Xie et al. [9], [10] also carried out in-depth research on this problem by using the linearized Updated Lagrange (UL) finite element method, but the coordinate transformation matrices at both ends of the element in the UL method were not consistent, and the influence of the bending deformation of the structure could not be analyzed. Carnie and Qin [11], [12] used FLUENT and ANSYS interface program to carry out the aeroelastic analysis of high aspect ratio wing, but there were some problems such as long calculation time. Zhang and Xiang et al. [13], [14] used aerodynamic ONERA method and completely nonlinear beam model to conduct in-depth analysis and research on the coupling problem of nonlinear aeroelasticity and flight dynamics of flexible aircraft with high aspect ratio. Ma and Zhang et al. [15] applied FLUENT/NASTRAN weak coupling solution to carry out the static aeroelastic analysis of high-aspect ratio wing, but the coupling analysis method has the problem that the calculation time is too long due to repeated iterations.
In view of the above problems of flexible UAV with high aspect ratio, the actual calculation error can't be avoid due to the simplified structure model is different from engineering structure model. This paper applies the analysis base on the finite element model of the whole structure verified by static test, and the influence of large structural deformation caused by aerodynamic loads is considered, too. The aerodynamic model is established by dividing the deformed aerodynamic surface into aerodynamic mesh, the aerodynamic mesh nodes are updated and iterated with the deformation of the structure through reasonable spline interpolation, and finally it reaches the convergence state to complete the static aeroelastic analysis of the UAV.
Tackle the above challenges, the contributions of this paper are summarized as follows: a) 1) The finite element model of the structure modified by static test is adopted, which is more accurate and can better reflect the nonlinear characteristics of the structure with large deformation; 2) The Thin Plate Spline-interpolation (TPS) method is applied to complete the coupled interpolation of aerostructure, and the influence of structural elasticity on the angle of attack, angle of elevator deflection and distribution of aerodynamic loads of UAV are analyzed. 3) To solve the problem of low efficiency of the elevator, we optimize the mounting angle of horizontal stabilizer to reduce the aerodynamic loads of the elevator surface, thus improving the efficiency of the elevator. The rest of this paper is arranged as follows: Section II describes the structural modeling and analytical design methods, as well as the design of analysis and calculation process. Then the proposed finite element modeling and validation is expound in Section III. According to the solution scheme, the static aeroelastic analysis and optimization scheme, as well as simulation results are presented in Section IV. The conclusions are summarized in Section V, finally.

II. STRUCTURAL MODELING AND AERODYNAMIC ANALYSIS DESIGN
In this chapter, the structural modeling and analytical design methods, as well as the design of nolinear static aeroelastic analysis and calculation process are created.

A. STRUCTURAL GEOMETRIC NONLINEARITY
The flexible wing UAV with high aspect ratio will have obvious bending and torsional deformation under external loads, so the linear solution method based on the traditional small deformation hypothesis is no longer applicable. However, due to the small structural strain, its constitutive relation is still linear, that is geometrically nonlinear problem, approximately [16], [17]. Geometrically nonlinear aeroelasticity theory considers the aeroelasticity problem of nonlinear factors from both structural mechanics and aerodynamics [18], [19]. This theory consists of three parts: geometric nonlinear theory of structure, the aerodynamic analysis method of curved surface after structural deformation and the corresponding coupling method of structure and aerodynamic interface.
Considering the engineering applicability of the structural model, the wing structure adopts a more general nonlinear finite element model. The incremental analysis method is usually used in geometric nonlinear finite element analysis. According to the different initial position referenced by displacement variable, it can be divided into Full Lagrang (FL) method and Updated Lagrang (UL) method. The equivalent virtual displacement principle of FL at time t + t is: where t s ij is the equilibrium stress at time t, t+ t Q is the external load virtual work at time t + t. The total strain is ε 0 ij = e 0 ij + η 0 ij , in which e 0 ij is the linear strain increment at the current time, and η 0 ij is the nonlinear strain increment. The virtual displacement principle of the UL can be achieved when all the variables take time t as the reference configuration: After the linearization of the above two equations, the virtual displacement principle of TL scheme can be obtained as follows: The virtual displacement principle of UL format is: (4) The above two Equations can be variated to obtain the linear equations of displacement increment. Thereinto, 0 D ijkl and t D ijkl are functions of time t, and refer to the tangent constitutive tensor of the configuration measures of time 0 and time t, and t τ ij represents the stress increment at time t.
By discretizing the structure with finite element method, the strain of the structure can be decomposed into linear part and nonlinear part. Then, the transformation relationship between structural strain and displacement increment is established through finite displacement functions, and the matrix equation of structural element can be obtained as: where K L is the linear stiffness term, K NL is the nonlinear stiffness term,Q is the external load of the new increment step, and F is the equivalent internal force of the structure.

B. CURVED SURFACE VORTEX LATTICE METHOD
As the SP-UAVs usually are large flexible aircraft with highaspect ratio, the large deformation of the structure will cause the aerodynamic loads to produce a lateral component and lead to the decrease of the effective lift area of the wing. Therefore, Vortex Lattice Method (VLM) is used to calculate the aerodynamic loads on the lifting surface [20]. A series of circular vortices are arranged on the middle arc of the wing to simulate the wing, and then the linear equations of vortex ring strength i can be obtained by using the boundary conditions of the impenetrable object surface of the wing surface: where Aij is the induced velocity of the ith ring vortex to the jth ring vortex under unit strength, then we have RHS k = U V W ∞ · n k , where n k is the outer normal vector of the kth plane element.
Based on i , the load on the wall vortex line can be obtained by using Joukowsky theorem. Then, the total lift, drag and torque characteristics of the wing can be obtained by integrating each circular vortex [21], [22].

C. AERO-STRUCTURAL COUPLED INTERPOLATION
As the aerodynamic mesh and nodes of flexible aircraft surface will form into spatial surface with structural deformation, the traditional two-dimensional wireless flat spline interpolation method IPS (Infinite Plate Splines) is no longer applicable. Thin Plate Spline (TPS) interpolation method is adopted in this paper, which has one dimension more than the IPS method. The surface described by simulation function w i (x i , y i , z i ) is used to replace the known data, and the equations required to be solved are: where it has For the displacement difference, to solve the above equation, four additional force and momentum equations are needed, those are: Then the displacement difference function δq a = Gδq s can be obtained, where q a is the aerodynamic displacement vector, q s is the structural node displacement vector, and the displacement interpolation matrix G can be obtained from the structure to the aerodynamic node. The interpolation form of the load can be obtained from the virtual work equivalence: As the wing and tail of the UAV studied in this paper contain multiple control surfaces, and the stiffness of the  control rudder differs greatly from that of the main wing structure, it is easy to produce scissors difference under the action of aerodynamic force. Therefore, the control surface and the wing are divided in aerodynamic interpolation and the difference is processed respectively, as shown in figure 1.

D. ANALYSIS AND CALCULATION PROCESS DESIGN
According to the above structural geometric nonlinear theory, curved vortex lattice theory and aero-structural coupling interpolation method, the aeroelastic analysis process (as shown in figure 2) can be established to conduct nonlinear static aeroelastic analysis for the structural geometric large deformation wing. Firstly, the steady aerodynamic loads of the structure after large deformation is interpolated to the structure, and the structural nonlinear solver is used to calculate the deformation of the updated structure. Then iterated repeatedly until the result of structural deformation satisfying the requirement of precision is obtained.

A. CALCULATION MODEL
The SP-UAV studied in this paper adopts conventional layout, single beam and multi-rib structure of carbon fiber composite material. The UAV has a wingspan of 35m, a length of 13.5m, a take-off weight of 308.5kg and a aspect ratio of 15.6. The engine compartment skin, wing ribs, frame and truss are simulated by SHELL elements, and truss bar on wing is simulated by CBAR elements. Similarly, the main beam of wing, main beam of horizontal tail and vertical tail, control surface rotating shaft and landing gear support rod are simulated by SHELL elements.
The finite element model of the whole UAV is shown in figure 3. Explicitly, main beam of the wing and other structures are designed with high modulus and high strength unidirectional tape materials for layering. Considering the huge difference of tensile and compressive strength of materials, the upper and lower surfaces of the wing beam are designed with variable thickness layering, and the layer dropping design is adopted in the length direction to reduce the weight of the structure as much as possible. The relevant material parameters are presented in table 1.

B. VALIDATION OF STRUCTURAL FINITE ELEMENT MODEL
The accuracy of the structural finite element model needs to be verified before the static aeroelastic analysis. Static tests are performed on the UAV structure prior to the first flight,   It can be concluded from figures 7 and 8 that the displacement error between the static test results and the simulation results is within 3%, indicating that the structural finite element model has high accuracy and can be used for the subsequent static aeroelastic analysis.

IV. STATIC AEROELASTIC ANALYSIS AND OPTIMIZATION SCHEME
This section presents a series of aerodynamic analysis and optimization scheme design, including static aeroelastic analysis, verification of static aeroelastic results, as well as optimization and improvement of tail wing.

A. STATIC AEROELASTIC ANALYSIS
According to the geometrical nonlinear aeroelastic calculation flow, the longitudinal static aeroelastic trimming problems of SP-UAV is analyzed. The calculation condition is set as: the altitude is 1.5km, true flight airspeed is 11m/s, takeoff weight is 308.5kg, and normal flight condition is 1.0G level flight. The static aeroelastic analysis results of the whole    The results show that the Angle of attack for elastic trimming is larger than that of rigid due to the elastic effect. Meanwhile, because of the poor elevator stiffness, the loss of aerodynamic efficiency is serious, resulting in a sharp increase in the trim angle of elevator compared with the rigid state. By extracting the aerodynamic loads on the wing, we can draw the shear force and bending moment of the wing along the span(as shown in figures 11 and 12). Additionally, the variation trend of torsion angle with wingspan of the wing and tail is shown in figure 13. The results show that the aerodynamic loads on the wing will be redistributed due to the elastic deformation of the structure.  The shear force of the wing in the elastic trimming state is lower than that in the rigid trimming state, but the shear force in the elastic state and the rigid state are basically the same at the symmetric plane of the wing. Secondly, the bending moment of the corresponding elastic state is smaller than that of the rigid state, and the bending moment of the elastic state is 6.1% less than that of the rigid state at the plane of symmetry of the wing. As the wing is a single-beam multi-rib structure, its torsional stiffness is poor, and the torsion angle of the wing tip of the elastic trimming state is 2.76 • which makes the wing bow, as shown in figure 13. Therefore the angle of attack required by the elastic trimming is larger.

B. VERIFICATION OF STATIC AEROELASTIC ANALYSIS
During the first flight, the rear camera is installed on the UAV, and the installation position of the camera has enough stiffness to keep the relative position with the body unchanged during the flight. The tail is located in the center of the field of view of the camera and is less affected by the distortion of the airborne camera, so it can be used as the basis for the measurement of the deformation of the tail. A stable flight segment is selected from the video obtained by the rearview camera and captured frame by frame. An image analysis script is used to identify the position of the vertical tail tip  in the image. By setting an appropriate recognition threshold, the error can be controlled within 1 pixel.
Taking the upper left corner as the origin of coordinates, the axis is positive in the right direction and positive in the downward direction, as shown in figure 15, then it can obtained: where y 1 and y 0 are the pixel coordinates of the displacement reference points in the parking state and the flying state, respectively, y is the length of the vertical tail pixel coordinate, and L 0 is the actual length of the vertical tail. Then, the displacement of the tail relative to the parking state can be calculated according to the pixel coordinates of the reference point at the top of the vertical tail at different times. When the UAV is parking on the ground, the displacement of the tail is L grav = 95mm under its own weight. In order to reduce the error, the average pixel coordinates of the tail reference points within a period of time are taken as the pixel coordinates of the reference points in the stable level flight state, and the average pixel coordinates of the test points are given as y 1 = 330.4, y 0 = 303. Then, the displacement of the tail fin relative to the parking state can be obtained as L = 256.9mm. Therefore, the displacement of the reference point of the tail fin during flight is L total = L grav + L = 351.9mm, where the image recognition error is 1 pixel, that is 1/(y 1 − y 0 ) = 3.65%.
By extracting the trimming load of the tail from the trimming result and substituting it into the load-displacement  curve of the static test of the fuselage and tail, the displacement of the tail under the trimming load can be obtained. table 3 shows the comparison between the static aeroelastic analysis results and static test results of tail displacement with the first flight results which is obtained by image recognition analysis. table 4 shows the comparison of simulation results of UAV's trimming angle of attack and trimming elevator deflection angle with the first flight. The calculation results show that the deviation between the simulation results of static aeroelastic analysis and the image recognition results of tail deformation in the first flight is 3.3%, and the error of the trimming angle of attack and the trimming elevator deflection angle is less than 5%. Therefore, the static aeroelastic analysis model we established has high accuracy and can be used for analysis and calculation.

C. TAIL WING OPTIMIZATION SCHEME
During the level flight of UAV, the elevator deflects upward (the trim rudder is negative), which leads to a larger angle of attack on the stabilizer surface of the horizontal tail, resulting in a torsion angle of head up, which further reducing the efficiency of the elevator surface. It is found that the installation mode of the horizontal stabilizer is not conducive to the transmission of torque, and its torsional stiffness is poor under the action of aerodynamic loads. In level flight, the maximum torsional angle of the stabilizer reaches 10.7 • and the loss of elevator efficiency even reaches 75.34% under high dynamic pressure. In order to improve the efficiency of elevator, it is necessary to optimize the connection mode of tail fins to improve its torsional stiffness, so as to improve the operating efficiency of elevator surface. Due to the large deflection angle of the elevator in trimming state, it is necessary to consider optimizing the mounting angle of the horizontal stabilizer to adjust the trimming pressure of the elevator surface, too.
By setting the analysis condition as 1.5km above sea level, the true air speed of the flight is 11m/s, and the optimization target is the mounting angle of the horizontal stabilizer. We adjust the installation angle of the horizontal stabilizer to obtain the deflection angle of the elevator. Trend of elevator deflection angle with mounting angle of horizontal stabilizer is shown in figure 15. When the elevator deflection is 0 • ,   the corresponding mounting angle of horizontal stabilizer is −5.87 • , which can be set as the optimized mounting angle of the horizontal stabilizer. We made a second flight after the optimization, and table 5 shows the second flight results compaired with the optimized simulation results. We can know that the simulation results are in good agreement with the second flight results.
The variation trend of elevator efficiency with dynamic pressure before and after optimization is shown in figure 16, and the trimming elevator deflection angle and angle of attack are shown in figure 17.
As can be seen from figure 16 and 17, after structural optimization, the efficiency of elevator is significantly improved.  At a low altitude of 1.5km and a flight speed of 11m/s, the efficiency of elevator is increased by about 15%. As the longitudinal trim of UAV is mainly carried out by the installation angle of the horizontal stabilizer, the elevator of the horizontal stabilizer has a small deflection angle, basically within 2 • , and the efficiency of the elevator is relatively high. Due to the small area of the tail wing, the lift provided by the tail is smaller than that of the wing, so the local optimization of the mounting angle of horizontal stabilizer has little influence on the trim angle of attack, and the change curves of the attack angle with the dynamic pressure of the UAV before and after the optimization are basically coincident.
For the flight speed of 11m/s, figure 18 shows the variation trend of the torsion angle of the wing and tail along with the wingspan before and after the optimization. The results indicate that the maximum torsion angle of the horizontal stabilizer is 10.7 • in level flight, and 4.9 • after optimization. Additionally, as the area of the tail is relatively small compared with the area of the wing, the change of local aerodynamic force of the tail has little influence on the wing, thus the torsion Angle of the wing before and after optimization does not change much, and the maximum torsion angle of the wing tip is about 2.76 • .

D. ANALYSIS OF TRIM AT 20km HEIGHT
In this subsection, the static aeroelastic analysis at 20km altitude is carried out for the optimized model. The conditions are set as follows: the flight overload is 1G, altitude is 20km, flight speed is 32m/s, dynamic pressure is 45.5Pa. The trimming results of rigid and elastic state are shown in table 6. The calculation results show that the trimming deflection angle of the elevator is −3.98 • (downward deflection) and the deflection angle of the rigid trim is 3.86 • (upward deflection) at the height of 20km.Therefore, the stiffness of the elevator surface itself has a great influence on the trim deflection angle, which should be paid more attention in the subsequent optimization and improvement. The deformation of the wing in trimming state is less than 5%, and the wing has strong stiffness.

V. CONCLUSION
This paper investigated a longitudinal static aeroelastic analysis and optimization scheme of HALE SP-UAV with high aspect ratio. The simulation results are well match with the real flight test results. Conclusions are as follows: 1) By analyzing the influence of elastic effect on the wing load distribution of high-aspect ratio SP-UAV, it is found that the shear force and bending moment in the elastic state will decrease compared with the rigid state, in which the bending moment decreases by about 6%.
2) We studied the trend of efficiency of the elevator vary with the pressure of the flight, and optimized the elevator efficiency problem, which can be controlled within 25% in level flight condition. However, the elevator efficiency loss of SP-UAV under large dynamic pressure is still up to 40%. Therefore, for high-aspect ratio UAV, the loss of rudder efficiency caused by elastic influence cannot be ignored and needs to be evaluated. Additionally, the elastic control method should be considered in flight control.
3) As the main structure of the wing adopts the single beam and multi-rib structure, the wing will produce positive torsional deformation to make the wing raise its head under the influence of elastic effect. Therefore, the angle of attack for the elastic trimming state will decrease correspondingly.
4) The optimized results of the mounting angle of horizontal stabilizer is about 4.9 • . Additionally, the elevator efficiency can be improved by improving the structural stiffness of itself, but this will bring additional structural weight gain, which needs to be evaluated comprehensively. 5) Considering the low efficiency of the elevator under the large dynamic pressure, it is necessary to reduce the flight speed to reduce the dynamic pressure, so as to improve the efficiency of the elevator during flight, but it also brings new challenges to the wind resistance of the SP-UAV. The SP-UAV is generally powered by propellers, which are distributed in parallel in front or rear of the wing. As to future research, it is worth researching the influence of propeller slipstream on the static aeroelastic analysis of SP-UAV.