Static Aeroelastic Beam Model Development for Folding Winglet Design

: Wing shape adaptability during ﬂight is the next step towards the greening of aviation. The shape of the wing is typically designed for one cruise point or a weighted average of several cruise points. However, a wing is subjected to a variety of ﬂight conditions, which results in the aircraft ﬂying sub-optimally during a portion of the ﬂight. Shape adaptability can be achieved by tuning the shape of the winglet during ﬂight. The design challenge is to combine a winglet structure that is able to allow the required adaptable shape while preserving the structural integrity to carry the aerodynamic loads. The shape changing actuators must work against the structural strains and the aerodynamic loads. Analyzing the full model in the preliminary design phase is computationally expensive; therefore, it is necessary to develop a model. The goal of this paper is to derive an aeroelastic model for a wing and winglet in order to reduce the computational cost and complexity of the system in designing a folding winglet. In this paper, the static aeroelastic analysis is performed for a regional aircraft wing at sea level and service ceiling conditions with three degree and eight degree angle of attack. MSC Nastran Aeroelastic tool is used to develop a Finite Element Model (FEM), i.e., beam model and the aerodynamic loads are calculated based on a doublet lattice panel method (DLM).


Introduction
Aeroelasticity is a discipline that studies the phenomena that involve the interaction between aerodynamic forces and elastic forces. This study is crucial due to two main reasons. First, the interaction can lead to static or dynamic instabilities that can result in failure of the structure. Secondly, the loads calculated using aeroelastic simulations are significantly different for flexible structures as compared to decoupled structural and aerodynamic simulations. This has significant influence on the design and weight of the structure, and its aerodynamic performance. Aeroelasticity is divided in to two main groups as Static and Dynamic aeroelasticity. Static aeroelasticity is the study of the deflection of flexible aircraft structures under aerodynamic loads, where the forces and motions are considered to be independent of time. Dynamic aeroelasticity studies the interactions among aerodynamic, elastic, and inertial forces. Examples of dynamic aeroelastic phenomena are Flutter, Buffeting, and Transonic aeroelasticity. General Form of the Aeroelastic Equations is Aq + (ρVB + D)q + (ρV 2 C + E)q = 0 (1) where q is the generalized coordinates, A is the structural inertia, B is aerodynamic damping, C is aerodynamic stiffness, D is structural damping, and E is the structural stiffness matrices. Equation (1) is one of the most important equations and it describes the fundamental interaction between the flexible structure and the aerodynamic forces [1].
Morphing is a capability to provide superior performance while in flight by tailoring the vehicle's state to adapt to the external operating environment and multivariable mission roles. Morphing aircraft are multi-role aircraft that change their external shape in both small-scale and large-scale substantially to adapt to a changing mission environment during flight. Bowman et al. identified two distinct classes of morphing: Morphing for Mission Adaptation, and Morphing for Control. Morphing for Mission Adaptation is a large-scale, relatively slow, in-flight shape change to enable a single vehicle to perform multiple diverse mission profiles. On the contrary, Morphing for Control is often a small-scale or component level in-flight physical or virtual shape change, being used to achieve multiple control objectives, such as noise suppression, flutter suppression, load alleviation, and active separation control [2].

