Nonlinear aeroelastic stability analysis of a two-stage axially moving telescopic wing by using fully intrinsic equations

During the process of span extension for an aircraft wing equipped with a telescopic morphing mechanism, the wing aspect ratio increases, and hence, the geometrical nonlinearities might become more significant. In this regard, this paper aims to investigate the effect of structural nonlinearity on the aeroelasticity of span morphing wings using the exact fully intrinsic equations for the first time. Furthermore, the effects of various parameters such as thrust force, engine location, chord size, flight altitude, initial angle of attack, and overlapping mass on the aeroelasticity of the wing are studied. The applied aerodynamic loads in an incompressible flow regime are determined using Peters’ unsteady aerodynamic model. In order to check the stability of the system, first the resulting nonlinear partial differential equations are discretized by using the central finite difference method and then linearized about the static equilibrium. Finally, by obtaining the eigenvalues of the linearized system, the stability of the wing is evaluated. It is observed that by using the fully intrinsic equations, the instability of the axially moving telescopic wing can be determined more accurately. Moreover, the results show that the morphing length and overlapping mass have significant effects on the aeroelastic stability of the telescopic wing.


Introduction
Morphing wings have received a great deal of attention over the past decade for their promising performance. Among all possible morphing concepts, the wing extension concept has been shown to be an effective way to improve aircraft range and endurance 1 which can be considered in the design of UAV wings. However, as a result of changes in the wing dimensions, the dynamic characteristics of the system can be altered. 2 Therefore, understanding the dynamic behavior of morphing wings is a critical step in the design process of such structures.
Due to their high aspect ratio, the structural dynamics of UAV wings can be modeled using beam models. 3,4 Many researchers have analyzed the dynamic behavior of axially moving beams. Wang and Wei 5 studied the vibration of a robot arm modeled by a moving slender prismatic beam. In their study, it was shown that increasing or decreasing the length of the flexible arm has destabilizing or stabilizing effects on the arm vibrations. Stylianou and Tabarrok 6 presented a finite element analysis of an axially moving beam. This study was then continued by considering the effects of physical damping, tip mass, tip support, and wall flexibility on the stability characteristics of the aforementioned beam using the eigenvalue analysis. 7 Raftoyiannis and Michaltsos 8 employed a modal superposition technique for dynamic analysis of telescopic cranes' boom based on a continuum approach. Chang et al. 9 used a finite element method to derive the equations of motion of an axially moving beam based on the Rayleigh beam theory. In this study, the stability of the beam, with constant extension speed, was sought using the eigenvalue analysis. Furthermore, the Floquet theory was employed to investigate the stability of the beam with a periodical back-and-forth motion. Duan et al. 10 studied the dynamic response of an axially moving nested beam theoretically and experimentally. The equations of motion were obtained using D'Alembert's principle and a good agreement between the experimental and numerical results was observed. The effects of a moving mass on the dynamic behavior of a two-stage telescopic mechanism, used in truss structures of a bridge inspection vehicle, was considered by Sui et al. 11 The structural dynamics of the telescopic mechanism was modeled using Euler-Bernoulli beam theory. This study was then continued to investigate the dynamic behavior of a 2-DOF telescopic mechanism. 12 As it was mentioned above, the structural dynamics of span morphing beam is dependent on its length, and hence, the aeroelastic characteristics of telescopic wings can also be affected. Huang and Qiu 13 studied the effects of span morphing velocities on the aeroelastic stability of a single variable-span with uniformity assumption (uniform time-invariant parameters). The established aeroelastic model was based on Euler-Bernoulli beam theory and the unsteady vortex lattice aerodynamic theory. It was shown that using a span morphing mechanism might improve the aeroelastic performance of the wing. Huang et al. 14 studied the effects of rigid-body motions on the aeroelastic response of span morphing wings. They combined Euler-Bernoulli beam theory with an unsteady strip aerodynamic theory and showed that the quasi-static stability of the morphing wing is dependent on the fuselage flexibility. Ajaj and Friswell 15 investigated the sensitivity of the flutter speed of a single variable-span morphing wing for various system parameters and morphing velocities with uniformity assumption. They combined wing bending and torsional equations with Theodorsen's unsteady aerodynamic model and showed that the morphing speed affects the aeroelastic stability of the wing and needs to be taken into account. The quasistatic aeroelastic behavior of telescopic, multi-segment, stepped, span morphing wings was studied by Ajaj et al. 16 Euler-Bernoulli beam theory and Theodorsen's unsteady aerodynamic model were combined to form the aeroelastic equations. The presented results showed that this mechanism can be used as a means for wing flutter suppression. It must be noted that in their study, the effect of morphing speed was ignored.
Although most of the previous studies concerned with the dynamics of axially moving beams were mainly focused on the linear behavior of the beam, several researchers developed nonlinear models to study these structures. Park et al. 17 determined the dynamic behavior of an axially moving beam modeled using the nonlinear von Karman strain theory. The results showed that the response of the system for both linear and nonlinear conditions was consistent. Furthermore, they showed that depending on the morphing acceleration and Young's modulus values, the differences between linear and nonlinear solutions might increase. Zhang et al. 18 investigated the nonlinear dynamic behavior of a single deploying-and-retracting wing modeled by a cantilever laminated composite shell in supersonic airflow. This study was then extended to subsonic airflow by combining the von Karman theory with Kutta-Joukowski lift theorem. 19 They characterized the effect of extension velocity on the nonlinear dynamic behavior and stability of the wing. Lu et al. 20 investigated the effects of piezoelectric material on the nonlinear dynamic behavior of a deploying laminated composite plate. They showed that the nonlinear vibration of the deploying cantilevered laminate can be suppressed by a suitable voltage and polarity. Guo et al. 21 studied the nonlinear vibration of a Z-shaped plate with inner resonance experimentally and numerically.
As understanding the nonlinear behavior of morphing wings is a key in their design, the present study is aimed at developing a novel nonlinear model to investigate the nonlinear dynamics of such structures in more detail. There are different nonlinear theories of beam available in the literature. Some theories are called displacement-based theories which are formulated only using displacement variables. 22,23 Some other beam theories are referred to as mixed formulation in which the rotation and displacement variables are related to the generalized velocity and strain measures by the kinematic partial differential equations. 24 Apart from these two categories, there is a third category of beam theories which are stress 24 or strain 25 based formulation. The present study is based on the fully intrinsic beam theory which was initially developed by Hegemier and Nair 26 and then revised by Hodges. 27 The important advantages of this theory in comparison with other structural beam theories are complete modeling without simplifying assumptions in large deformations, low-order nonlinearities, and thus less complexity. This beam theory has been used by several researchers to investigate the aeroelastic stability of aircraft wings. For example, Sotoudeh and Hodges 28 studied the effects of the joint position and sweep angle on the static deformation and dynamic stability of an aircraft wing subjected to a follower force. They showed that the static and dynamic stability of the joint-wing is completely different from conventional wings. Mardanpour et al. 29 Investigated the effect of engine placement and wing sweep on the flutter characteristics of the wing. This study was then continued by considering the effects of gust on the dynamic behavior of the wing. 30 Moravej Barzani et al. 31 studied the aeroelastic stability of swept wings using exact beam formulation. The results showed that by utilizing the geometrically exact fully intrinsic beam equations, the instability of the swept wings can be determined more accurately. Amoozgar et al. 32 investigated the effect of bend-twist elastic coupling and wing taper ratio in combination with the pretwist angle on the aeroelastic stability of the wing. In another work, Amoozgar et al. 33 studied the effect of the out-of-plane curvature on the flutter velocity and flutter frequency of the wing by using the geometrically exact beam equation. They highlighted the importance of the initial curvature and the length of the curved segment on wing aeroelastic stability.
In the present study, the nonlinear aeroelastic behavior of a two-stage span morphing wing is investigated quasistatically using the geometrically exact fully intrinsic beam model. 27 The important advantages of the geometrically exact fully intrinsic beam equation in comparison with other structural beam equations are complete modeling without simplifying assumptions in large deformations, low-order nonlinearities, and thus less complexity. The aerodynamic loads on the wing are determined using Peters' unsteady aerodynamic model. 34 The governing aeroelastic equations are discretized using a central finite difference method, 27 and the stability of the wing is determined using the eigenvalues of the linearized equations about the nonlinear steady-state condition. It is noted that in none of the previous studies, the geometrically exact fully intrinsic beam equation was used for the aeroelastic analysis of telescopic morphing wings, and hence, this is the main novelty and the main purpose of this study. Moreover, investigation of the effects of some important parameters such as thrust force, engine location, chord size, and flight altitude in the presence of angle of attack and overlapping mass on the aeroelasticity of the wing is another goal of this article.

