Multi-body dynamics modeling and transient characteristics analysis for a folding wing

Abstract

The investigation of modeling the time evolution of a folding wing during the morphing process and the dynamic characteristics analysis is carried out. The governing equations with uniform form are developed from the integration of floating frame method in multi-body dynamics and component mode synthesis in structural dynamics. The time-dependent aerodynamic force is taken into the governing equations. The equation achieves the time-dependent coupling between structure and aerodynamics and avoids the data transmission and low efficiency, which holds true for the multi-segmented folding wing. The relative parameters in constraint equations are easily modified to be applied to both slow and fast-varying processes for a folding wing. Also, the influence of the velocity and attack angle on transient responses can be investigated. Transient response analysis shows that slower morphing means more stable transient responses. The flexibility of the folding wing has the significant influence on transient responses. To some extent, the aerodynamic force can be beneficial to the morphing process.

Keywords: Folding wing, multi-body dynamics, structural dynamics, aerodynamics, time-dependent coupling

1. Introduction

Some parameters such as position, mass center and inertia for different segment will be time-varying while the folding morphing concept is to change an aircraft’s wing configuration, resulting in significant changes of the dynamics for the design of the folding wing. Meanwhile, the investigation on dynamic characteristic plays an important role in the design and dynamic response improvement of a folding wing as well as the design of a control system.

Obradovic et al. [1-2] presented a simulation methodology from the integration of rigid-body dynamics and unsteady vortex-lattice method due to complex configuration for a gull-wing configuration, resulting in a six degree of freedom nonlinear equation of motion in which the changes in both the aerodynamic and inertial properties were taken into account. They analyzed the dynamics through nonlinear simulations. Selitrennik et al. [3] proposed an aeroelastic simulation of rapid morphing air vehicles, in which the time-varying coordinate transformation matrix was used to perform the fictitious-mass modal synthesis and the aerodynamic force was based on CFD. Ameri et al. [4] developed a multi-body modeling tool, in which aerodynamic loads based on quasi-steady linear vortex-lattice method was implemented. They investigated the dynamic and shape changing response of morphing aircraft and analyzed the effect of the variation in the shape on the inertia properties, the center of gravity location and modal parameters. Seigler et al. [5] carried out the dynamic analysis for a large-scale morphing aircraft based on the integration of the multi-body dynamic, aerodynamic theory and flight control theory. They focused on the effect of the time-varying aerodynamic force, inertial force and moment on the transient response. Zhao et al. [6] established a set of differential-algebraic equations that governed the time evolution of the folding wing during the morphing process by the integration of the floating frame method and modal component synthesis. They obtained the transient responses of the folding wing without aerodynamic force during a morphing process. Wickenheiser et al. [7] demonstrated the transient and steady-state behavior of a perching aircraft undergoing changes in flight condition and vehicle reconfiguration through analytical and empirical aerodynamic analysis, trim and stability analyses and flight simulation. Reich et al. [8] and Scarlett et al. [9] proposed an integated simulation approach of a unmanned combat aerial vehicle through the integation of finite element modelling, aerodynamic loading and continuous-time multi-body simulation. They mainly used the commerical software platform that was not convenient to integrate the time-varying aerodynamic forces into the structural dynamic system for a folding wing. Yue et al. [10] deduced the six degree of freedom nonlinear equation for a tailless folding wing morphing aircraft in wing folding process. They investigated the key factors which affected the dynamic charateristics of morphing process by quantity.

What have been obvious in the current dynamic research for a morphing wing is that there are two kinds of methods. One is that the six degree of freedom equations are obtained by some modifications in traditional dynamic modeling such as the inclusion of the variations in inertial properties, the center of gravity location and modal parameters, then achieving the dynamic response. However, this method is not accurate and comprehensive for the morphing process. Moreover, time-varying terms will result in the complex governing equations, implying the difficulties in finding the solution. The other is the establishment of the multi-body dynamic model from the integration the finite element method, multi-body dynamic theory and coupling the varied aerodynamic force. Nevertheless the time-varying aerodynamic is achieved by the interpolation of steady-state aerodynamic force at different configurations, implying that the effect of aerodynamic force on the transient response cannot be taken into consideration. Meanwhile, the structure dynamic, the aerodynamic force and the coupling between the structure and aerodynamic all is time-varying, indicating that it is not beneficial to the data transferring between commercial softwares and adds to the difficulty in modeling and analysis.

This paper describes the development of the transient dynamics modeling for a folding wing and the analysis of displacement and mechanic response is carried out. Firstly, the dynamic formulation for each component is established by integration of the floating frame approach in multi-body dynamics and component mode synthesis in structural dynamics. Then, the configuration for folding wing is expressed as the generalized coordinates in floating frame approach. According to current configuration, the unsteady aerodynamic force is achieved on the base of rational function approximation to Doublet Lattice Method, in which the angle of attack is taken into consideration. The equivalent structural force on the node is obtained by the virtual principle. Finally, the equations that govern the time evolution of the folding wing during the morphing process are established in MATLAB. The constraint equations are used to present the motion relation between each component. Investigations of time histories are presented for changing morphing velocity, free flow velocity and angle of attack.

2. Multi-body dynamic equations for folding wing