Morphing Winglet with Aerodynamic Actuation
Morphing wings have large potential to improve the overall aircraft performances, by dynamically optimizing the shape to various flight conditions. The morphing winglet is a promising solution among the morphing technologies due to its small size and large effect on the aerodynamics. In this paper, the morphing winglet has to change its cant angle, and its stiffness has to be large enough to carry loads. While increasing the total stiffness,the allocation of the stiffness can be unsymmetrical, introducing deformation under a linear actuation force. A simplified beam model is developed to estimate the induced deformation. The deformation is also visualized using finite element analysis. To demonstrate the application of the morphing structure, the baseline wing design of a regional twin turboprop airliner ATR 72 is used [3,4].
There are different suggested concepts to actuate the folding winglet mechanism. Some which were addressed in previous works are: hinge line actuator, system of pulley gears, and a worm gear system. Each of these systems have their own advantages and disadvantages. The advantage of the hinge line actuator is that the direct connection between the actuators and the hinge line mitigates additional transmission losses, making the system more efficient. The disadvantage is the loss in the winglet area due to the mounting position of the actuator. An additional disadvantage due to its location is the actuator mass will maximize the root bending moment and is susceptible to damage should the aircraft have a tip strike upon landing or take off. The advantage of the system of pulley gears is ability to bring actuation component inboard by maximizing the winglet surface area. The advantage of a worm gear system is the power off locking capability. Power is not required for the system when remaining in one position, because the wheel cannot rotate unless the screw is actively driven. The worm drive acts to gear the system, enabling larger torques to be developed. This could be advantageous by allowing a less powerful (smaller) actuator to be used, but it will also make the system slower. Figures 1 and 2 graphically represent the proposed concepts [5,6].
In this project, the topic of investigation is applying a secondary shape changing mechanism that will alter the aerodynamic loads before the folding angle shape change. The alteration of loads should be such that the aerodynamic loads themselves change the cant angle of the winglet. The envisioned shape changing mechanism to achieve this is a twisting winglet. In summary, the change in folding angle of the winglet will be achieved by first twisting the winglet, changing the local aerodynamic loads such that they will work with the folding angle change rather than work against it. The advantage of aerodynamic actuation is a reduction in the actuator cost and weight. The design challenge is to combine a winglet structure that is able to allow the required adaptable shape while preserving the structural integrity to carry the aerodynamic loads. The shape-changing actuators have to work against the structural strains and the aerodynamic loads. This contradictory requirement often leads to rather heavy actuators that are needed to effectuate the desired shape change [7].

Materials and Methods
The flow chart in Figure 3 shows the approach taken while performing the modeling and analysis of the problem. The first step is to import the wing and perform clean up and prepare it to the structural modelling. In the structural model, the geometry is divided into elements and the material properties are distributed. Later, the lifting surface definition and panel box divisions are performed in the aerodynamic modeling. Once we have the structural and aerodynamic grids, the interpolation is carried out using theory of infinite beam splines. Finally the analysis is performed and the desired results are viewed at the post processing. In Figure 4, a three-dimensional (3D) base wing CAD model and two-dimensional (2D) wing with winglet is represented. All CAD models are done on SolidWorks.

Structural Model
The structure is modelled as a beam element. Generally, beam models are developed based on various assumptions that lead to different level of accuracy. Several beam theories have been developed; the Euler-Bernoulli beam theory, first described by Euler and Bernoulli, is one of the simplest and most useful theories. A fundamental assumption of this theory is that the cross-section of the beam is infinitely rigid in its own plane, i.e., no deformations occur in the plane of the cross-section. Consequently, the in-plane displacement field can be simply represented by two rigid body translations and one rigid body rotation. This fundamental assumption deals only with in-plane displacements of the cross-section. Two additional assumptions deal with the out-of-plane displacements of the section: during deformation, the cross-section is assumed to remain plane and normal to the deformed axis of the beam. While these assumptions work fine for thin beams under static loading, for thicker beams and higher frequencies they are an oversimplification. In order to overcome these restrictions, the more complex Timoshenko Beam theory was formulated in 1921 and extends the Euler-Bernoulli theory by taking shear deformation and the rotary inertia into account [8].

Beam Section
There are different techniques for reducing detailed finite elements models, such as static condensation of mass and stiffness matrices, stick models, etc. The stick models are simplified beam finite element models commonly used in civil aircraft design and multidisciplinary design optimization (MDO). The accuracy of the stiffness characteristics in its model highly affects the accurate prediction of bending and twisting deformations of the aircraft structure in flight.
The stick model is usually generated by extracting the stiffness properties of the main structure and applying it to a set of beam elements extending along the structure's elastic axis. This work exposes a methodology to extract section properties for beam elements while using the Femap surface property generator. Femap calculates all the surface properties by selecting the desired surface. Subsequently, those properties are transferred into the beam element. There are different levels of structural modelling, 3D finite element model is always prepared once the aircraft's layout and structural details are obtained. The three dimensional finite element model is commonly used in structural design validation and optimization. Imposing a detailed 3D FEM in a multidisciplinary iterative process or in loads analyses is quite expensive. Alternatively, a simplified beam FEM of the wing can be generated and employed in the iterative process and dynamic analyses.