Governing equations
As shown in Figure 1, a two-stage beam model is employed to simulate the structural dynamics of an axially moving telescopic wing. The wing has an axially moving speed (morphing speed) of _ κ, and the length of the fixed and moving parts are denoted by l 0 and l m , respectively. Also, as shown in Figure 2, X denotes the engine location in the x-direction, and p is the thrust force.
The following simplifications are considered in driving the equations of motion: 1. It is assumed that the morphing speed is very slow, and hence, a quasi-static condition is considered. 2. The amount of mass change of the fixed part is proportional to the size of the change in its length.
3. The engine's location is always on the fixed part without any offset from the reference axis. 4. Effects of engine mass are ignored in the calculations.
The total time-dependent length, l, can be expressed as As it was mentioned earlier, in this study, the geometrically exact fully intrinsic beam equations are used to simulate the nonlinear dynamic behavior of the wing. The fully intrinsic equations can be written as ( 27 ) F where ð Þ 0 denotes the derivative with respect to the beam reference line and _ ð Þ denotes the absolute time derivative, M B and F B are the internal moment and force measures, H B and P B are the sectional angular and linear momenta, V B and V B are the angular and linear velocity measures, and κ γ denote the moment and force strain measures. Furthermore, Κ B ¼ k b þ κ is the curvature vector in which k b is the initial twist and curvature of the beam. Also m B and f B are the external moments and forces such as aerodynamic moments and forces.
The generalized strains (γ, κ) and the generalized forces (F B , M B ) are related together using the crosssection stiffness matrix as follows where S, R, and T are the cross-sectional flexibility matrices. Furthermore, the generalized velocities (V , V) and the generalized moments (P, H) can be converted to each other using the cross-sectional inertia matrix as where Δ is the identity matrix, ξ is the cross-sectional mass centroid offset from the beam reference axis, μ is the mass per unit length, and I is the inertia matrix per unit length.
Wing aerodynamic loads are determined using Peters' unsteady aerodynamic model. This model is based on the assumption of incompressible potential flow, and consequently, the effects of viscosity have been neglected. 35 The aerodynamic moment and force can be written as 36   f a1 f a2 f a3 where a is referred to the aerodynamic reference axis, b is the semi-chord, and ρ is the air density. Also, λ is a column matrix which includes inflow states, and ½A inflow , fb inflow g, fc inflow g are constant matrices derived in Peters and Karunamoorthy 34 work. Also, C a is the rotation matrix, and y mc is a row matrix (position vector) from the beam reference axis to the mid-chord and can be written as Aerodynamic loads on the aerodynamic coordinate are transferred to the beam reference frame using the following equations The boundary conditions are written as equation (9) V ¼ ½0 To take into account any kind of nodal discontinuities such as external mass or force, the following equation is considered.
where b Cl r is the slope discontinuity. 36

