A Parametric Study on the Aeroelasticity of Flared Hinge Folding Wingtips

: This paper presents a parametric study on the aeroelasticity of cantilever wings equipped with Flared Hinge Folding Wingtips (FHFWTs). The ﬁnite element method is utilized to develop a computational, low-ﬁdelity aeroelastic model. The wing structure is modelled using Euler–Bernoulli beam elements, and unsteady Theodorsen’s aerodynamic strip Theory is used for aerodynamic load predictions. The PK method is used to estimate the aeroelastic boundaries. The model is validated using three rectangular, cantilever wings whose properties are available in literature. Then, a rectangular, cantilever wing is used to study the effect of folding wingtips on the aeroelastic response and stability boundaries. Two scenarios are considered for the aeroelastic analysis. In the ﬁrst scenario, the baseline, rectangular wing is split into inboard and outboard segments connected by a ﬂared hinge that allows the outboard segment to fold. In the second scenario, a folding wingtip is added to the baseline wing. For both scenarios, the inﬂuence of fold angle, hinge-line angle (ﬂare angle), hinge stiffness, tip mass and geometry are assessed. In addition, the load alleviation capability of FHFWT is evaluated when the wing encounters discrete (1-cosine) gusts. Finally, the hinge is assumed to exhibit cubic nonlinear behavior in torsion, and the effect of nonlinearity on the aeroelastic response is assessed and analyzed for three different cases.