As shown in Figure 1, the folding wing consists of three separate components: the fuselage, the inboard wing and the outboard wing, denoted A, B, C, respectively. $oxyz$ is the inertial frame, which is used to describe the motion of the system. $o_{A}x_{A}y_{A}z_{A}$， $o_{B}x_{B}y_{B}z_{B}$， $o_{C}x_{C}y_{C}z_{C}$ are the floating frame to describe the elastic deform. ${\theta }^{A}$, ${\theta }^{B}$, ${\theta }^{C}$ are the folding angles of component A, B, C, respectively, which are positive in counter clockwise direction. All components are modeled by plates with the Young’s modulus of $7.1\times {10}^{10}$ Pa, Poisson’s coefficient of 0.33 and a density of $2.7\times {10}^{3}$kg/m$^{3}$, which are discretized using CQUAD4 elements with $2$mm thickness in MSC.NASTRAN. Each node of this element has six degrees of freedom: three translations, namely, $u_{x},u_{y},u_{z}$ and three rotations, i.e., ${\theta }_{x},{\theta }_{y},{\theta }_{z}$. For each component, all nodes at the interface are used to coordinate the boundary displacement. In this paper, the interfaces are assumed to be rigid. The finite element model for a folding wing is shown in Figure 2 .

(a)	(b)
Figure 1. Schematic of folding wing and coordinate systems. (a) Top view. (b) Rear view

Figure 2. Finite element model of folding wing

By integration of the component modal synthesis and the floating frame approach, the equations of motion of the flexible component can be written as：