Solution methodology
Since the governing equations are composed of algebraic and partial differential equations, a finite difference discretizing scheme is utilized for solving these equations (as seen in the work of Hodges 27 as well as Chang 36 ). The resulting discretized equations can be rewritten in a compact form as where fXg is a vector and it contains the structural and aerodynamic variables.
To determine the stability of the wing, first, the nonlinear steady-state of the system is obtained by dropping all time-dependent variable from equation (11) and solving the resulting nonlinear equation by using the Newton-Raphson method. Then, the dynamic equations are linearized about this steady-state condition, and the stability of the wing is sought by investigating the eigenvalues of this linearized system.
It should be noted that after determining the length change in each step by equation (1), the amount of mass change of the fixed part is proportional to the size of the change in its length. Also, the number of finite elements is fixed, and hence, the size of the elements changes at each step.

Results and discussion
To check the validity of the developed aeroelastic model, first, the flutter speed and frequency of two wings with low and high aspect ratios are obtained and compared with those reported in the literature. The wing parameters used are listed in Table 1. According to Table 2, it is clear that the present results are in very good agreement with those presented by Patil 37 with a maximum difference of 0.3%. It is noted that Patil considered this problem by solving the mixed beam formulation using a finite element approach. It is noted that in the present study, six aerodynamic inflow states have been used. Furthermore, Figure 3 shows the convergence of the flutter speed for different numbers of elements. It is clear that by using 25 elements, the flutter speed can be predicted accurately, and hence, from here on, 25 elements are used for all case studies.
Next, the effect of wing angle of attack on the flutter speed of the Goland wing is investigated and the results are depicted in Figure 4 in accordance with those reported by Patil and Hodges. 38 This clearly shows that using a linear beam model results in inaccurate prediction of flutter speed when the wing has an initial angle of attack.
The effect of engine thrust (p) on the flutter speed of a high aspect ratio wing with the parameters presented in Table 1 is considered next for various nondimensional thrust ratios Ψ ¼ It should be noted that the    thrust is applied at 15 m from the root. Figure 5 compares the present results with those calculated by Hodges et al. 39 which shows a good agreement.
It should be noted that v f ¼ u f bω α denotes reduced flutter speed, where ω α is the first uncoupled torsional frequency.
By considering all the above studies, it can be concluded that the present model is able to capture the nonlinear aeroelastic stability of wings with acceptable accuracy.
In what follows, it is assumed that the length of the Goland wing can be extended up to 50% due to the moving part in the telescopic wing. This implies that although the Goland wing is considered as a low aspect ratio wing, after extension, it becomes a higher aspect ratio wing, and hence, the effect of structural nonlinearities might become more effective. In this regard, the nonlinear trend of aeroelastic stability of a two-stage telescopic wing is investigated quasi-statically. To this aim, the following nondimensional parameters are used In the first step, the effect of wing length and overlapping mass, due to the motion of the moving part, on the reduced flutter speed is determined and shown in Figure 6. This figure shows that by considering overlapping mass, a 25% increase in span length leads to a difference of 1% to 8% in reduced flutter speed. On the other hand, as the length increases, the reduced flutter speed is reduced by up to 22%. This clearly shows that the morphing length and overlapping mass significantly affect the flutter speed of the wing.
It should be noted that in the following sections, an angle of attack of α ¼ 7 deg has been considered, and the effects of various system parameters on the wing stability are investigated. It is noted that in this study, the effects of stall have been ignored.
The effect of thrust force on the dimensionless critical length, (L f ), the length at which the wing gets unstable, at location of χ ¼ 0:5 without any offset of reference axis is determined and shown in Figure 7. It is noted that unless otherwise stated, from here on, all results are presented at a flight speed of 140 m/s, as this is slightly lower than the original wing flutter speed (157.4 m/s). Also, the effects of engine mass are ignored in the calculations.
This figure shows when the thrust force increases (decreases), the morphing wing critical length decreases (increases). Therefore, increasing the thrust force in a morphing wing can accelerate the onset of flutter. Also, the wing critical length decreases (increases) by decreasing (increasing) angle of attack, and it is more sensitive to the lower angle of attack. Moreover, as shown in Figure 7, frequency is directly related to the thrust force.
The effect of engine location on the dimensionless critical length of the wing is determined and shown in Figure 8. By moving the thrust force from the root to the tip of the fixed part and also increasing the morphing length, the critical length decreases, while the critical frequency increases. Therefore, moving the engine location from the root toward the tip of the wing can accelerate the onset of flutter and can change the critical length (critical frequency) up to about 12% (14%). This shows the importance and the effect of engine location on the performance of the morphing wing.
As the chord of the moving part is usually smaller than the fixed part chord, therefore it is necessary to investigate how this might affect the stability of the wing. Figure 9 shows the effect of moving part chord on the critical length for the various angles of attack, assuming the bending and     torsional stiffness values remain constant. As the chord of the moving part decreases, the wing critical length increases. A 40% decrease in chord length leads to variation of critical length of up to about 7%. Also, as the moving part chord decreases, the critical length difference increases for higher angle of attacks. Moreover, the wing critical frequency and moving part chord are directly related to each other, and hence, an increase in the moving part chord results in an increase in the critical frequency. A 40% increase in moving part chord length leads to variation of critical frequency of up to about 6%.
One of the goals of changing the morphing span (change in aspect ratio) is the possibility of changing the missions 40 and the lift to drag ratio. 41 Thus, a change in wing aspect ratio would result in a change in both range and endurance. 42 To change the range and endurance, the flight altitude is important. Figure 10 shows that the wing critical length increases when the flight altitude increases. Also, by increasing the length by 50%, the altitude can increase up to more than 7 km, which depends on the ability to increase the length of the morphing wing. Moreover, the critical frequency and flight altitude are inversely related to each other, and hence, an increase in flight altitude results in reduction of the critical frequency.

Conclusions
In this study, the nonlinear aeroelastic stability of a two-stage spanwise morphing wing has been studied. For structural and aerodynamic modeling, the geometrically exact fully intrinsic beam equations and Peters' unsteady aerodynamic model have been used, respectively. The resulting nonlinear equations have been discretized using the finite difference method. The stability of the system with different parameters has been evaluated by investigating the eigenvalues of the linearized system. The obtained results were compared with those available in the literature, and a good agreement was observed. Moreover, the effects of different parameters such as thrust force, engine location, chord size, and flight altitude in the presence of angle of attack and overlapping mass effect have been investigated. The results of this study are summarized as follows: 1. The morphing length and overlapping mass have significant effects on the aeroelastic stability of the telescopic wing. 2. The wing moving part chord, thrust force, and engine location inversely (directly) change the critical length (critical frequency) of the wing. Despite, wing angle of attack directly (inversely) affects the critical length (critical frequency). Therefore, the amount and point of effect of each parameter are effective on the increase in length. 3. By directly affecting the altitude change on the length change, the ability to increase the altitude can be strengthened. 4. By properly selecting each of the morphing parameters (due to the effects of length changing), the flight performance of the morphing wing can be improved.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.