Introduction
Folding wingtips have gained significant attention for various applications such as load alleviation, flight performance and storage [1]. A representative example for utilizing folding wingtips for flight performance (in-flight) and storage (on-ground) is the Boeing B777X aircraft, which employs a folding wingtip to be activated only on the ground during taxing to and from the gates, allowing the aircraft to operate from smaller airports and fit within gate limitations. During flight, the wingtip is unfolded to increase the wingspan, reducing the induced drag. Recently, significant research efforts, led by Airbus, have been dedicated toward utilizing folding wingtips for load alleviation [2]. The load alleviation capability associated with folding wingtips is mainly achieved by flaring/sweeping the hinge-line relative to the airflow as shown in Figure 1. Due to folding, the wingtip geometric incidence angle becomes negative due to hinge-line flare. In other words, an upward wingtip motion results in an aerodynamic downforce which can be used to alleviate loads and to ensure structural stability of the wing-wingtip system.
Since this paper focuses on the aeroelasticity of flared hinged folding wingtips for load alleviation, only the most recent studies utilizing them for load alleviation are highlighted and discussed. For more details on the aeroelasticity of folding wingtips (dihedral morphing), the reader is advised to refer to Ajaj et al. [2]. In terms of load alleviation, Wilson et al. [3] confirmed that a zero stiffness, flared hinge folding wingtip helps in gust and maneuver load alleviation with weight-saving opportunity for short-range aircraft. This concept might induce flutter if not eliminated by adding tip mass. They showed that the choice of hinge flare and location had a small effect on the root bending moment. Similarly, Castrichini et al. [4] utilized folding wingtips to reduce static and dynamic loads during flight for a transport aircraft. The wingtip was connected to the main wing using elastic hinges, and the effect of stiffness, wingtip weight, damping and hinge orientation on the response was assessed. The static analysis showed that the hinge orientation affected the response. For low values of wingtip mass and hinge spring stiffness (including zero stiffness), an increase in the hinge angle with respect to the freestream direction allowed better load alleviation capability. Furthermore, it was also determined that the flutter analysis showed that the low wingtip mass at 25 • hinge angle was beneficial for the aeroelastic stability. It was concluded that an increase in wingspan by 25%, using the folding wingtip, caused almost no increase in the load level for the cases considered. Castrichini et al. [5,6] took the idea of flared hinge folding wingtip further by considering a passive nonlinear hinge spring to connect the folding wingtip to the main wing. The nonlinear passive spring allowed wingtip deflection only when a load threshold was reached. The analysis showed that the load alleviation capabilities were influenced by the hinge moment threshold to release the wingtip and by the hinge damping value. Low hinge damping, coupled with low threshold hinge moment, allowed rapid deflection of the folding device due to trim loads and positive gusts. However, this approach caused an increment in the wing root bending moment. On the other hand, an increase in the hinge moment threshold caused higher wing root bending moment and delayed the wingtip rotation. Similarly, Cheung et al. [7] performed a series of low-speed steady and dynamic wind-tunnel tests using a flared hinge folding wingtip and found that such a device could provide gust load alleviation. The level of load alleviation varied with hinge spring stiffness and lifting condition, with the best scenario achieving a 56% reduction in peak loads. Castrichini et al. [8] studied the impact of the Semi Aeroelastic Hinge (SAH), which is a floating flared hinge folding wingtip, on the flight dynamics of a narrow-body transport aircraft. The integrated model of flight dynamics and aeroelasticity was built upon the simplified version of the formulation proposed by Saltari et al. [9]. The Doublet Lattice Method (DLM) was used to model the unsteady aerodynamic loads. The results showed that regardless of the 25% increment in wingspan, the aircraft equipped with the floating wingtip had the same dynamic response and handling qualities as the baseline aircraft (with no wingtip extension). They also concluded that the SAH could be used to alleviate the increase in roll damping associated with the longer wingspan, and it can be used as a load alleviation device as well. Similarly, Ajaj [1] studied the impact of flared hinge folding wingtips on the flight mechanics characteristics of a narrow-body transport aircraft using quasi-static aeroelastic tool consisting of FE structural model and tornado vortex lattice method. He concluded that both the hinge-line/flare angle and wingtip span had significant effect on the static flight loads causing large reductions in the root bending moment and root shear force for a rigid wing. However, the effectiveness of flared hinge folding wingtips dropped sharply when accounting for structural flexibility. Flared hinge folding wingtips had little effect on the longitudinal flight dynamic modes regardless of the flare angle, fold angle and wingtip size. They slightly improved the damping for the Phugoid mode but had noticeable impact on the lateral-directional flight dynamic modes, especially the spiral mode. Similarly, the influence of structural flexibility was of minor importance except on the spiral mode causing large reduction in its time constant. Finally, the structural flexibility reduced the effectiveness of the ailerons and the roll rate of the aircraft, but this disappeared as the fold and flare angles increased. Wilson et al. [10] advanced the AlbatrossONE SAH small-scale demonstrator aircraft to represent a fullscale transport aircraft. The AlbatrossONE was a scaled-down short-range aircraft with a fuselage based on Airbus A321 as shown in Figure 2. They conducted extensive windtunnel and flight tests. It should be noted that in their study, the dynamic scaling for either handling qualities and/or aeroelasticity was not performed. The flight testing showed that the wingtips were both statically and dynamically stable. The effect of wing load alleviation from the free wingtips was confirmed by comparing the strain gauge measurements between two flights. They concluded that a nonlinear relationship existed between the angle of attack and the coasting (free folding) angle of the wingtip due to effective sideslip of the hinge flare angle. Finally, a near-linear variation of wingtip flapping frequency with the airspeed existed, and this was confirmed by wind-tunnel testing. As stressed above, the aeroelasticity of flared-hinge folding wingtips has been studied from a computational perspective; however, all the studies utilized high-fidelity models which are computational expensive to conduct a comprehensive parametric study and account for nonlinearities. This paper aims to fill this gap and establish a low-fidelity aeroelastic model that facilitates fast and robust aeroelastic assessments of flared-hinge folding wingtips. Due to its low-fidelity, the model allows a deeper understanding of the various aeroelastic phenomena and allows conducting a wide range of parametric studies to determine the influence of various design parameters on the behavior of the system. Section 2 of the paper presents the FE aeroelastic toolbox development and its validation. Section 3 presents the parametric study on the impact of folding wingtips on the aeroelastic behavior of the wing. In Section 4, nonlinear aeroelastic analysis is performed assuming the hinge to exhibit cubic structural nonlinearity in torsion. Finally, Section 5 presents the conclusions and main findings of this paper.

Modelling
A flared hinge folding wingtip experiences a change in its geometric incidence angle (∆θ FWT ) as when it folds according to [1]: where Λ hinge is the hinge-line angle (or flare angle) and Γ FWT is the fold angle of the wingtip. The hinge flare angle is defined as the angle between the longitudinal axis of the aircraft (or the freestream velocity vector at zero sideslip angle) and the hinge rotation axis as shown in Figure 1a. According to the coordinate system shown in Figure 3, a negative fold angle occurs when the tip folds upward. Similarly, a negative flare angle is when the leading edge of the hinge is outboard of its trailing edge. It should be noted that for load alleviation purposes, only negative flare and fold angles are considered.