Extracting Beam Sectional Properties
This method extracts the stiffness properties of the beam model from the complete 3D FEM model. The method requires some processing in Femap with Nx Nastran and requires the user to be knowledgeable of the software package. The method assumes that the stiffness extraction process begins by locating the reference coordinate system at the root section of the detailed FEM model; there are no constraints and this coordinate system could be located anywhere in the section. To make the description of the procedure easier, the coordinate reference system is the same for all the parameters, defined as the one of Figure 5; the general axis of the beam model is the z-axis [9].
where F is the vector of forces and moments, u is the displacements vector, [K] is the stiffness matrix, and [FF] is the flexibility matrix. The wing and winglet geometric properties and beam characteristics are mentioned in Tables 1 and 2. Aluminum 2024 was the material of choice and the material properties are listed in Table 3.    CBEAM entry defines a beam element in Nastran and its properties are defined with a PBEAM, PBCOMP, or PBEAML entry. The beam element includes extension, torsion, bending in two perpendicular planes, and the associated shear.
Tensile stresses are given a positive sign and compressive stresses a negative sign. Only the longitudinal stresses are available as complex stresses. The stress recovery coefficients on the PBEAM entry are used to locate points on the cross section for stress recovery. Because there is no torsional stress recovery for the CBEAM element, the margin-of-safety computation does not include the torsional stress [10].
The blue curve in Figure 6 shows the beam element model of the wing and winglet.

Doublet-Lattice Subsonic Lifting Surface Theory (DLM)
The DLM is an aerodynamic finite element method for modelling oscillating lifting surfaces that reduces to the Vortex-Lattice Method at zero reduced frequency. Each of the interfering surfaces (or panels) is divided into small trapezoidal lifting elements ("boxes"), such that the boxes are arranged in strips parallel to the free stream with surface edges, fold lines, and hinge lines lying on box boundaries. Uniform concentration of the unknown lifting pressures across the one-quarter chord line of each box is assumed. On the three-quarter chord line of the box there is one control point per box centered span-wise and the surface normal wash boundary condition is satisfied at each of these points. The number of finite elements (boxes) required for accurate results depend on different parameters, such as the aspect ratio and reduced frequency. At high reduced frequency, the chord-wise dimension of the boxes must be small. However, an approximation in the method that the variation of the numerator of the incremental oscillatory kernel function is parabolic across the span of the box bound vortex, restricts the box aspect ratio to approximately 3. If a higher order approximation is used for the span-wise variation of the numerator of the incremental oscillatory kernel, the limitation on box aspect ratio can be relaxed and the number of span-wise divisions required in high-frequency analyses will be reduced significantly, leading to a reduction in the total number of boxes. The DLM is a FEM for the solution of the oscillatory subsonic pressure-normal-wash integral equation for multiple interfering surfaces.
where w is the complex amplitude of the dimensionless normal wash, p is the complex amplitude of the lifting pressure coefficient, (x,s) is the orthogonal coordinates on the nth surface S, such that the undisturbed stream is parallel to the x-axis, and K is the complex acceleration potential kernel for oscillatory subsonic flow. The original DLM algorithm was presented at the same time (1968) as the Lifting Line Element Method (LLEM) of Stark R . Although numerous comparisons@ with experiments were shown at the time, the complete details of the LLEM were never published. A refinement to the expressions for the kernel given by Rodemich R and Landahl was presented by Rodden, Giesing, and Kalman in the form K = (K 1 , T 1 , /r 2 + K 2 T * 2 /r 4 )exp(−iωx/U) (4) in order to analyze the non-planar interference correctly. K 1 , and K 2 , are the planar and non-planar parts of the kernel numerator, respectively. Here, ω is the frequency, x 0 is the distance between the sending and receiving points parallel to the freestream, and U is the velocity of the freestream, T * 2 = (z 0 cosγ r − y 0 sinγ r )(z 0 cosγ s − y 0 sinγ s ) where y 0 and z 0 are the Cartesian distances between sending and receiving points perpendicular to the freestream, and γ r and γ s are the dihedral angles at the receiving and sending points, respectively. The description of K 1 and K 2 as the planar and non-planar parts of the kernel numerator is a convenience, because both are obviously non-planar in general. The refinement in the equation was found necessary so that the DLM could predict the interference between a nearly planar wing and horizontal tail, The refinement retained the original primary approximation, i.e., that the incremental oscillatory normal wash factors are obtained by integrating the difference between the oscillatory and steady kernels over the length of the bound vortex assuming a quadratic (parabolic) variation in the numerator of the difference. The total normal-wash factor is then the sum of the incremental oscillatory normal-wash factor and the steady normal-wash factor obtained from the expressions for a horseshoe vortex, e.g., the Vortex-Lattice Method (VLM) of Hedman. In this way, the DLM converges to the VLM at zero reduced frequency, and the error in the parabolic approximation of the kernel numerator difference is small at low reduced frequencies and increases with the reduced frequency [11,12]. The Doublet-Lattice method (DLM) can be used for interfering lifting surfaces in subsonic flow. The following general remarks summarize the essential features of the method.

