Aerodynamic Load Of An Aircraft With A Highly Elastic Wing

doi:10.14311/AP.2017.57.0272 Acta Polytechnica 57(4):272–281, 2017 ©Czech Technical University in Prague, 2017 available online at http://ojs.cvut.cz/ojs/index.php/ap AERODYNAMIC LOAD OF AN AIRCRAFT WITH A HIGHLY ELASTIC WING Pavel Schoř Institute of aerospace Engineering, Faculty of Mechanical Engineering ,Brno University of Technology , Technicka 2896/2, 616 69 Brno Czech Republic correspondence: [email protected] Abstract. In this article, a method for calculation of air loads of an aircraft with an elastic wing is presented. The method can predict a redistribution of air loads when the elastic wing deforms. Unlike the traditional Euler or Navier-Stokes CFD to FEM coupling, the method uses 3D panel method as a source of aerodynamic data. This makes the calculation feasible on a typical recent workstation. Due to a short computational time and low hardware demands this method is suitable for both the preliminary design stage and the load evaluation stage. A case study is presented. The study compares a glider wing performing a pull maneuver at both rigid and and elastic state. The study indicates a significant redistribution of air load at the elastic case. Keywords: aerodynamics, panel method, finite element method, wing load, fluid structure interaction. 1. Introduction 1.1. Problem overview A flight of a modern glider or a High Altitude Long Endurance (HALE) aircraft is a representative example of the fluid-structure interaction problem. These flying vehicles are a result of a search for the most aerodynamically effective shape. This leads to a need for a wing with a large wingspan and an aspect ratio. Most recent gliders have a wingspan ranging from 18 to 30 meters with an aspect ratio approximately from 30 to 50. A NASA Helios HALE aircraft had a wing with a wingspan of 75 m and aspect ratio of 30 [ 1 ]. During a flight of such aircraft, the wing undergoes large deformations, still in the elastic regime, as shown on Figure 1. The Helios HP-03 experimental HALE aircraft crashed on June 26, 2003 after encountering atmospheric turbulence. The aircraft developed a wing Figure 1. The Helios aircraft undergoing large deformation [1]. tip displacement of more than 12 m and consecutive pitch oscilations, resulting in structural failure of the leading edge structure, but leaving the main spar undamaged. After the crash of Helios HP03 prototype, see Figure 2, a demand rose for new methods, which are capable to compute loads of highly elastic, “morphing” aircrafts [ 1 ]. In this article, a method for a solution of quasi-static load cases of such aircraft is presented. Moreover, recent airworthiness regulations explicitly require to calculate the load redistribution in case “if deflections under load would significantly change the distribution of external or internal loads” [2, 3]. 1.2. Methods overview Various methods are traditionally used in order to capture the change of the aerodynamic load generated in rigid and elastic state. The simplest option is to combine a vortex lattice method [ 4 ] with a simple beam finite element model [ 5 ] . The most sophisticated approach consists of Navier-Stokes CFD solver Figure 2. The Helios aircraft after structural failure during flight [1] 272 vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing for aerodynamic input combined with a detailed FEM structural model with solid elements [ 6 ]. The aerodynamic load is then directly interpolated to the FEM model. Both of these approaches need to iterate between the aerodynamic and structural model to obtain the equilibrium. The simple model with the lifting line theory provides the fastest solution times, however, errors may be introduced because of an improper pitching moment distribution due to the simplified VLM aerodynamic model. The detailed model with the Euler or Navier-Stokes CFD solver with a detailed FEM structural model can provide the most accurate prediction of the load redistribution. Typically, the CFD volume mesh has to be updated according to the results from the FEM solution. This may cause issues with the deformation of the CFD volume mesh when the deformation of the wing is large. Recently an Immersed Boundary Method is used to solve this kind of fluid-structure interaction [7] When assuming the costs of the computational time of the CFD solver and also assuming that up to 10 iterations are usually required to obtain the converged solution for one particular angle of attack, the use of the CFD solver is restricted to resend high performance computer clusters, therefore making it very expensive and not feasible for a typical glider manufacturer neither at a design stage nor at a load evaluation stage. In this article, a 3D panel method is used, combined with a geometrically nonlinear finite element method with beam elements to model the structure of the wing. The savings of computational time when compared to the CFD N-S solver and the differences in computed air loads justify the use of this low fidelity aerodynamic tool, making it possible to compute large number of load cases on an average workstation. 1.3. Eligibility Based on assumptions made in Sections 2.1 and 2.2, the method presented in this article is suitable for analysis of quasi-static maneuvers of an aircraft with elastic wings, operating at a low Mach number ( M < 0 . 3) and sufficiently high Reynolds number ( Re > 1000000). This includes the HALE aircraft during its low altitude flight or a glider performing a pull maneuver. 2. Methods 2.1. Aerodynamic model A fluid motion is traditionally described by set of equations called Navier-Stokes equations [ 8 ]. The number of the exact solution of Navier-Stokes equations is small. For practical problems, these equations are solved numerically [9]. Solving the Navier-Stokes equations provides a detailed description of the flow, however the computational costs are enormous. Therefore some assumptions are made to reduce the complexity of the problem. If the flow is assumed to be steady, inviscid, Figure 3. Potential flow domain. incompressible and irrotational, is called a potentital flow and it can be described by Laplace’s equation [ 8 ] : ∇2Φ=0 (1) Where Φis a scalar function called velocity potential and velocity −→ q at each point can be obtained by its derivation: −→ q=∇Φ(2) An approach from [ 10 ] is used to solve the Laplace’s equation. Assume a cross section of a wing as shown in Figure 3. The surface S∞ encloses the problem at infinity, surface S represents a body (wing) and surface W represents its wake. The surface S + W divides the domain into two regions: the external region with the flow field of interest and velocity potential Φand the internal region with the fictitious flow and velocity potential Φ i . The surfaces are modeled by a doublet and source singularities. Green’s theorem is applied to both the outer and inner region and velocity potential Φ P is obtained by combining both expressions. [ 10 ] The velocity potential Φ P gives a velocity potential anywhere in the two regions. It is expressed in terms of surface integrals in terms of velocity potential and its normal derivative over the boundary surface as ΦP=1 4πZZS+W+S∞ (Φ −Φi)−→ n· ∇ 1 rdS −1 4πZZS+W+S∞ 1 r −→ n·(∇Φ− ∇Φi) dS. (3) Where r is the distance from the point P to the element d S on the surface, and −→ n is a unit normal vector of the element d S . [ 10 ] The first integral in (3) represents the disturbance potential from a surface 273 Pavel Schoř Acta Polytechnica distribution of doublets with a density of (Φ − Φ ∞ ). The second integral represents the contribution of sources with a density of −−→ n·(∇Φ− ∇Φi).[10] Next, the internal Dirichlet boundary condition is introduced. This boundary conditions sets the Internal flow equal to the Onset flow by Φi= Φ∞,(4) ΦP= Φ∞.(5) Following the procedure from [ 8 , 10 , 11 ], the surface is discretized into n surface panels and nw panels in the wake region with a constant distribution of singularities. The (3) is also discretized and the Dirichlet boundary condition is evaluated at centroid of each surface panel. These points are called Collocation points. The result of this procedure is a linear equation [8] n X k=1 Ckµk+ nw X l=1 Clµl+ n X k=1 Bkσk= 0,(6) where Ck , Cl are doublet influence coefficients of body and wake panels, Bk is the source influence coefficient, µk is the doublet strength of a body panel, µl is the doublet strength of a wake panel, σk is the source strength of a body panel. The σkis set as σk=−→ nk·−−→ Q∞,(7) where −→ nk is the panel normal unit vector and −−→ Q∞ is the freestream velocity. The method described above is known as the “first order panel method” [ 8 ]. It is used as the only source of aerodynamic data for the fluid-structure interaction. The formulation of the panel method is identical to well-known and validated codes such as the VSAERO [10] and PMARC [11]. Boundary layer. The potential flow may give useful results for some applications, however for applications where viscous effect become more important, the results obtained under assumption of potential flow may become inaccurate in terms of over predicted lift curve slope. A typical example of a flow with a significant viscous effect is a glider wing at a low Reynolds number or a flow over a multi-element wing. The viscous effects can be introduced into the potential flow model by displacing the actual body surface by a displacement thickness of the boundary layer [ 12 ]. The actual displacing of the body surface is performed by a modification of panels source strength in (7) to σk=−→ nk·−−→ Q∞+d ds(ueδ∗),(8) where sis a distance on the body, ueis a panel edge velocity and δ∗ is a displacement thickness computed by the boundary layer analysis [12]. 2.2. Structural model The structure of the wing is modeled using the FEM with beam elements only. Each node of the element has six degrees of freedom. The stiffness properties are assumed constant for each element. The stiffness matrix of such element can be found in many FEM related textbooks [ 13 ]. However such elements are usually derived for isotropic materials. In fact, wings of modern gliders are made of a carbon fiber reinforced composite orthotropic material. A wing made of such material exhibits coupling between bending and torsion. This causes issues, when one tries to reduce stiffness properties of a composite wing into single stiffness matrix of a beam element. A method from [14] was adopted to solve this issue. As the wing undergoes large deformations it is necessary to capture the geometrical non linearity. To solve this issue a Newton/Rhapson iterative method is used. The geometrical non-linearity, however, causes difficulties when calculating elements stiffness matrices — one must determine the derivation of the global stiffness matrix (tangent stiffness matrix). A method from [ 15 ] was adopted, which provides a numerically generated tangent stiffness matrix for truss and beam elements. 2.3. Aero-structural interaction Work flow. A brief overview of the interaction between the aerodynamic and internal forces in the wing structure is shown on Figure 4. The interaction is done by using a modified Newton-Rhapson method. This method and steps called the “Compute aerodynamic forces” and “Compute wing displacement” from Figure 4 are explained in detail in the paragraph Load control scheme. Computing the rigid lift curve as a first step is not necessary, the main reason for this is to find the zero lift angle of the attack α0 , where the computation of elastic lift curve starts. Starting from α0 makes the wing almost unloaded and the Newton-Rhapson method is likely to converge. Also, the increment of the angle of attack ∆ α should be small enough to assure the convergence, it is usually from 0.5°to 2°. Interacting element. The interacting element is explained in Figure 5. It consists of aerodynamic panels, but there is only one element in the span-wise direction and one FEM beam element. Forces and moments are always integrated from the wingtip to plane of symmetry, at a point which is also a node of the finite element model. The load applied at the FEM node is a sum of forces moments from all panels belonging to the interacting element as shown on Figure 5: −→ f= n X i=1 −→ fi,(9) −→ m= n X i=1 −→ fi×ri,(10) 274 vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing Figure 4. Aero-structural interaction — work flow. where fi is the force on i th panel in the global coordinate system x, y, z , as shown in Figure 5; ri is the vector between the reference node and control point of ith panel. The result of the FEM analysis is a displacement vector u at each node. Coordinates of aerodynamic panels are updated at each load step using nodal translational and rotational displacements, according to (11). Also, the nodal coordinates are updated as explained in [15]. PTSnew =Rx(∆φx)Ry(∆φy)Rz(∆φz)PTSold + ∆U, (11) where ∆ φx ,∆ φy ,∆ φz, are rotational increments in the global coordinate system calculated from the last step of the Newton-Rhapson iteration, see bellow; ∆ U is the displacement increment; Rx ( φx ) , Ry ( φy ) , Rz ( φz ) are rotational matrices defined by Rx(∆φx) =   1 0 0 0 cos φx−sin φx 0 sin φxcos φx  ,(12) Ry(∆φy) =   cos φy0 sin φy 0 1 0 −sin φy0 cos φy  ,(13) Rz(∆φz) =   cos φz−sin φz0 sin φzcos φz0 0 0 1  .(14) Figure 5. Panel method — FEM interaction. Load control scheme. In the Newton/Rhapson method, the total load is divided into m load steps and for each load step, a Newton iteration is done until an equilibrium between the internal forces and external load is found. After that, the external load is increased by the load step. However, in the case of an elastic wing, the total external load is unknown because the load (lift) also depends on how much the wing is deformed. Therefore, the standard load control in the Newton/Rhapson method was changed in the following way: The load is controlled by increasing the angle of attack as explained in Figure 4. Usually the lift generated by the wing (sailplane) is the dependent variable in the u – F curve. To compute the load generated by the elastic wing, a following scheme is used: (1.) The first iteration starts at an angle of attack of zero lift of a rigid wing. The equilibrium between the external load computed from the panel method and internal forces is computed using the Newton iteration. (2.) The angle of attack is increased by a small step, e.g., 1 ° and the equilibrium between the aerodynamic load and the internal forces is computed by the Newton iteration. (3.) The angle of attack is further increased until the total load (lift) is equal to the desired load (lift). The result is a nonlinear lift curve. For each load step, the vector of internal forces is also a nonlinear function of previous load steps. Therefore when computing a load of an elastic wing it is convenient to compute a wing lift curve as described above with very small increment of the angle of attack, e.g., 0 . 5 ° and to save the vector of deformation and internal forces for each converged load step. When computing large number of load cases, significant savings in the computational time are achieved when the pre-computed deformation vector and internal forces vector are loaded for the closest lower angle of attack. One then needs to compute only the difference between the pre-computed angle of attack and the desired angle of attack. 275 Pavel Schoř Acta Polytechnica Figure 6. Aero-structural interaction — load control in Newton-Rhapson method. Figure 7. Coordinate systems. 3. Numerical study 3.1. Overview A numerical study was conducted to demonstrate the importance of the elasticity of the wing for its load. Assuming a flight envelope from [ 2 , par. CS22.333], and considering the point A of the envelope to be critical, the aerodynamic load at this point is evaluated using three different methods: (1.) panel method with elastic wing; (2.) panel method with rigid wing; (3.) lifting line theory. The lifting line theory is used as a validation tool for the panel method. 3.2. Coordinate systems A global right-handed coordinate system OXY Z defined in Figure 7 is used. Additional local coordinate systems oixiyizi are used for evaluation of the local load. As the wing deforms, these coordinate systems Property Symbol Value Area A 11.55 m2 Span L 21.0 m Root chord cr0.7 m Tip chord ct0.4 m Dihedral Γ0 rad Sweep Θ0 rad Table 1. Wing geometry. Property Symbol Unit Cross section area ASm2 Second moment of area Ixm4 Second moment of area Izm4 Torsional constant Jym4 Offset of neutral axis ∆Xnm Offset of neutral axis ∆Znm Offset of elastic axis ∆Xtm Offset of elastic axis ∆Ztm Young modulus EN/m2 Shear modulus GN/m2 Table 2. Wing cross sectional properties. follow the deformed geometry. The free stream velocity V with positive angle of attack is also shown in Figure 7. 3.3. Wing specimen A simple tapered wing is examined in this study. The wing is designed as a planar wing in undeformed state. The wing has one straight main spar located at 0 . 25 c , where c is local chord. The main reason for this geometry is the evaluation of the load by the lifting line theory, which can not be used for a non-planar wing geometry and nor swept wings. 276 vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing Figure 8. Spanwise variation of Iz,Ix,Jy. Figure 9. Airfoil E603 inviscid aerodynamic data. Geometry. Geometry of the wing is described in Table 1. An Eppler E603 airfoil is used for all sections of the wing. Stiffness properties. When analyzing the deformation of the wing by the FEM analysis described above, an instrumental step is the determining of all cross-sectional characteristics defined in Table 2 for all cross sections of the wing. In this example, the Young’s modulus E and shear modulus G are assumed as constant for the whole wing with following values: E= 60 GPa, G = 7 GPa. The wing has only one spar at 0 . 25 c and the nodes are coincident with the spar. Therefore the offset of the elastic axis is almost negligible and both values are assumed to be zero. The neutral axis is located at the local centroid with non-zero offsets, but these offsets are not shown in this paper. The spanwise variation of the Second moments of area and torsional constant is shown on Figure 8. Finite element model. The wing is divided into 162 spanwise sections. Each section has one twonoded beam element described in Section 2. Nodal coordinates always correspond to a location of a main spar. The node located in the plane of symmetry has all degrees of freedom prescribed to zero. All coupling terms in the element stiffness matrices related to the bending-torsion coupling are intentionally set to zero for clarity. Airfoil data. The aerodynamic characteristics of the Eppler E603 airfoil were obtained from the XFOIL inviscid analysis, as shown in Figure 9. This aerodynamic data are used only for the lifting line analysis. 3.4. Loadcase VA Assuming a maximum lift coefficient of the wing CLmax = 1 . 4, from Figure 10, the flight conditions at point A of the maneuvering envelope are shown in Table 3 The total mass of the glider is MTOW = 750 kg, but inertia forces acting on the wing are not assumed in the analysis for clarity. 277 Pavel Schoř Acta Polytechnica Figure 10. Comparison of rigid and elastic wing — lift. Property Symbol Value Airspeed V64.0 m/s Wing lift coefficient CL 1.4 Load factor n5.3 Table 3. Load case VA. Deformation sequence. Deformation sequence for various wing lift coefficient is shown in Figure 11. It must be noted that this deformation sequence is valid only for airspeed V = 64 . 0m/s. As the deformation of the wing depends on forces and moments applied on the wing surface, for a different airspeed, the corresponding deformed shape will be different. Integral aerodynamic characteristics. Integral aerodynamic characteristics of the rigid and elastic wing are compared on Figures 10 and 12. A close agreement for rigid cases can be observed. For the elastic case a significant decrease of the slope of the wing lift curve can be seen in Figure 10. Additionally, the increase of pitching moment of the elastic wing can be seen in Figure 12. The differences between rigid and elastic cases are discussed in Section 4 Load comparison. Actual comparison of the load for the rigid and elastic wing performing a steady maneuver with the same load factor n = 5 . 3is shown in Figures 13–15. Only bending and torsional moments are shown. The differences between the rigid and elastic case are discussed in Section 4 4. Discussion 4.1. Validation of results There are two main options to validate the result computed by the presented method: Figure 11. Deformation sequence at VA. Numerical validation. A numerical validation means to prerform the case study using a RANS CFD solver as a source of aerodynamic data coupled to a wing modeled by shell elements or by 3D solid elements. Experimental validation. Experimental validation means to actually fly a glider, perform the pull maneuver and measure the displacement of the wing. Conclusion on validation. Unfortunately both numerical and experimental validation would require much more effort and funds than what was spent for this whole work. A literature review also did not give any example, which could be used as a validation benchmark. Therefore, this work is presented only as a numerical study without any other validation. 4.2. Decrease of wing lift curve slope As the lift generated by the elastic wing increases, the wing bents more, as shown in Figure 11. As the wing bents more, the slope of the lifting curve drops. This is demonstrated in Figure 10, which shows a reduced slope of the lift curve of the elastic wing, compared to the rigid wing. This phenomena could possibly be explained by the reduction of the projected area of 278 vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing Figure 12. Comparison of rigid and elastic wing — pitching moment. Figure 13. Comparison of rigid and elastic wing — Normal bending moment Mx. the wing and the change of direction of the lift force. This has following consequences: (1.) When the prescribed total lift is reached, the wing operates at a higher angle of attack. In the particular case presented in Section 3 the difference of the critical angle of attack is 6°. (2.) The normal bending moment is slightly reduced, see Figure 13. The tangential bending moment is, however, increased due to the increased angle of attack, see Figure 15. In the particular case presented in Section 3 the tangential bending moment of the elastic wing is 1 . 34 times higher, than in the rigid case. The normal bending moment is, however, 0 . 87 times lower. 4.3. Increase of wing pitching moment The explanation of the increase of the wing pitching moment requires more labor than the explanation of the decrease of the wing lift curve slope. First, assume a wing cross section located at the wing tip. The wing is of an undeformed shape and operates near the critical angle of attack. The section generates a pitching moment around the main spar MY 1 = 1 2ρv2S1c1cm1 ( α ). The resultant of the aerodynamic force is assumed to act at the main spar (0 . 25 c ). The direction of the force is upward and forward (negative X , positive Z ). When the wing has zero dihedral and zero sweep, the contribution of this force to the pitching moment is zero, since the vector of the section reference point related to the wing reference point ~r1 = [0 ,− 10 . 37 , 0]. The wing reference point is defined as a point (node) on the main spar in the XZ plane. Next, assume a highly deformed wing shown in Figure 16. The wings tip section is displaced by a vector ~ d1= [−0.12,2.5,5.15], the vector of the section reference point related to the wing reference point is now ~r1= [−0.3,−8,5.15]. The resultant of the aerodynamic force at the wing 279 Pavel Schoř Acta Polytechnica Figure 14. Comparison of rigid and elastic wing — Torsional moment My. Figure 15. Comparison of rigid and elastic wing — Tangential bending moment Mz. reference point is computed as ~ f1= [−6.07,49.38,9.31], then, the resulting moment at the wing reference point generated by the force ~ f1is ~m1=~r1×~ f1= [−329.05,−28.51,−63.32]. When the pitching moment of the section is evaluated, the contribution of ∆ MY 1 = − 28 . 51 Nm causes to increase the pitching moment coefficient of the section by a value of − 0 . 55 mostly due to the vertical displacement of the section and the force resultant directing forward. 4.4. Wing maximum lift In this article, the maximum lift coefficient of the rigid wing was computed using the nonlinear lifting line theory as CLmax = 1 . 4. The airfoil sectional data used in the analysis were computed by the XFOIL inviscid analysis, assuming to be close to the experimental data. Since the purpose of this article is to present a Figure 16. Pitching moment of highly deformed wing. potential problems of an aircraft with a highly elastic wing, rather than to make a prescription based on practical results, the maximum lift coefficient of the elastic wing is also assumed to be CLmax = 1.4. 280