Structural Dynamics Solver
A finite element (FE) structural solver is developed in Matlab TM . The lifting surface is discretized into Euler-Bernoulli beam elements as shown in Figure 3. The elastic potential energy (U) of a beam element, including the effects of the geometric warping, can be expressed as [11]: where EI is the bending rigidity, GJ is the torsional rigitiy, C w is the warping constant, E is the Young's modulus, θ is the twist angle, w is the bending deflection, y is the position along the beam span and l is the wing length. Mei [11] provided expressions for the warping constant for different beam cross-sections. Similarly, the kinetic energy (T) of a beam element can be expressed as: where m is the mass per unit length, t is the time, r α is the radius of gyration and e m is the distance between shear center and center of gravity.
To obtain the stiffness and mass matrices for each beam element, Hermite's cubic shape functions, which satisfy the boundary conditions, are used to model the bending deflection: where N b represents the bending shape functions (whose expressions are detailed in the Appendix A), and it can be written as: and w e represents the bending displacements (w 1 and w 2 ) and rotations (φ 1 and φ 2 ) at the ends of a beam element. It can be expressed as: On the other hand, the following linear shape functions are used to model the torsion deflection of the beam element as: where N T t represents the torsion shape functions (whose expressions are given in the Appendix A) and it can be expressed as: and θ e represents the twist angles at the ends of the beam element. It can be expressed as: Using transformation matrices, a beam element can be rotated from its local coordinate system to the wing/global coordinate system. This approach allows the structural solver to handle different wing geometries with different twist, sweep, and dihedral distributions. Since a linear shape function is used for torsion, the effect of torque-induced warping cannot be captured and is ignored. It should be noted that the solver can easily account for warping by using higher order shape functions.