•
The theoretical basis of the DLM is linearized aerodynamic potential theory.

•
The undisturbed flow is uniform and is either steady or varying (gusting) harmonically.

Interconnection of the Structure with Aerodynamics
The connection between the structural grids and the aerodynamic grids is performed by an interpolation. Therefore, the grids for both aerodynamic and structural elements can be selected independently adapted to the particular theory. The structural model for a wing may involve a one-, two-or three dimensional array of grid points. A general interpolation method is available that will interconnect the various combinations. Any aerodynamic panel or body can be subdivided into sub-regions for interpolation, while using a separate function for each. The interpolation method is called splining. The theory involves the mathematical analysis of beams.
While solving aeroelastic problems, the structural degrees of freedom are the independent and the aerodynamic degrees of freedom are dependent. A matrix is derived that relates the dependent degrees of freedom to the independent ones. Two transformations are required: the interpolation from the structural deflections to the aerodynamic deflections and the relationship between the aerodynamic forces and the structurally equivalent forces acting on the structural grid points. The splining methods lead to an interpolation matrix that relates the components of the structural grid point deflections to the deflections of the aerodynamic grid points. This interpolation matrix [G kg ] defined below relates the components of structural grid point deflections u g to the deflections of the aerodynamic grid points u k , The derivation of the elements of [G kg ]is discussed in [11]. A general conceptualization of all the spline methods is that displacements are available at a set of points and it is desired to derive an interpolant that passes through all these points.
The general form of these interpolants is where α µ , µ = 1, Mandβ ν , ν = 1, N are unknown parameters and ψ u andψ v are given interpolation functions. This provides a spline matrix for a single spline. The final matrix that is required is one that combines all of the splines in the model into an overall spline matrix that relates all the aerodynamic degrees of freedom to all degrees of freedom.
The aerodynamic forces F k and their structurally equivalent values F g acting on the structural grid points therefore do the same virtual work in their respective deflection modes, where δu k and δu g are virtual deflections. substituting in the above equations and rearranging yields due the arbitrariness of the virtual deflections, the required force transformation is obtained.
to complete the formulation of aeroelastic problems Equations (8) and (13) are both required in which the aerodynamic and structural grids do not coincide. The transpose of the deflection interpolation matrix is all that is required to connect the aerodynamic forces to the structure; therefore, it can be said that the interconnection problem is simply a problem of interpolating from the structural to the aerodynamic grid points.

Theory of Infinite Beam Splines
A linear spline is a "beam" function, w(x), which passes through the known deflections, w = w(x i ) with twists φ i = φ(x i ). Bending deflections are easily solved by the three-moment method, which is appropriate for simple beams. However, an extension of the method is required for splines with torsion, rigid arms, and attachment springs. The blue spline on Figures 9 and 10 show the interpolation of the structural and aerodynamic grids.

Static Aeroelasticity
Static aeroelastic problems deal with a steady state response to the interaction of aerodynamic and structural forces on a flexible vehicle that results in a redistribution of the aerodynamic loading as a function of airspeed. Nastran computes the aircraft trim conditions, with the subsequent recovery of structural responses, aeroelastic stability derivatives, and static aeroelastic divergence dynamic pressures.
For static aeroelasticity, the downwash relation becomes: where w j is aerodynamic degrees of freedom, u k is aerodynamic displacements (deformation), u x is extra aerodynamic points used to describe, for example, control surface deflections and overall rigid body motions, w g j is an initial static aerodynamic downwash, [D jk ] is substantial derivative matrix for the aerodynamic displacements, and [D jx ] is substantial derivative matrix for the extra aerodynamic points.
The aerodynamic forces can be written as

Static Aeroelastic Equations of Motion
The aerodynamic forces are transferred to the structure using the spline matrix reduced to the a-set to form an aerodynamic influence coefficient matrix, Q aa , which provides the forces at the structural grid points due to structural deformations and a second matrix, Q ax , which provides forces at the structural grid points due to unit deflections of the aerodynamic extra points: The general equations can be written as: where K aa is the stiffness matrix, M aa is the mass matrix, and P a is vector of applied loads [11]. The flight condition parameters for the static aeroelastic analysis are presented in Table 4.