$\left\{{\begin{array}{l}{{\bf {M}}^{A}{\ddot {\bf {q}}}^{A}+{\bf {D}}^{A}{\dot {\bf {q}}}^{A}+{\bf {K}}^{A}{\bf {q}}^{A}={\bf {Q}}_{c}^{A}+{\bf {Q}}_{g}^{A}+{\bf {Q}}_{v}^{A}+{\bf {Q}}_{a}^{A}}\\{{\bf {M}}^{B}{\ddot {\bf {q}}}^{B}+{\bf {D}}^{B}{\dot {\bf {q}}}^{B}+{\bf {K}}^{B}{\bf {q}}^{B}={\bf {Q}}_{c}^{B}+{\bf {Q}}_{g}^{B}+{\bf {Q}}_{v}^{B}+{\bf {Q}}_{a}^{B}}\\{{\bf {M}}^{C}{\ddot {\bf {q}}}^{C}+{\bf {D}}^{C}{\dot {\bf {q}}}^{C}+{\bf {K}}^{C}{\bf {q}}^{C}={\bf {Q}}_{c}^{C}+{\bf {Q}}_{g}^{C}+{\bf {Q}}_{v}^{C}+{\bf {Q}}_{a}^{C}}\end{array}}\right.$

(1)

where Failed to parse (syntax error): \textbf{M}^i，\, \textbf{D}^i，\, \textbf{K}^i

（$i=A,B,C$） are the generalized mass matrix, generalized damping matrix and generalized stiffness matrix, respectively,  ${\textbf {Q}}_{c}^{i}$ is the vector of constraint force,  ${\textbf {Q}}_{g}^{i}$ refers to the vector related to the gravity, ${\textbf {Q}}_{v}^{i}$ represents the quadratic velocity vector that arises from the differentiating the kinetic energy with respect to time and with respect to the generalized coordinates, ${\textbf {Q}}_{a}^{i}$ denotes the generalized aerodynamic force and  ${\textbf {q}}^{i}$ is the generalized displacement of each component and it is written in detail.

${\textbf {q}}^{i}={\left[{\begin{array}{ccc}{\textbf {R}}^{iT}&{\theta }^{iT}&{\textbf {q}}_{f}^{iT}\end{array}}\right]}^{T}$

(2)

where ${\textbf {R}}^{i}$ is the position vector of the origin of each floating frame, ${\textbf {q}}_{f}^{i}$ is the modal displacement, ${\textbf {D}}_{ff}$ and ${\textbf {K}}_{ff}$, which are the sub matrices in generalized damping matrix and generalized stiffness matrix, will be different when taken the aerodynamic force into consideration. The displacement of each component consists of large rigid body motions and small elastic deformations. Therefore, the mass matrices are function of the generalized displacement.

3. Unsteady aerodynamic forces

An et al. [11] compared the quasi-steady and unsteady aerodynamic charateristics of variable swept wing aircraft in wing sweep angular velocity by 8 deg/s, 10 deg/s, 15 deg/s, 20 deg/s, and the aircraft dynamic response results were quite similar. Therefore, the dynamic response affected by unsteady aerodynamic force during morphing process can be ignored within a certain range of morphing velocity. Although Doublet Lattice Method is an aerodynamic finite element method for modeling oscillating lifting surfaces in subsonic flows, the unsteady aerodynamic force resulting from the small oscillation of lifting surface is mainly taken into consideration while the morphing velocity is rather small and the change of aerodynamic force is not significant. As the preliminary research of the aerodynamic force, the Doublet Lattice Method is adopted here. The equivalent structural force on each node can be achieved on the base of virtual principle. While taking account of the effect of angle of attack, the equivalent force can be written:

${\textbf {Q}}_{a}={\frac {1}{2}}\rho V^{2}{\textbf {Q}}{\textbf {q}}_{f}+$${\frac {1}{2}}\rho V^{2}{\textbf {G}}\alpha$

(3)

where ${\textbf {Q}}$ is the generalized aerodynamic influence coefficient matrix, ${\textbf {G}}=q_{d}{\textbf {G}}_{n}^{T}\left({\textbf {G}}_{s1}^{T}-\right.$$\left.{\textbf {G}}_{s2}^{T}{\textbf {C}}_{\Delta }\right){\textbf {S}}{\textbf {D}}_{aic}^{-1}$, ${\textbf {G}}_{n}$ is a matrix used to extract the normal displace from structural node and ${\textbf {G}}_{s1}$ and ${\textbf {G}}_{s2}$ are the spline matrices generated by the Infinite Plate Spline (IPS) method, $\alpha$ is the angle of attack. Hence ${\textbf {Q}}_{a}$ and ${\textbf {q}}$ are relative with $\alpha$ and are time-dependent vectors. ${\textbf {S}}$ is the area matrix of aerodynamic panels.

4. Equations of motion of the folding wing during morphing process

Here the equivalent force on structural node cannot be directly used to the response analysis in time domain, because ${\textbf {Q}}$ in Equation (3) is in frequency domain. In order to convert ${\textbf {Q}}$ to time domain, there are many methods of rational function approximation [12-14], but the Minimum-State Approximation is used here for simplicity and effectiveness. The approximation form is as follows:

${\textbf {Q}}({\overline {s}})\approx {\textbf {A}}_{0}+$${\textbf {A}}_{1}{\overline {s}}+{\textbf {A}}_{2}{\overline {s}}^{2}+$${\textbf {D}}_{s}{\left({\textbf {I}}{\overline {s}}-{\textbf {R}}_{s}\right)}^{-1}{\textbf {E}}_{s}{\overline {s}}$

(4)

Here, ${\textbf {A}}_{0}$, ${\textbf {A}}_{1}$, ${\textbf {A}}_{2}$, ${\textbf {D}}_{s}$ and ${\textbf {E}}_{s}$ are unknown matrices in the minimum state approximation, which can be calculated from the least square fit. ${\textbf {R}}_{s}$ is the aerodynamic lag matrix.

Introducing aerodynamic state variable ${\textbf {q}}_{a}$ yields：

$\left\{{\begin{array}{l}{\textbf {M}}_{ff}^{i}{\ddot {\textbf {q}}}_{f}^{i}{}_{f}+{\textbf {D}}_{ff}^{i}{\dot {\textbf {q}}}_{f}^{i}{}_{f}+{\textbf {K}}_{ff}^{i}{\textbf {q}}_{f}^{i}={\frac {1}{2}}\rho V^{2}{\textbf {D}}_{s}{\textbf {q}}_{a}+{\frac {1}{2}}\rho V^{2}{\textbf {G}}\alpha \\{\dot {\textbf {q}}}_{a}={\textbf {R}}_{s}{\frac {V}{b}}{\textbf {q}}_{a}+{\textbf {E}}_{s}{\textbf {q}}_{f}\end{array}}\right.$

(5)

where $b$ is the reference length, ${\textbf {M}}_{ff}^{i}={\textbf {M}}_{f}^{i}-{\frac {1}{2}}\rho V^{2}{\left(b/V\right)}^{2}{\textbf {A}}_{2}$, ${\textbf {D}}_{ff}^{i}={\textbf {D}}_{f}^{i}-{\frac {1}{2}}\rho V^{2}\left(b/V\right){\textbf {A}}_{1}$, ${\textbf {K}}_{ff}^{i}={\textbf {K}}_{f}^{i}-{\frac {1}{2}}\rho V^{2}{\textbf {A}}_{0}$, ${\textbf {M}}_{f}^{i}$, ${\textbf {D}}_{f}^{i}$, ${\textbf {K}}_{f}^{i}$ are modal mass matrix, modal damping matrix and modal stiffness matrix, which are still the sub matrices in generalized damping matrix and generalized stiffness matrix. There exist two kinds of constraints for morphing motion, namely the kinematic constraints and the driving constrains. They can be expressed as

${\textbf {C}}({\textbf {q}},t)=\left[{\begin{array}{c}{\textbf {R}}^{A}\\{\theta }^{A}\\{\textbf {R}}^{B}-{\textbf {f}}_{c}^{B}\\{\theta }^{B}-f_{d}^{B}\\{\textbf {R}}^{C}-{\textbf {f}}_{c}^{C}\\{\theta }^{C}\end{array}}\right]={\textbf {0}}$

(6)

where

${\textbf {R}}^{A}=\left[{\begin{array}{c}R_{x}^{A}\\R_{y}^{A}\\R_{z}^{A}\end{array}}\right]$ ， ${\textbf {R}}^{B}=\left[{\begin{array}{c}R_{x}^{B}\\R_{y}^{B}\\R_{z}^{B}\end{array}}\right]$ ， ${\textbf {R}}^{C}=\left[{\begin{array}{c}R_{x}^{C}\\R_{y}^{C}\\R_{z}^{C}\end{array}}\right]$ ， ${\textbf {f}}_{c}^{B}=\left[{\begin{array}{c}l_{A}\tan(\Lambda )\\l_{A}\\0\end{array}}\right]$ ， ${\textbf {f}}_{c}^{C}=\left[{\begin{array}{c}R_{x}^{B}+l_{B}\tan(\Lambda )\\R_{y}^{B}+l_{B}\cos({\theta }^{B})\\R_{z}^{B}+l_{B}\sin({\theta }^{B})\end{array}}\right]$

$f_{d}^{B}$ refers to the desired folding angle. Components A and C are kept horizontal. Considering the constraint equation, ${\textbf {Q}}_{c}$ can be replaced with $-{\textbf {C}}_{q}^{T}{\boldsymbol {\lambda }}$, introducing the vector of Lagrange multipliers $\lambda$. ${\textbf {C}}_{q}$ is the constraint Jacobian equation. The resulting equation for morphing process can be as follows:

$\left\{{\begin{array}{l}{\textbf {M}}{\ddot {\textbf {q}}}+{\textbf {D}}{\dot {\textbf {q}}}+{\textbf {Kq}}+{\textbf {C}}_{q}^{T}{\boldsymbol {\lambda }}={\textbf {Q}}_{g}+{\textbf {Q}}_{v}+{\textbf {Q}}_{a}\\{\dot {\textbf {q}}}_{a}={\textbf {R}}_{s}{\frac {V}{b}}{\textbf {q}}_{a}+{\textbf {E}}_{s}{\textbf {q}}_{f}\\{\textbf {C}}({\textbf {q}},t)={\textbf {0}}\end{array}}\right.$

(7)

where

${\textbf {M}}=blockdiag({\textbf {M}}^{A},{\textbf {M}}^{B},{\textbf {M}}^{C})\quad ,\quad$ ${\textbf {D}}=blockdiag({\textbf {D}}^{A},{\textbf {D}}^{B},{\textbf {D}}^{C})\quad ,\quad$ ${\textbf {K}}=blockdiag({\textbf {K}}^{A},{\textbf {K}}^{B},{\textbf {K}}^{C})$

Evidently, the equation above is a set of differential-algebraic equations. It can be solved by the method proposed by Baumgarte [15]. Differentiating Equation (6) with respect to time yields:

${\dot {\textbf {C}}}({\textbf {q}},t)={\textbf {C}}_{q}({\textbf {q}},t){\dot {\textbf {q}}}+$${\textbf {C}}_{t}({\textbf {q}},t)={\textbf {0}}$

(8)

Differentiating Equation (7) with respect to time again yields:

${\ddot {\textbf {C}}}({\textbf {q}},t)={\textbf {C}}_{q}({\textbf {q}},t){\ddot {\textbf {q}}}+$${\left({\textbf {C}}_{q}({\textbf {q}},t){\dot {\textbf {q}}}\right)}_{q}{\dot {\textbf {q}}}+$$2{\textbf {C}}_{q}({\textbf {q}},t){\dot {\textbf {q}}}+$${\textbf {C}}_{tt}={\textbf {0}}$

(9)

Combining Equations (8) and (9),One can rewrite Equation(7):

$\left[{\begin{array}{ccc}{\textbf {M}}&{\textbf {0}}&{\textbf {C}}_{q}^{T}\\{\textbf {0}}&{\textbf {I}}&{\textbf {0}}\\{\textbf {C}}_{q}&{\textbf {0}}&{\textbf {0}}\end{array}}\right]\left[{\begin{array}{c}{\ddot {\textbf {q}}}\\{\textbf {q}}_{a}\\{\dot {\boldsymbol {\lambda }}}\end{array}}\right]=\left[{\begin{array}{c}-{\dot {\textbf {D{q}}}}-{\textbf {Kq}}+{\textbf {Q}}_{g}+{\textbf {Q}}_{v}+{\textbf {Q}}_{a}\\{\textbf {R}}_{s}{\frac {V}{b}}{\textbf {q}}_{a}+{\textbf {E}}_{s}{\textbf {q}}_{f}\\{\textbf {Q}}_{c}-2c_{1}{\dot {\textbf {C}}}({\textbf {q}},t)-c_{2}^{2}{\textbf {C}}({\textbf {q}},t)\end{array}}\right]$

(10)

where ${\textbf {Q}}_{c}=-{\textbf {C}}_{q}({\textbf {q}},t){\ddot {\textbf {q}}}+$${\left({\textbf {C}}_{q}({\textbf {q}},t){\dot {\textbf {q}}}\right)}_{q}{\dot {\textbf {q}}}+$$2{\textbf {C}}_{q}({\textbf {q}},t){\dot {\textbf {q}}}+$${\textbf {C}}_{tt}$. $c_{1}$ and $c_{2}$ are positive and play the rule of control terms which achieves the stability.

5. Transient response analyses for a folding wing

To predict the transient response for a folding wing, the folding angle ${\theta }^{B}$ can meet the following rule.

${\theta }^{B}=\left\{{\begin{array}{c}{120}^{\mbox{°}}({\frac {t}{T_{s}}}-{\frac {1}{2\pi }}sin{\frac {2\pi t}{T_{s}}}),(0\leq t\leq T_{s})\\{120}^{\mbox{°}},(t>T_{s})\end{array}}\right.$

(11)

where ${120}^{\mbox{°}}$ stands for the steady state angle, $T_{s}$ is the time that the folding angle ${\theta }^{B}$ goes from $0^{\mbox{°}}$ to ${120}^{\mbox{°}}$. The whole program is established in MATLAB, which can be solved by the computational algorithm in Figure 3 . The first six modes of each component are used and the modal damping ratio is 0.02. Therefore, there are 46 equations in Equation (9). $c_{1}$ and $c_{2}$ can be set to 5 and 10, respectively. To determine the effect of aerodynamic force, the transient response of the rigid and flexible wing with gravity can be performed first, then the one with gravity and aerodynamic force.

Figure 3. Computational algorithm for the dynamic analysis

5.1. Transient response of rigid folding wing with gravity

To determine the effect of aerodynamic force on the transient response, the response of rigid folding wing with gravity can be carried out firstly. Figure 4 shows the vertical displacement of the wingtip. The final displacement is 118.2mm for different $T_{s}$. When $T_{s}$ is set to 2s, 4s, 6s, ${\theta }^{B}$ equals ${90}^{\mbox{°}}$, implying that component B is perpendicular to components A and C and the displacement is the maximum. The greater $T_{s}$, the slower the morphing process.

Figure 4. Displacement response of the wingtip of rigid folding wing with gravity

During the morphing process, the inertia tensors of each component vary relative to the inertial coordinate. Superimposing the tensor of each part, which is obtained in the inertial coordinate, the inertia tensor of the folding wing can be obtained. As shown in Figure 5, it can be seen that for different $T_{s}$ , the inertia tensor remains constant at the end of the folding motion as the time changes. And the rotational inertia Izz, Iyy, Ixx in the three directions sequentially decrease and remain positive, which is consistent with the definition of the moment of inertia, while the inertial product Ixy, Ixz, Iyz are negative. The more slowly the folding motion is, the slower each component of the inertia tensors varies.

(a)	(b)
Figure 5. Time history of inertial tensor. (a) $Ts=2s$. (b) $Ts=6s$

The mechanical responses of the rigid folding wing with gravity are shown in Figure 6.

(a)	(b)
(c)	(d)
Figure 6. Mechanical response of rigid folding wing with gravity. (a) Resulting reaction force acting on rigid component B. (b) Resulting reaction moment acting on rigid component B. (c) Reaction force acting on rigid component C. (d) Reaction moment acting on rigid component C

In Figure 6, the resulting reaction force and resulting reaction moment acting on the rigid component A are not shown, since the resulting reaction force is the gravity of the component A, i.e. 6.542 $N$, the resulting reaction moment resulting from gravity, namely 427.6 $Nmm$. The motion of the rigid component A in the floating frame is the same with that of the inertial coordinate. $T_{s}$ has no effect on the mechanical response of component A. So the resulting reaction force and resulting reaction moment acting on the rigid component A are not shown in Figure6. It can be seen, however, that the stable steady value of the resulting reaction force on the rigid component B is equal to the gravity of the component B at the beginning and the end of the folding motion, i.e., 4.89 $N$. And the resulting reaction moment acting on the rigid substructure B is about 341.36 $Nmm$ at the beginning and -157.18 $Nmm$ in the end. The negative value represents the moment in the opposite direction. The initial and steady values are equal to the moments resulting from the gravity of the component B. At 1.27 $s$, 2.5 $s$ , 3.81 $s$, the folding angle ${\theta }^{B}$ equals ${90}^{\mbox{°}}$, implying that the center of gravity passes through the floating frame, and the gravity produces a moment of zero. In Figure6, both the initial reaction force and the steady one acting on the rigid component C are the same, i.e. 11.564 $N$. The initial and steady values of the reaction moment are the same, namely 2885 $Nmm$. The reason is that the component C is always kept in horizontality during the morphing process. Hence the center of gravity of the component C does not change much, relative to its floating frame. Either the resulting reaction force and moment acting on the rigid substructure B or the reaction force and the moment acting on the rigid substructure C vary with the increase of time. Moreover, the smaller $T_{s}$, the more significant the variation. Due to the smaller $T_{s}$, the deformation velocity and the quadratic velocity vector that arises from the differentiating the kinetic energy with respect to time and with respect to the generalized coordinates is greater, so the resulting force is greater. It is obvious that the slower morphing help to reduce the force acting on the rigid components B and C. Nevertheless the decrease of moment is not obvious, because ${\textbf {Q}}_{v}$ points to the axis, producing no additional moment.

5.2 Transient response of rigid folding wing with gravity and aerodynamic

Then, the transient responses of the rigid folding wing with gravity and aerodynamic force is compared with ones in Section 5.1 to reveal the effect of aerodynamic force on the transient response. When the morphing process is relatively slow, the unsteady effect of the aerodynamic force can be ignored. Hence, here $T_{s}$ is set to $6s$ , indicating that the process is relatively slow. Since the folding wing is rigid, the displacement response of the wingtip is consistent with one of the rigid folding wing with the gravity, implying that the vertical displacement response of the wingtip with gravity and aerodynamic force is shown in Figure 4, the corresponding displacement response corresponds to $T_{s}=6s$. Therefore, only the mechanical response is analyzed at different velocities and angles of attack, respectively, as shown in Figures 7 and 8. The free flow velocity is set to 20 $m/s$, 25 $m/s$ and 30 $m/s$, and the angle of attack equals to $1^{\circ }$, $2^{\circ }$ and $3^{\circ }$. Compared with the mechanical responses of the rigid folding wing with gravity, namely, the corresponding curves in Figure 6, it can be found that the aerodynamic force changes the mechanical response.

(a)	(b)
(c)	(d)
(e)	(f)
Figure 7. Mechanical response of rigid folding wing at different velocities. (a) Resulting reaction force acting on rigid component A. (b) Resulting reaction moment acting on rigid component A. (c)Resulting reaction force acting on rigid component B. (d) Resulting reaction moment acting on rigid component B; (e) Reaction force acting on rigid component C. (f) Reaction moment acting on rigid component C

As shown in Figure 7, within the constant angle of attack and a certain free flow velocity range, as the velocity increases, the resulting reaction force and moment on the components A and C as well as the reaction forces and moments acting on the component C are reduced with gravity and aerodynamic force when the folding motion stabilizes. For the component B, the resulting reaction force and moment decrease as the velocity increases with the folding angle less than ${90}^{\circ }$. On the contrary, those increases with the folding angle more than ${90}^{\circ }$. The reason for this phenomenon is as follows. The components A and C always keep horizontal during the morphing process, in which the direction of the aerodynamic force is opposite to the direction of gravity. Hence the direction of gravity is always opposite to the one of aerodynamic force. And as the velocity increases, the aerodynamic force is also increasing, so the resulting force and moment acting on the components A and C are reduced. The direction of the aerodynamic force of the component B is normal along the wing surface in the folding process. There is a component of the aerodynamic force in the positive direction of the axis Z of the inertial coordinate with the folding angle less than ${90}^{\circ }$, implying that it can weaken the influence of gravity. However, the aerodynamic force has a component in the negative direction of the Z axis of the inertial coordinate with the folding angle more than ${90}^{\circ }$. Therefore, at the beginning, the resulting reaction force and moment of the component B are reduced with the increase of the velocity. At the end, the resultant force and moment increase. For the component A, the resulting reaction force and moment are not constant. This is due to the fact that the aerodynamic forces of the various components of the entire folding wing are interrelated, and the aerodynamic force on the component A will change as the configuration of the wing changes. Therefore, the resulting reaction force and moment on component A are slowly reduced and then remain stable. For the component C, it can be seen that the reaction force and moment increase first and then decrease from the partial magnification of (e) and (f). The reason is partly similar to the one for the mechanical response of a rigid folding wing with gravity. Another part of the reason is due to the folding wing in different configurations, the aerodynamic changes on component C.

The mechanical response of the rigid folding wing with different angle of attack is shown in Figure 8 with the constant free flow velocity. The mechanical response is similar to ones in Figure 7. The angle of attack increases and the aerodynamic force effect increases. The resultant force and moment are reduced.

(a)	(b)
(c)	(d)
(e)	(f)
Figure 8. Mechanical response of rigid folding wing at different angle of attack. (a) Resulting reaction force acting on rigid component A. (b) Resulting reaction moment acting on rigid component A. (c) Resulting reaction force acting on rigid component B. (d) Resulting reaction moment acting on rigid component B. (e) Reaction force acting on rigid component C. (f) Reaction moment acting on rigid component C

From Figure 8, it can be seen that with the increase of the angle of attack, the resulting reaction force and moment acting on the rigid component A and the reaction force and moment acting on the rigid component C are reduced at the constant free flow velocity. For the component B, with the increase of the angle of attack, the resulting reaction force and moment are reduced with the folding angle less than ${90}^{\circ }$. Meanwhile, the resulting reaction force and moment is increased with the folding angle more than ${90}^{\circ }$. The cause of this phenomenon is similar to the reason for the mechanical response of a rigid folding wing at different velocity. It is worth noting that for a rigid component C, the reaction force and moment become negative with the angle of attack $3^{\circ }$. From the aerodynamic theory, within a certain angle of attack range, the lift increases with the increase of the angle of attack. The aerodynamic force acting on the component C is greater than its gravity, so the reaction force is downward, and the reaction moment is opposite to the X-axis direction of the inertial coordinate.

5.3 Transient response of flexible folding wing with gravity

In order to analyze the effect of aerodynamic force on the transient response of a flexible folding wing, the aerodynamic effect is not taken into account, implying that, the aerodynamic relation is neglected in Equation (10). And transient response of the flexible folding wing is analyzed with gravity. The vertical displacement response of the wingtip is shown in Figure 9. When $T_{s}$ is set to 2s and 6s, respectively, the component B is perpendicular to the components A and C in the vicinity of 1.27s, 3.81s and the vertical displacement of wingtip reaches the maximum. Obviously, the flexibility has a great effect on the displacement response of the wingtip. During the morphing process, the displacement of the wingtip is prone to convergence, and stability value is smaller than that of the rigid wing in the end. Due to the structural damping in Equation (10), the displacement response tends to converge. The larger $T_{s}$, the less the oscillation of normal displacement of the wingtip. The oscillation of the displacement of the flexible wing is obvious with gravity.

(a)	(b)
Figure 9. Displacement response of the wingtip of flexible wing with gravity. (a) $Ts=2s$. (b) $Ts=6s$

As shown in Figure 10 the flexibility of the wing has a great influence on the resulting reaction force and moment acting on the components A and B, and the reaction force and moment acting on the component C with gravity. Moreover, ${\theta }^{B}$ varies from $0^{\mbox{°}}$ to ${90}^{\mbox{°}}$, the mechanical response shows a trend that is first oscillation convergence and then divergence. After ${\theta }^{B}$ is greater than ${90}^{\mbox{°}}$ the mechanical response is convergent.

(a)	(b)
(c)	(d)
(e)	(f)
Figure 10. Mechanical response of flexible folding wing with gravity. (a) Resulting reaction force acting on flexible component A. (b) Resulting reaction moment acting on flexible component A. (c) Resulting reaction force acting on flexible component B. (d) Resulting reaction moment acting on flexible component B. (e) Reaction force acting on flexible component C. (f) Reaction moment acting on flexible component C

For the rigid component A, the resulting reaction force and moment are constant, while the resulting reaction force and moment has a substantial oscillation, treated as the flexible one. Moreover, the resulting reaction force and moment vibrate near 6.542 $N$, 427.75 $Nmm$, which are the force and moment arising from the gravity of the component A, indicating that the deformation of the flexible substructure A is small and its center of gravity position does not change significantly.

Compared with (c) and (d) in Figure 6, the resulting reaction force acting on the flexible component B has a large oscillation, which has changed the trend in the rigid state, but shown a convergence trend as time increases. For the resulting reaction moment, the trend is basically similar to the one in the rigid state, just with a little oscillation. Compared with (e) and (f) in Figure 6, the reaction force and moment acting on the flexible component C are tend to converge. Therefore, the flexibility of the wing folding has a significant effect on the response.

5.4 Transient response of flexible folding wing with gravity and aerodynamic force

Finally, the transient response of the flexible folding wing with gravity and aerodynamic force is obtained and analyzed at different free flow velocity and angle of attack, respectively, as shown in Figures 11 to 14. Compared with the results in Figures 9 and 10, the effect of aerodynamic force on the transient response of flexible folding wing is revealed. Here the selection of velocity and angle of attack is consistent with section 5.2. When the angle of attack is constant, the displacement response of the wingtip at different free flow velocity is shown in Figure11. Apparently the displacement response has a large oscillation due to the flexibility of the wing. But the oscillation at the end of morphing process decreases as the velocity increases in the partial enlarged view. This is mainly due to the existence of aerodynamic damping. With the increase in the free flow velocity, the aerodynamic damping is enhanced. Hence, the aerodynamic force is beneficial to the morphing process.

Figure 11. Displacement response of the wingtip at different velocities

(a)	(b)
(c)	(d)
(e)	(f)
Figure 12. Mechanical response of flexible folding wing at different velocities. (a) Resulting reaction force acting on flexible component A. (b) Resulting reaction moment acting on flexible component A. (c) Resulting reaction force acting on flexible component B. (d) Resulting reaction moment acting on flexible component B. (e) Reaction force acting on flexible component C. (f) Reaction moment acting on flexible component C

From Figure 12 , the resulting reaction force and moment acting on the flexible components A and B, the reaction force and the moment acting on the flexible component C are convergent with gravity and aerodynamic force. With the increase of velocity, the aerodynamic damping effect is enhanced and the oscillation is reduced. Compared with Figure 7, for the flexible component A, at the end of the morphing process, the resulting reaction force and moment are smaller under the same conditions. For example, with the free flow velocity 30 m/s, the resulting reaction force is 5.314 $N$, the resulting reaction moment is 340.4 $Nmm$ for the rigid condition, while the resulting reaction force and moment are 5.221 $N$ and 321.7 $Nmm$ for the flexible condition, respectively. For the flexible component B, only the resulting reaction force substantially centers on the results under the rigid condition throughout the morphing process and oscillates. The resulting reaction moment centers on the results under the rigid condition with the folding angle less than ${90}^{\mbox{°}}$ , while the resulting reaction moment is less than the one for the rigid wing with the folding angle more than ${90}^{\mbox{°}}$. For the flexible component C, the reaction force centers on the result of the rigid condition during the morphing process.

With gravity and aerodynamic force, the vertical displacement response of the wingtip at different angles of attack is shown in Figure 13. As the angle of attack increases, the oscillation of the wing tip is significantly reduced. When the angle of attack is set to $3^{\mbox{°}}$, the initial value of the displacement is positive, indicating that the aerodynamic force acting on the folding wing is greater than its gravity and the deformation is positive. When the angle of attack is $1^{\mbox{°}}$ and $2^{\mbox{°}}$, the aerodynamic force acting on the folding wing is smaller than its gravity, the deformation is negative. So the initial value of the wingtip is negative. With the increase of the angle of attack, the aerodynamic force acting on the wing is also increasing, the resultant force is reduced.

Figure 13. Displacement response of the wingtip at different angles of attack

Similar to the mechanical response of the flexible folding wing at different free flow velocity, with aerodynamic force, the resulting reaction force and moment acting on the flexible components A and B, the reaction force and moment acting on the flexible component C has the trend to convergence at different angles of attack in Figure 14. And with the increase of the angle of attack, the oscillation of the mechanical response is reduced due to the decrease in the elastic deformation. Compared with Figure 8, under the same conditions, at the end of morphing process, the resulting reaction force and moment on the flexible component A are smaller than those under the rigid condition. Taken the flexible component A as an example, with the angle of attack $3^{\mbox{°}}$ the resulting reaction force and moment are 3.676 $N$ and 206.4 $Nmm$ while the resulting reaction force and moment are 3.977 $N$ and 245.5 $Nmm$ under the rigid condition. For the flexible component C, the reaction force substantially centers on the one of the rigid wing and oscillates up and down. The reaction moment has a little different. With the angle of attack $1^{\mbox{°}}$, it centers on the results in the rigid condition. However, it is less than the results of the rigid wing with the angle of attack $3^{\mbox{°}}$. In the case of the rigid wing, the moment is -1117 $Nmm$, while the one is -843.2 $Nmm$ for the flexible wing. So to a certain extent, the aerodynamic force is beneficial to the morphing motion.

(a)	(b)
(c)	(d)
(e)	(f)
Figure 14. Mechanical response of flexible folding wing at different angles of attack. (a) Resulting reaction force acting on flexible component A. (b) Resulting reaction moment acting on flexible component A. (c) Resulting reaction force acting on flexible component B. (d) Resulting reaction moment acting on flexible component B. (e) Reaction force acting on flexible component C. (f) Reaction moment acting on flexible component C

6. Conclusions

In order to analyze the dynamic characteristics of the wing during the morphing process, the dynamic governing equation of the folding wing in the morphing process is established and the response is analyzed in detail. The equation is established on the base of the floating frame method in multi-body dynamics and the component modal synthesis method in the structural dynamics, in which time coupling of the structure and the aerodynamics is achieved, and the relevant parameters in the constraint equation can be changed to analyze the transient response of the rigid and flexible folding wing in the slow and fast-varying process. Introducing the aerodynamic force in time domain, the effect of the free flow velocity and angle of attack on the response is investigated. The conclusions can be drawn as follows:

1. The governing equation of the folding wing is established by the integration of the floating frame method and the component modal synthesis method to realize the time-varying coupling between the structure and the aerodynamics, which avoids the data transfer between the commercial software, the computational complexity and the computational efficiency. The proposed approach is applicable to multi-segmented wing, which is similar to folding wings.
2. In the case of rigid folding wing, the slower the morphing process, the more stable the transient response. For the flexible folding wing, the slower the morphing process, the smaller the oscillation of the transient response.
3. For the rigid wing, the resulting reaction force and moment acting on the components, as well as the reaction force and moment, are reduced with gravity and aerodynamic forces. So to a certain extent, the aerodynamic force can assist the morphing of the wing.
4. With gravity and aerodynamic force acting on the flexible wing, with the constant angle of attack, the free flow velocity is larger within a certain range, the smaller the oscillation of displacement response of the wingtip. With the constant free flow velocity, the greater the angle of attack, the smaller the deformation of the wing in the morphing process within a certain range. Moreover, the oscillations of the resultant force and the moment are reduced.

Author Contributions

Y.N. (Yingge Ni) and W.Z. (Wei Zhang) conceived the whole numerical simulations; Y.N. (Yingge Ni) wrote the paper. Y. L. (Yi Lv) reviewed and edited the paper.

Conflicts of Interest: None.

References

[1] Obradovic B., Subbarao K. Modeling of flight dynamics of morphing-wing aircraft. J. Aircr., 48(2):391-402, 2011.

[2] Obradovic B., Subbarao K. Modeling of dynamic loading of morphing-wing aircraft. J. Aircr., 48(2):424-435, 2011.

[3] Selitrennik E., Karpel M., Levy Y. Computational aeroelastic simulation of rapidly morphing air vehicles. J. Aircr., 49(6):1675-1686, 2012.

[4] Ameri N., Lowenberg M.H., Friswell M.I. Modelling the dynamic response of a morphing wing with active winglets. AIAA Atmospheric Flight Mechanics Conference and Exhibit, Hilton Head, South Carolina, AIAA 2007-6500, 20-23 August 2007, .

[5] Seigler T.M., Neal D.A., Bae Jae-Sung, Inman D.J. Modeling and flight control of large-scale morphing aircraft. J. Aircr., 44(4):1077-1087, 2007.

[6] Zhao Y.H., Hu H.Y. Prediction of transient responses of a folding wing during the morphing process. Aerosp Sci. Technol., 24:89-94, 2013.

[7] Wickenheiser A.M. Garcia E. Longitudinal dynamics of a perching aircraft. J. Aircr., 43(5):1386-1392, 2006.

[8] Reich G.W., Bowman J.C., Sanders B.. Development of an integrated aeroelastic multibody morphing simulation tool. 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Newport, Rhole Island, AIAA 2006-1892, 1-4 May 2006.

[9] Scarlett J.N., Canfield R.A. Multibody dynamic aeroelastic simulation of a folding wing aircraft. 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Newport, Rhole Island, AIAA 2006-2135, 1-4 May 2006.

[10] Yue T., Wang L.X. Multibody dynamic modeling and simulation of a tailless folding Wwng morphing aircraft. AIAA Atmospheric Flight Mechanics Conference, Chicago, Illinois, AIAA 2009-6155, 10-13 August 2009.

[11] An J.G., Yan M., Zhou W.B., Sun X.H., Yan Z. Aircraft dynamic response to variable wing sweep geometry. J. Aircr., 25(3):216-221, 1988.

[12] Jones R.T. The unsteady lift of a wing of finite aspect ratio. NACA Report-681, 1941.

[13] Roger K.L. Airplane math modeling methods for active control design. AGARD-CP 228, pp. 1-11, 1977.

[14] Karpel M. Extensions to the minimum-state aeroelastic modeling method. AIAA J., 29(11):2007-2009, 1991.

[15] Baumgarte J. Stabilization of constraints and integrals of motion in dynamical system. Comput. Method Appl. Mech. Eng., 1(1):1-16, 1972.