Modelling the Flared Hinge
As stressed earlier, the lifting surface is discretized into elements modelled as Euler-Bernoulli beams. However due to hinge flaring/sweeping, the two elements inboard and outboard of the hinge have swept edges as shown in Figure 4. Therefore, they cannot be modelled as standard beam elements and developing analytical expressions for their stiffness matrices proved to be a very difficult task and beyond the scope of this paper. Consequently, a condensation-based computational procedure is adopted here in order to compute numerically the 6 × 6 stiffness matrices associated with these two nonstandard elements [12]. We believe that the choice of 3 degrees of freedom per node is suitable to capture properly the aeroelastic response of each nonstandard beam. These are the plunge displacement, w, in the z-direction, the bending rotation, φ, about the x-axis; and the torsional rotation, θ, about the y-axis as shown in Figure 3. The present procedure consists of modelling each of the nonstandard elements shown in Figure 4 by means of full 3D quadratic shell type elements. Furthermore, rectangular beam elements are added to the leading and trailing edges of the nonstandard elements in order to capture the asymmetrical geometry of the airfoil with the aid of an equivalent cross-section. The thickness of the shell elements and the cross-section of the beam elements are chosen such that the total sum of their cross-section areas and moments of inertia coincide with those values of the actual airfoil geometry. A typical finite element mesh contains 2590 shell elements, 370 beam elements and 8169 nodes. This results in a total number of 49,014 degrees of freedom and, therefore, in an associated 49014 × 49014 stiffness matrix.
As the central idea of this approach is to condensate the large number of degrees of freedom of the above full 3D shell structure into a 6 × 6 stiffness matrix of a mechanically equivalent beam, six cases of prescribed nodal displacements and rotations in the 3D model are analyzed. The corresponding values for each degree of freedom in each loading case are listed in Table 1.
The above-prescribed displacements/rotations are applied to the nodes 1 and 2 located along the elastic axis, as shown in Figure 3, leaving the rest of the degrees of freedom free to move. Nevertheless, all the degrees of freedom of those nodes contained in the cross-sections located at nodes 1 and 2 of the shell model are coupled to reproduce rigid body motion. This is compatible with the structural behavior of the other standard Euler-Bernoulli beam elements when modelling the full wing geometry.
For each displacement/rotation vector {w 1 φ 1 θ 1 w 2 φ 2 θ 2 } T applied to structure, the following standard relationship can be established: is the corresponding reaction forces and moments vector, and [K] is a 6 × 6 stiffness matrix which can be associated with that of the equivalent nonstandard beam element. For the first loading case, the following expression can be written: where the terms K ij are the stiffness components of [K], with i and j referring to the numbers of degree of freedom from 1 to 6. By observing the above expression, we can conclude that the reaction vector obtained from the solution of the problem is equal to the first column of the equivalent nonstandard beam element stiffness matrix. In a similar fashion, the remaining five columns of [K] can be obtained when the other five loading cases are analyzed.
It should be noted that the numerical condensation procedure described previously is carried out only twice, just to produce the stiffness matrices of the two nonstandard beam elements shown in Figure 4. Therefore, the CPU times associated with the computation of these two matrices are low when compared to the solution of the aeroelastic system of the full wing. In addition, the condensation procedure is applied to reduce the large number of degrees of freedom existing in a complex 3D quadratic shell element model containing swept edges into the 6 degrees of freedom present in the beam element shown in Figure 3. This numerical strategy is implemented in the commercial finite element analysis software Ansys ® . It should be noted that the stiffness matrices of these nonstandard beam elements contain bending-torsion coupling terms that do not exist in the matrices of the other standard Euler-Bernoulli beam elements. The magnitude and sign of these coupling terms depend on the flare angle and the dimensions of the element. The hinge region (shown in Figure 4) is modelled as a swept standard beam element whose sweep angle is equal to the flare angle of the hinge.
For the sake of illustration and to evaluate the magnitude of the coupling terms associated with these nonstandard beam elements, the flared hinge is added to Goland wing whose properties are listed in Table 1. Three scenarios, corresponding to flare angles of 0 • , 35 • and 70 • are considered. The fold angle is set to 0 • for all the scenarios. Figure 5 shows the stiffness matrices of the inboard beam element (assuming 3 degrees of freedom per node corresponding to bending and torsion) for the three different scenarios. It should be noted that the transposition of the degrees of freedom vector is The beam element length is set to 1 m in the three scenarios. It can be seen from Figure 5 that the magnitudes of the bending-torsion coupling terms are zero at 0 • flare angle, and they increase as the flare angle increases. It should be noted that the magnitude of these coupling terms is much smaller when compared to the other terms even at a flare angle of 70 • . Due to their small magnitude, it is anticipated that these coupling terms do not have significant influence on the aeroelasticity of a wing equipped with FHFWT.

Aerodynamic Solver
The aerodynamic solver used here is based on Theoderosen's unsteady aerodynamic strip method. Theodorsen's unsteady aerodynamic theory, developed in 1935, has a circulatory component to account for the effect of the wake on the airfoil, and it contains the main damping and stiffness terms and a noncirculatory component to account for the acceleration of the fluid surrounding the airfoil [13]. The unsteady lift and pitching moment equation (around elastic axis) per unit span are given as: where ρ is the air density, b is the semichord,â is the normalized pitch axis location with respect to half chord, V is the true airspeed, k is the reduced frequency and C(k) is the Theodorsen's function, which is used to model the changes in amplitude and phase of the unsteady forcing of harmonically oscillating airfoils [13]. The function is essentially a Fourier transform of the Wagner function, and it acts as a frequency domain filter for oscillating motion input to give a reduced frequency dependent output [14,15]. Theodorsen's theory ignores the three-dimensional effects; however, it allows modifying the lift curve slope to account for the effect of aspect ratio and sweep. The aerodynamic stiffness and damping matrices are then obtained using the structural shape functions utilized for structural modelling. Assuming that the wing is subjected to linear harmonic forcing (at the flutter speed), the response in plunge and pitch can be expressed for each coordinate respectively as: where w o and θ o represent the amplitudes of oscillation in plunge and pitch, respectively. The lift and pitching moment per unit span can be expressed as: where L w , L . w , L θ , L . θ , M w , M . w , M θ and M . θ are the nondimensional oscillatory aerodynamic derivatives. Their definitions are given in Wright and Cooper [15]. The aerodynamic loads on a strip can be rearranged in matrix format as follows: where and To obtain the aerodynamic loads on a strip element, the work conducted by these loads on a strip can be expressed as: where N e represents the shape functions in bending and torsion and can be written as: This allows one to obtain the aerodynamic damping and stiffness matrices as follows: As stressed earlier, the flared hinge is modelled as a swept beam from a structural perspective; however, from an aerodynamic perspective, its sweep is the same as that of the wing. In the case of Goland wing (a rectangular wing), the hinge has zero aerodynamic sweep.