Results
In this section, the effect of changing different variables, like mesh size, number of panel divisions, angle of attack, and dynamic pressure on the static aeroelastic result is presented. The total translation, total rotation and stress values are obtained. These values are normalized in order to be able to compare the effect of change of variables. The non-dimensional stability and control derivative coefficients, aerodynamic force and pressure on the aerodynamic elements, structural monitor point integrated loads, and aerodynamic monitor point integrated loads are presented as an output f.06 dat file.

Change in Mesh Size
FEMAP has robust tools for mesh sizing and meshing. In addition to manually creating nodes and elements, it can extrude, revolve, and sweep meshes. The meshing toolbox has powerful tools for editing geometry and meshes. In this paper, three mesh sizes were used for the beam model. First, the mesh should be sized, therefore, the size on curve command is selected from the mesh control tool. After choosing the curve to mesh, the desired parameters are entered in the Automatic Mesh Sizing tool. The number of nodes and elements for different mesh sizes are presented in Table 5.

Change in Division of the Panel
Here, two panel divisions are compared. The first one has 600 boxes, i.e., the Chord is divided into 10 boxes and the span into 60 boxes for the base wing and 200 boxes for the winglet. The second one has 1350 boxes, i.e., the Chord is divided into 15 boxes and the span into 90 boxes for the base wing and 450 boxes for the winglet. Increasing the number of divisions increases the accuracy of the results, however it also increases the computational time.

Angle of Attack
Two different angle of attack values were chosen and inserted in the trim variables. For the first analysis, a 3 • AoA is chosen. In Femap, it has to be entered in radian. 3 • ×π 180 = 0.05235987. For the second case, α = 8 • is chosen, and it is also entered in radian. 8 • ×π 180 = 0.139626.

Dynamic Pressure
Dynamic pressure is the increase in a moving fluid's pressure over its static value due to motion. In compressible fluid dynamics, impact pressure (dynamic pressure) is the difference between total pressure (also known as pitot pressure or stagnation pressure) and static pressure. In aerodynamics notation, this quantity is denoted as q c or Q c .
The change in altitude causes change in the density, pressure, and temperature, which, in turn, changes the dynamic pressure.

Sea Level
Referring to the above Equation, the dynamic pressure at Sea Level can be calculated as:

Service Ceiling
The Service Ceiling of ATR 72-600 is at 7600 m. The dynamic pressure is calculated near the service ceiling at 7000 m as: Figures 11 and 12 show the deflection and rotation of the model respectively. Figures 13 and 14 represent the effect of dynamic pressure and angle of attack on the model. It can be vividly seen that the total translation and total rotation are higher for a high dynamic pressure (sea level) and higher angle of attack. (α = 8 • ). Figure 15 shows the convergence of the model for different element size.      Table 6 represent the computational time for all mesh sizes and number of panel divisions (boxes) analyzed. As expected, it can be seen that, for a smaller element size and higher number of panel divisions, the computational cost increases.

Conclusions
Detailed modelling of the aeroelastic behaviour is computationally expensive, therefore, it is required to develop a less complex model that can be used during the preliminary design stage. In this paper, an aeroelastic beam model was developed using Femap with Nx Nastran.The beam model is computationally less expensive; therefore, a beam model for the wing with winglet was derived. The static aeroelastic analysis was performed for different mesh sizes, panel division, and flight conditions. Four cases were analyzed, 1. α = 3 • at Sea Level with dynamic pressure = 11,348.4 Pa 2. α = 3 • at 7000 m with dynamic pressure = 4604.32 3. α = 8 • at Sea Level with dynamic pressure = 11,348.4 Pa 4. α = 8 • at 7000 m with dynamic pressure = 4604.32 and each of these cases was implemented on six models with different element size and panel division (Table 6). In total, 24 analysis were performed and, from all the analysis performed, it can be concluded that decreasing the mesh size and increasing the number of panel division increase the accuracy of the numerical results. However, as the number of elements and the number of panel division increases, the cost (time and memory) increases.
Out of the four cases analyzed, case 3 (sea level with α = 8 • ) has the higher total translation and higher stresses, as expected, due to larger angle of attack and higher dynamic pressure.
In this paper, only the static aeroelastic case was studied; however, since it is known that flutter is the most important instability, the dynamic aeroelastic analysis and optimization will be performed in future works.
All of the works here are performed using commercially available software; however, in future works an in built code will be developed and the numerical accuracy of this code will be compared with the results of commercially available software and with the experimental results. A prototype for the wind tunnel test will also be produced and experimental results will be obtained.