Aeroelastic System
The lifting surface is split into a number of elements. For each element, the mass and stiffness matrices are computed and transformed to the global coordinate system. Similarly, the aerodynamic stiffness and damping matrices are computed for each element and transformed to the global coordinate system of the wing using direct cosine matrix. The matrices are then assembled and arranged such that the aeroelastic system can be expressed as:Â ..
whereÂ,D andÊ are the global structural inertia, global structural damping (set to zero here) and global structural stiffness, respectively. On the other hand, X represents the generalized coordinates. The aerodynamic stiffness Ĉ and damping B matrices depend on the reduced frequency (k). Thus, a frequency matching method is necessary to determine the aeroelastic boundaries. The PK (British Method) developed by Hassig [16] is used. The main assumption of the PK method is that the aerodynamic behavior is dependent upon a harmonic response. This is true at the flutter condition but not true below or above the flutter speed. More details on the PK method can be found in Wright and Cooper [15].

Validation
Three rectangular, cantilever wings are considered to validate the FE aeroelastic toolbox. The properties of these wings are listed in Table 2 [17]. The validation results are summarized in Table 3. It should be noted that 6 beam elements are used for the Goland wing, 12 for the HALE wing and 12 for the representative wing. It can be seen that the aeroelastic model developed and used here is quite accurate in predicting the aeroelastic boundaries when compared to similar methods and tools available in the literature.

Parametric Study: Linear Aeroelasticity
The Goland wing is used here for further aeroelastic analysis where two scenarios are considered. In the first scenario, the baseline Goland wing is split into inboard and outboard segments connected by a flared hinge and the outboard segment is allowed to fold. This scenario is studied in Section 3.1. In the second scenario, a folding wingtip is added to the baseline wing. This is presented in Section 3.2. The main difference between Sections 3.1 and 3.2 is that in Section 3.1, the baseline wing is split into two segments but the span of the wing does not change. However, in Section 3.2, a wingtip is added to the baseline wing so the effective span of the wing increases. It should be noted that the divergence speed, flutter speed and flutter frequency, in some of the figures below, are normalized using their corresponding values of the baseline wing (listed in Table 2).

Folding the Baseline Wing
The baseline wing is split into two segments: a fixed inboard segment and a folding outboard segment. At zero-fold angle, the semispan of the wing is 6.096 m. The impact of folding the outboard segment is studied for different flare angles. Two different lengths for the outboard segment, corresponding to 1.016 m (16.66% of the semispan) and 2.032 m (33.33% of the semispan).

Flare Angle
The variations of flutter speed, flutter frequency and divergence speed with fold angle for different flare angles and outboard segment lengths are shown in Figures 6 and 7. It can be deduced from Figures 6 and 7 that increasing the fold angle always increases the divergence speed regardless of the flare angle and outboard segment length. However, larger flare angles result in a higher increase in divergence speed. The flutter frequency has the same trend as the divergence speed.

Quasi-Steady Response to Discrete Gusts
To evaluate the load alleviation capability of FHFWT, it is essential to determine the variations in root bending moment and shear force when the wing encounters gusts. For the sake of illustration, discrete, 1-cosine gusts are considered here. The 1-cosine gust profile can be defined, according to FAR Part 25, Section 25.341, as: where H, the gust gradient, is the distance parallel to the airplane's flight path for the gust to reach its peak velocity, and it varies from 9.144 to 106.7 m. U ds is the design gust velocity, and it can be expressed as: where U re f , the reference gust velocity, has a magnitude of 17.07 m/s at sea level and reduces linearly from 17.07 to 13.4 m/s EAS at 15,000 feet. F g is the flight profile alleviation factor, and it is set to 1. The airspeed and the angle of attack are set to 80 m/s and 0 • respectively. Figure 8 shows the variations of the maximum root bending moment (RBM) and maximum root shear force (RSHF) for outboard segment lengths of 1.016 and 2.032 m. It should be noted that in this section, quasi-steady aerodynamic loads are used where the gust velocity is assumed to result in an instantaneous change in the angle of attack. It can be clearly seen from Figure 8 that for the positive gust cases, folding the wingtip reduces the RBM and RSHF regardless of the flare angle due to the reduction in the effective wingspan (moment arm). However, the combination of fold and flare angles offers significant reductions in the root loads for all the gust length considered. It should be noted that larger outboard segment length offers more reductions in root loads. For negative gusts, the opposite is true where the FHFWT increases the roots loads mainly due to the aerodynamic downforce generated due the hinge flare angle.

Hinge Stiffness
To determine the influence of hinge stiffness on the aeroelastic boundaries, the spanwise bending rigidity of the hinge (beam element connecting the main wing and the wingtip) is varied from 10 to 10 8 Nm whilst keeping its torsional rigidity as the baseline wing. Then, its torsional rigidity is varied from 160 to 10 8 Nm whilst keeping its bending rigidity as the baseline wing. The outboard segment length is fixed at 1.016 m. The variations of flutter speed, flutter frequency and divergence speed are shown in Figures 9 and 10.  Figure 9 shows that divergence speed is not very sensitive to the bending rigidity of the hinge. Increasing the bending rigidity has no effect on the divergence speed except for a flare angle of 15 • and fold angle of 40 • where the divergence speed increases slightly from 277 at 10 Nm 2 to 285 m/s at 200 Nm 2 , and then it remains constant even when the bending rigidity increases further. On the contrary, flutter speed and frequency are very sensitive to the bending rigidity of the hinge especially for values between 10 2 and 10 4 Nm 2 . In this region, they increase sharply and then drop suddenly regardless of the flare and fold angles. This sudden change in the flutter behavior is associated with a change in the flutter mode from the 2nd bending to 1st torsion. For example, consider the 0 • fold angle and 15 • flare angle case. At a bending rigidity of 10 2 Nm 2 , the flutter speed is 121 m/s, and it rises to 189 m/s at 310 Nm 2 . Then, it drops sharply to reach 138 m/s at 10 4 Nm 2 and remains almost constant for higher values of bending rigidity. Similarly, the flutter frequency is 81 rad/s at 10 2 Nm 2 , and it drops sharply to reach 62 rad/s at 310 Nm 2 . Then, it increases again to reach 69 rad/s at 10 4 Nm 2 and remains almost constant for higher values of bending rigidity. It should be noted that the curve of the 0 • fold angle and 15 • flare angle case coincides with the 0 • fold angle and 0 • flare angle case. It can be deduced from Figure 9 that increasing the fold angle alone causes the flutter speed to peak at higher values of bending rigidity. On the other hand, increasing the flare angle alone tends to increase the peak values without affecting the values of bending rigidity at which the peaks occur. It should be noted that the effect of Coriolis force can become important when the wingtip achieves a high angular velocity, and this is most likely to happen at low bending stiffness and zero flare angles (no aerodynamic stiffness). However, this effect was not modelled here because in practice there is always a flare angle which prevents the tip from reaching high angular velocities even at very low bending stiffness (due to the presence of aerodynamic stiffness). In addition, this paper focuses on the onset of instability and this allows ignoring the effect of Coriolis force. Figure 10 shows that for 80 • fold angle, the divergence speed is not sensitive to torsional rigidity of the hinge regardless of the flare angle. For 40 • fold angle, divergence speed is very sensitive to torsional rigidity at 0 • flare angle and insensitive at 25 • flare angle. For 0 • fold angle, divergence speed is very sensitive to torsional rigidity regardless of the flare angle. For 0 • and 40 • fold angle, the flutter speed drops as the torsional rigidity increases and then it starts increasing steadily until the torsional rigidity reaches 10 4 Nm 2 regardless of the flare angle. For 80 • fold angle and 0 • flare angle, the flutter speed is not sensitive to torsional rigidity, whilst for a 25 • flare angle, the flutter speed reduces as the torsional rigidity increases until 490 Nm 2 where a further increase in the torsional rigidity does not affect the flutter speed. The flutter frequency tends to increase with the torsional rigidity until 10 4 Nm 2 , after which it remains constant. The only exception is the 80 • fold angle at 0 • flare angle where the flutter frequency is insensitive to the torsional rigidity of the hinge. It also worth noting that at 250 Nm 2 torsional rigidity, the flutter speed for the 25 • flare angle and 80 • fold angle case becomes higher than that for the 0 • flare angle and 80 • fold angle case. The opposite is true for the flutter frequency. Regardless of the torsional rigidity, the flutter mode does not change and remains the first torsion mode. It can be deduced from this subsection that the torsional rigidity of the hinge affects static and dynamic aeroelastic boundaries, whereas its bending rigidity mainly affects the dynamic aeroelastic boundaries. In addition, the variation in torsional rigidity does not change the flutter mode (first torsion), whereas the variation in bending rigidity results in a change in the flutter mode (from second bending to first torsion).

Adding Folding Wingtip
In this subsection, a rectangular wingtip device is added to each side of the baseline Goland wing. The span, mass per unit length, sweep angle, flare angle, and fold angle associated with the wingtip are varied, and their effects on the aeroelastic boundaries are studied.

Wingtip Size
The fold angle is varied from 0 • to 60 • , and four flare angles, corresponding to 0 • , 15 • , 30 • and 45 • are considered. Two wingtip spans, corresponding to 1.016 and 2.032 m, are studied. The variation of flutter speed, flutter frequency and divergence speed are shown in Figures 11 and 12.

Wingtip Mass
To determine the impact of wingtip's mass on the aeroelastic stability of the wing, three mass per unit length m wingtip are considered. These values correspond to 17.86, 35.71 and 71.42 kg/m. The wingtip span is set to 1.016 m. It can be seen, from Figure 13, that as the mass per unit span of the wingtip increases, the flutter speed becomes more sensitive to the flare angle but less sensitive to the fold angle, especially for 25 • flare angle. In fact, at large values of m wingtip and flare angle, the flutter speed decreases with increasing the fold angle. On the other hand, flutter frequency becomes less sensitive to flare angle and fold angles as m wingtip increases. It should be noted variations in m wingtip do not affect divergence speed, and therefore, it is not shown here.

Wingtip Sweep
The effect of wingtip sweep angle, shown in Figure 14, is assessed here. It should be noted that for swept wings, Theodorsen's unsteady formulation detailed above should be adjusted. There are two popular formulations to account for sweep: the streamwise formulation and the chordwise formulation. The streamwise formulation is used here due to its simplicity, and this is justified by assuming that the wing ribs are pointing in the streamwise direction preventing streamwise camber change as the wing deforms. Different aft sweep angles, ranging from 0 • to 25 • , are studied as shown in Figure 15. The wingtip span is set to 1.016 m.    Figure 15 shows that divergence speed increases with increasing the aft sweep angle of the wingtip. It should be noted that combining flare angle with aft sweep angle tends to increase divergence speed significantly. This is mainly because the flared hinge is modelled as a beam with geometric sweep. At 0 • aft sweep, the flutter speed increases with the fold angle regardless of the flare angle. However, the rate of increase is much higher with 0 • flare angle. The same variation trend exists for 12.5 • aft sweep. On the contrary, for 25 • aft sweep, flutter speed increases with fold angle at 0 • flare angle, and it reduces at 25 • flare angle. Finally, flutter frequency is not very sensitive to aft sweep or to flare angles. However, it is quite sensitive to the fold angle, and it increases significantly with increasing the fold angle regardless of the sweep and flare angles.

Nonlinear Aeroelasticity
Nonlinearities in aircraft are unavoidable and arise in the structure, aerodynamics and the control system [23]. In this paper, the focus is on structural nonlinearities, and these may appear as a result of worn hinges, loose control linkages, nonlinear material behavior, cracks, notches and/or fatigue. Structural nonlinearities can be classified into either distributed or concentrated nonlinearities based on the region of their action. Distributed nonlinearities arise from general deformations of the entire structure, whereas concentrated structural nonlinearities act locally and are commonly found in control mechanisms or connecting parts. Concentrated structural nonlinearities can be approximated by cubic, bilinear and hysteresis structural nonlinearities.

Nonlinear Hinge
The hinge (beam) connecting the wingtip to the main wing is assumed to exhibit a cubic nonlinearity in the elastic torsional moment. The cubic nonlinearity in the elastic torsional moment can be represented by a quadratic nonlinearity in the torsional rigidity of the hinge as: GJ hinge (∆θ) = a 1 + a 3 ∆θ 2 = a 1 1 + γ∆θ 2 (28) where a 1 and a 3 are constants, γ = a 3 /a 1 and ∆θ represent the change in twist angle across the length of the hinge. Positive γ indicates a hardening hinge, whilst negative γ indicates a softening hinge. The above nonlinearity model assumes the hinge to exhibit antisymmetric loading and unloading behavior. The nonlinearity represented in Equation (28) directly affects the structural stiffness matrix of the wing.

Limit Cycle Oscillations
The aeroelastic equations of the wing, coupled with the nonlinear hinge model, are solved numerically in Matlab TM using the ODE23 solver. The ODE23 is a single-step solver based on an explicit Runge-Kutta formulation [24,25]. As in Section 3.1, the baseline Goland wing is split into an inboard and an outboard segment. The length of the outboard segment is set to 1.016 m. The hinge properties are a 1 = 987 Nm 2 and γ = 500 (hardening).
To determine the influence of hinge nonlinearity on the aeroelastic response of the wing, three cases are studied here: The airspeed is set to 121 m/s, which is above the linear flutter speed for the three cases (when ∆θ = 0), whilst the angle of attack is fixed at 5 • . Figure 16 shows the time history and phase response of the wingtip for the three cases (blue curve for Case 1, red curve for Case 2, and black curve for Case 3). It can be clearly seen that at the start, the wing diverges (oscillations grow), causing an increase in ∆θ, and therefore the torsional rigidity of the hinge rises according to Equation (28). This increase in torsional rigidity prevents the amplitude of oscillation from growing further and induces a limit cycle oscillation (bounded harmonic oscillation). The amplitude of oscillation is higher for Case 1 (where the fold angle is zero). This is mainly because the linear flutter speed of Case 1 is lower (108.7 m/s) than that of Case 2 (113.8 m/s) which implies that larger pitch oscillation (and therefore higher ∆θ) is required to achieve higher increase in torsional rigidity and therefore stop the diverging behavior and achieve an LCO. For Case 3, the flare angle reduces the linear flutter speed to 112.4 m/s when compared to Case 2. This implies that larger amplitude of oscillations is required to achieve the required stiffness. For all the cases, the LCO commences almost after 1 s from the start. In summary, the nonlinear hinge allows the wing to operate at speeds above its linear flutter speed without catastrophic failure. Increasing the fold angle reduces the amplitude of the LCO, whilst increasing the flare angle increases the amplitude of the LCO.

Conclusions
A low-fidelity finite element aeroelastic model was developed and utilized to study the aeroelasticity of wings with flared hinge folding wingtips. The wing structure was modelled using Euler-Bernoulli beam whilst unsteady Theodorsen strip theory was used for aerodynamic prediction. The toolbox was validated using three rectangular, cantilever wings. Then, the flared hinge folding capability was integrated on Goland wing, and two scenarios were studied. In the first scenario, Goland wing was split into an inboard segment, and an outboard segment connected by flared hinge where the outboard was capable of folding. In the second scenario, a flared hinge folding wingtip was added to Goland wing. For each scenario, the influence of flare angle, fold angle, wingtip mass, wingtip sweep, and hinge stiffness on flutter speed, flutter frequency and divergence speed were assessed. For both scenarios, the flared hinge folding wingtip expands the aeroelastic envelope for the majority of fold and flare angles. The variation of root loads with folding wingtips for 1-cosine gusts was studied, and it showed promising reductions in root bending moment and shear force due to the flare angle. Finally, the flared hinge was assumed to exhibit cubic nonlinearity in torsion with hardening effect, and the wing (from the first scenario) was set at an airspeed above the linear flutter speed. The nonlinearity in the torsion behavior of the hinge caused limit cycle oscillations allowing the wing to fly above its critical speed without diverging oscillations.  Data Availability Statement: Some or all data, models or code that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.
first derivatives with respect to time .. second derivatives with respect to time Subscripts e ends of the beam element o amplitude