- Email: [email protected]
Aeroelastic study for folding wing during the morphing process
Journal of Sound and Vibration 365 (2016) 216–229
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
Aeroelastic study for folding wing during the morphing process Wei Hu, Zhichun Yang n, Yingsong Gu School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China
a r t i c l e in f o
abstract
Article history: Received 15 April 2015 Received in revised form 14 October 2015 Accepted 29 November 2015 Handling editor: W. Lacarbonara Available online 14 December 2015
This paper focuses on the aeroelastic characteristics of a folding wing during the morphing process. The folding wing structure is modeled by using the flexible multi-body dynamics approach, and an efficient method is proposed to calculate the aerodynamic force of the folding wing during the morphing process. The aerodynamic influence coefficient (AIC) matrices at different folding angles are obtained by the Doublet Lattice based aerodynamics theory, and then the orders of these AIC matrices are reduced by the spline interpolation technique. Through the minimum state approximation, the reduced AIC matrices are described as rational functions in the Laplace domain. Then the Kriging agent model technique is used to interpolate the coefficient matrices of the rational functions obtained from several different folding angles and to build the aerodynamics model in the time domain. At some different folding angles, the element values of the coefficient matrices before and after the interpolation are compared to verify the accuracy of the aerodynamics model, and then the aeroelastic responses of the folding wing during its morphing processes are simulated. The results demonstrate that the folding and unfolding processes have opposite influences on the dynamic aeroelastic stability of the folding wing, and the influences become much more significant with the increasing of the folding and unfolding rates. When the folding wing is morphing with a very slow rate, the dynamic aeroelastic stability will be similar to that obtained by the quasi-steady aeroelastic analysis. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Folding wing Aeroelasticity Flexible multi-body Time varying system Morphing wing
1. Introduction Since the beginning of the Morphing Aircraft Structures (MAS) Program, some novel morphing aircraft concepts have been proposed. Morphing aircraft holds several advantages over conventional aircraft designs including multi-mission capability, stemming from its ability of varying wing geometries during the flight. One of the famous morphing aircraft concepts is the folding wing aircraft, proposed by Lockheed Martin Corporation [1]. In recent years, many studies have been conducted on the aeroelastic characteristics of the folding wing aircraft. Several investigations [2–4] were performed to study the influences of structural parameters, such as the folding angle and the hinge spring stiffness, on the aeroelastic characteristics. Weisshaar and Lee [5] built a high fidelity aeroelastics model of the folding wing, and the results showed that the flutter dynamic pressure increases with the increasing of the inboard wing folding angle. For the extended wing configuration, the flutter dynamic pressure is much more sensitive to the changes in n
Corresponding author. Tel./fax: þ86 29 88460461. E-mail addresses: [email protected], [email protected] (Z. Yang).
http://dx.doi.org/10.1016/j.jsv.2015.11.043 0022-460X/& 2015 Elsevier Ltd. All rights reserved.
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
217
the inboard hinge stiffness, while for the folded wing configuration, the flutter dynamic pressure is much more sensitive to the changes in the outboard hinge stiffness. Zhao and Hu [6] proposed a parameterized aeroelastic modeling procedure, and the flutter characteristics of a folding wing with different configurations were investigated. Some studies on the nonlinear aeroelastic characteristics of the folding wing were also carried out. Lee and Chen [7] investigated the nonlinear aeroelastic characteristics of a folding wing with free-play hinge nonlinearity. In their study, even with a small free-play nonlinearity, the limit cycle oscillation (LCO) was observed when the folding angle is 0–301. And the LCO phenomenon would disappear when the folding angle is larger than 301. Tang, Dowell and Attar [8,9] investigated the influence of geometric nonlinearity on aeroelastic behaviors of a folding wing through experimental and computational methods, and the responses obtained by these two methods agreed well with each other. The result showed that the LCO appeared at air speed beyond the linear flutter speed of the folding wing, and the LCO phenomenon was changed at different folding angles. It is noted that all of these studies mentioned above were analyzed in the quasi-steady condition with fixed folding angles, so the results were only suitable for the folding wing during very slow morphing process. Few studies have explored whether these results are suitable for the folding wing with a rather rapid morphing rate in normal morphing process. Zhao and Hu [10] have developed the flexible multi-body dynamics formulation by combining the Craig–Bampton synthesis technique with the flexible multi-body dynamics approach, and the accuracy of this method was verified by simulating the morphing process according to a certain morphing schedule. Reich, Bowman [11,12] and Scarlett [13] built a time-varying aeroelastic simulating system of the folding wing using the flexible multi-body dynamics method and their inhouse vortex lattice code, and the details about the development of this system were shown in [11,12]. The flight performances during the morphing process of a folding wing aircraft were investigated in [13], and the results indicated the influence of the folding angle on the structural load paths and the fold hinge moments. In this paper, the study is focused on the aeroelastic behavior of a folding wing during the morphing process, and the aerodynamics model is built based on the Doublet Lattice based aerodynamics theory, which is more accurate than the vortex lattice based aerodynamics. The structure model is built by incorporating the Craig–Bampton mode with the flexible multi-body dynamics approach, and the aerodynamics model is built by the Kriging agent model technique with the interpolation of the reduced AIC matrices in the time domain at several different folding angles. The aeroelastic responses of the folding wing during the morphing process are simulated, and then the influences of some morphing parameters, such as the morphing mode (folding or unfolding) and the morphing rate, on the dynamic aeroelastic stability of the folding wing are examined.
2. Structure model 2.1. Substructure mode orthonormalization The folding wing structure can be treated as a flexible multi-body structure, and each flexible body is a substructure, i.e. the central wing (I), the inboard wing (II) and the outboard wing (III) as shown in Fig. 1. The folding angle is defined as the angle between the inboard wing and the x–y plane, and to guarantee the lift performance of the folding wing, the outboard wing remains parallel to the x–y plane at all time. In order to build the dynamic model for the flexible multi-body structure, the Craig–Bampton synthesis technique [14] is employed. For each substructure, the equation of motion is written as follows: " ι #( ι ) " ι ( ) ι #( ι ) €i kii kij ui mii mιij 0 u þ ¼ (1) ι ι ι ι ι € mιji mιjj u f uj kji kjj j j
Fig. 1. Sketch of a folding wing.
218
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
where the subscript i and j denote the interior and interface coordinates of the substructure, respectively. The superscript ι ι ¼ I, II and III denotes the three different substructures. f j is the constraint force applied to the interface coordinates by adjacent substructure. The substructure mode calculated by Craig–Bampton synthesis technique is a combination of normal mode φιk and constraint mode ψιc as follows: " ι # h i φik ψιic ι ι ι φ ψ (2) Φ ¼ c ¼ k 0ιjk Iιjc where φιik is obtained from the eigenvalue Eq. (3). The constraint mode ψιc , which includes the rigid body mode, is calculated from the static equilibrium Eq. (4) as follows: ι ι ι (3) kii ωι2 i mii φik ¼ 0 "
ι
ι
kii
kij
ι kji
ι kjj
#"
ψιic Iιjc
#
" ¼
0 ι f jc
# (4)
Using the Craig–Bampton mode in Eq. (2), the structural displacements in Eq. (1) can be transformed into modal coordinate: uι ¼ Φι qι
(5)
ι
where q is the substructure modal coordinate. Substitution of Eq. (5) into Eq. (1) and premultiplying Eq. (1) by ΦιT yields ( ι ι ι ι m q€ þ k qι ¼ f ι
ι
ι
mι ¼ Φι T mι Φι ; k ¼ Φι T k Φι ; f ¼ Φι T f ι
(6)
ι
ι
where mι , k and f are the mass matrix, stiffness matrix and applied force vector in Eq. (1), respectively. When formulating the dynamics equation of the folding wing by using the flexible multi-body dynamics approach, the rigid body coordinates of each substructure will be defined to describe the rigid body motion of the substructure, and the substructure mode is used to describe the elastic deformation of the substructure. So the rigid body mode included in the substructure mode is redundant and must be eliminated. To eliminate the rigid body mode, the eigenvalue problem in Eq. (7) should be solved: ι
ι
k qι ¼ λ mι qι
(7) ι
ι
ι
From Eq. (7) a transform matrix S is obtained, which is orthogonal with respect to m and k as follows: ι ι SιT mι S ¼ I
(8)
ι
SιT k Sι ¼ Ωι2
where Ωι is a diagonal matrix with natural frequencies of the substructural system obtained by Craig–Bampton synthesis technique. The terms on the right-hand side of Eq. (8) are the corresponding general mass and general stiffness matrices, respectively. The first 6 columns of the transform matrix Sι are the rigid body displacements and rotations, so the first 6 columns of matrix Sι are deleted to eliminate the rigid body mode. Then, the displacement vector in Eq. (1) can be expressed as ι
uι ¼ Φι Sι qι ¼ Φ qι
(9)
ι
where Φ is referred as the orthogonal substructure mode, which is required by the flexible multi-body dynamics method, qι is the corresponding general coordinate. 2.2. Flexible multi-body dynamics model In order to build the flexible multi-body dynamics model, two types of coordinate systems are required. A local coordinate is assigned to each substructure, which translates and rotates with the substructure, and a ground coordinate is defined to describe the motion of the local coordinates. The displacement of any point p on a substructure of the flexible multi-body system can be expressed as follows: (10) rιp ¼ rι0 þ Aι sιp þuιp where rι0 is the displacement vector from the origin of the ground coordinate to the origin of the local coordinate, and Aι is the coordinate transformation matrix from the local coordinate to the ground coordinate. sιp is the position vector from the origin of the local coordinate to point p, and uιp is the deformation vector of point p in the local coordinate.
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
219
According to Eq. (9), Eq. (10) can be rewritten as ι rιp ¼ rι0 þAι sιp þ Φ qι The velocity of point p can be obtained by differentiating Eq. (11) with respect to time _ ι sι þ Φι qι þAι Φι q_ ι vιp ¼ r_ ι0 þ A p
(11)
(12)
With Eq. (12) the kinetic energy of the substructure can be obtained, and then the Lagrange equation can be applied to derive the dynamic equation of the substructure [15]: ι T d ∂T ι T ∂T þ CTXι λ ¼ Q ι (13) ι _ dt ∂X ∂Xι T where T ι is the kinetic energy of the substructure, Xι ¼ Rι ϑι qι is the generalized coordinate of the substructure, R ι and ϑι are a set of Cartesian coordinate and rotational coordinate of the local coordinate, respectively. CXι is the constraint Jacobian matrix. Q ι ¼ Kι Xι þ Q ιe is the generalized force, where Q ιe is the generalized applied force. Expanding the first two terms on the left-hand side of Eq. (13), we obtain ι T
T d ∂T ι T ∂T _ ιT Mι X € ι þM _ ι ∂ 1X _ι _ ιX ¼ Mι X ι ι ι _ dt ∂X ∂X 2 ∂X
(14)
where the last two terms on the right-hand side of Eq. (14) can be called quadratic velocity vector [15], which is resulting from the differentiation of the kinetic energy with respect to time and with respect to the generalized coordinate. Building the dynamic equations of all the substructures according to Eqs. (13) and (14), and considering the constraint conditions, the flexible multi-body dynamics equation can be expressed as follows: h iT _ X _ þKX þCX T λ ¼ Q € þM _X _ 1 ∂MX MX e 2 ∂X C¼0 (15) h iT h iT I II III where M ¼ blockdiag MI MII MIII , K ¼ blockdiag KI KII KIII , X ¼ XI XII XIII , Q e ¼ Q e Q e Q e , C ¼ 0 is the constraint equation. Compared with the dynamics equation of the time invariant system, several terms in Eq. (15) are special: M is a nonlinear mass matrix with respect to the generalized coordinate X, which indicates a typical inertia nonlinearity feature of the flexible multi-body dynamics system. The second and third terms on the left-hand side of Eq. (15) contain the gyroscopic and Coriolis force components [15], and they are also the special terms in the dynamics equation of a time variant system.
Fig. 2. The flow chart of building the aerodynamics model.
220
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
3. Aerodynamics model In order to obtain the aeroelastic response of the folding wing efficiently in the time domain, the aerodynamic force is calculated using the Doublet Lattice method. As the geometry of the folding wing changes during the morphing process, the aerodynamics model must be rebuilt rapidly and efficiently. So the Kriging agent model technique is applied to build the aerodynamics model for the folding wing during morphing process. Before applying the Kriging technique to build the agent model, the order of aerodynamics model should be reduced at first. When calculating the substructure mode by the Craig–Bampton synthesis technique, the order of the substructure mode is larger than the number of the interface coordinate. So if the aerodynamics model is reduced by the substructure modes, the reduced order of the aerodynamics model will be larger than the number of the interface coordinate. Besides, by using this reduced aerodynamics model, the physical velocity and acceleration caused by the rigid body motion must be transformed into substructure modal coordinate firstly, before applying them to calculate the aerodynamic force induced by rigid body motion. In the present study, the reduction of aerodynamics model based on the spline interpolation technique is proposed, and by using this method, the reduced aerodynamics model is independent of the number of the interface coordinate, and the velocity and acceleration can be used directly to calculate the aerodynamic force induced by rigid body motion. The flow chart of building the aerodynamics model is shown in Fig. 2. Several AIC matrices of the folding wing at different folding angles θi are obtained by the Doublet Lattice method, and reduced by the spline interpolation technique (Sections 3.1 and 3.2). Then the minimum state approximation is used to transform these reduced AIC matrices into the time domain (Section 3.3). Finally, the Kriging method is used to build the agent model of the aerodynamics model (Section 3.4). 3.1. Aerodynamics model reduction For the configuration of the folding wing with a specific fixed folding angle θ, the aerodynamics model can be formulated by the Doublet Lattice method [16] as follows: f ¼ q1 AðωÞz
(16)
where q1 is the dynamic pressure, ω is the oscillating frequency, ΑðωÞ is the unsteady AIC matrix, f and z are the force vector and the displacement vector for aerodynamic panels, respectively. The spline interpolation technique, such as the plate spline interpolation technique [17,18] and the beam spline interpolation technique, can be adopted to interpolate the displacement or force through a spline matrix G. Interpolating the
Fig. 3. The flow chart of the genetic optimization process.
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
221
aerodynamics model in Eq. (16) into a set of Ns structural nodes, the relationship in Eq. (17) can be obtained as follows: ( f s ¼ GTs f (17) z ¼ Gs zs where f s and zs are the force vector and the displacement vector for the N s structural nodes, respectively. Substitution of Eq. (17) into Eq. (16) yields the reduced aerodynamics model based on the spline interpolation technique as follows: f s ¼ q1 Ass ðωÞzs where
(18)
Ass ðωÞ ¼ GTs AðωÞGs .
3.2. Optimization of reduced aerodynamics model Although the aerodynamics model can be reduced by the method proposed in Section 3.1, the locations of the interpolated structural nodes must be selected carefully to satisfy the accuracy of the reduced aerodynamics model. In the present study, the genetic algorithm is used to optimize the locations of the interpolated structural nodes. The design variables P are the identification numbers of structural nodes, which can be used to describe the location of each structural node. The optimization problem can be defined mathematically as follows: θn X V f af t V f þ ωf af t ωf Minimize: V ω f f θ θ ¼ θ0 ( node_IDmin oP i o node_IDmax ð1 ri aj r Ns Þ Subject to P aP (19) i
j
where θ is the folding angle, V f and V f af t are the flutter speeds before and after the reduction, respectively. ωf and ωf the flutter frequencies before and after the reduction, respectively. The flow chart of the genetic optimization process can be described as shown in Fig. 3:
af t
are
3.3. Reduced aerodynamics model in the time domain After obtaining the reduced aerodynamics model in frequency domain, the minimum state approximation [19] is employed to transform the aerodynamics model into the time domain, and the unsteady aerodynamics model can be expressed as
2 b b V 1 Ass ðt Þ ¼ A0 þ A1 sþ 2 A2 s2 þ D sI R Es (20) V b V where s is the Laplace variable, V is the air speed, b is the reference length, A0 , A1 , A2 , D, E and R are the coefficient matrices obtained by the minimum state approximation. Substitution of Eq. (20) into Eq. (18) yields the reduced unsteady aerodynamics in the time domain as follows: " #
2 b b V 1 f s ¼ q1 A0 þ A1 s þ 2 A2 s2 zs þ q1 D sI R Eszs (21) V b V 3.4. Agent model for interpolation of aerodynamics Even though the reduced aerodynamics model in the time domain is obtained by the procedures described in Sections 3.1 and 3.3, it will be of poor efficiency to repeat these procedures at each different folding angle during the morphing process. As shown in Eq. (21), the reduced aerodynamics model in the time domain is represented by the coefficient matrices, therefore we can calculate these coefficient matrices at a series of different folding angles, and then apply the Kriging technique [20] to build the agent models for each of the coefficient matrices. So, the reduced aerodynamics model can be rebuilt rapidly by these agent models of the coefficient matrices. Taking the matrix A0 as an example, the procedures to create its agent model are as follows. The coefficient matrices A0θi at different folding angles θi are calculated by the procedures in Sections 3.1 and 3.3, where i ¼ 1; 2; :::; N and N is the number of samples of the folding angle. Thus, the observed data set Y ¼ h iT A0θ1 ðm; nÞ A0θ2 ðm; nÞ ::: A0θN ðm; nÞ can be obtained, where A0 ðm; nÞ is the element value in the mth row and the nth θi column of matrix A0θi . The mathematical model of the Kriging technique can be expressed as T
yðθÞ ¼ f ðθÞβ þzðθÞ
(22)
222
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
h iT h iT where f ðθÞ ¼ f 1 ðθÞ f 2 ðθÞ ::: f p ðθÞ , β ¼ β1 β2 ::: βp . f ðθÞ is the regression function, which is a polynomial of θ, β is the coefficient of f ðθÞ, zðθÞ is a stochastic error with mean zero and its covariance is Cov½zðθm Þ; zðθn Þ ¼ σ 2 Rðξ; θm ; θn Þ
ð 1 r m; n rNÞ
(23)
where σ 2 is the variance of zðθÞ, Rðξ; θm ; θn Þ is the correlation function to represent the correlation between the sample point θm and θn , ξ is the parameter of this function. And the Gaussian correlation function is used in the present study, which is expressed as follows: (24) Rðξ; θm ; θn Þ ¼ exp ξ U ðθm θn Þ2 Based on the assumption of Kriging technique, the response at point θ is evaluated by the linear combination of the observed data set Y y^ ðθÞ ¼ cT ðθÞY
(25)
E y^ ðθÞ yðθÞ ¼ E cT ðθÞY yðθÞ ¼ 0
(26)
Subject to the unbiasedness constraint
Combining Eqs. (22) and (26), cðθÞ can be obtained. Substituting cðθÞ into Eq. (25), the following response can be obtained: T y^ ðθÞ ¼ f ðθÞβ^ þrT ðθÞR 1 Y Fβ^ (27) where 1 FT R 1 Y β^ ¼ FT R 1 F rðθÞ ¼ Rðξ; θ; θ1 Þ 2 6 R¼4
Rðξ; θ; θ2 Þ
Rðξ; θ; θN Þ
:::
R ðξ; θ1 ; θ1 Þ
:::
Rðξ; θ1 ; θN Þ
⋮
⋱
⋮
R ðξ; θN ; θ1 Þ
:::
Rðξ; θN ; θN Þ
h F ¼ f ðθ 1 Þ
f ðθ2 Þ
:::
f ðθ N Þ
T
3 7 5
iT
^ 0 ðθÞ can be obtained. After calculating the response of every element of matrix A0 by Eq. (27), the estimated matrix A
Fig. 4. The geometries of the folding wing.
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
223
The agent models of the other coefficient matrices can be obtained through the same procedure. Then the estimated reduced aerodynamics model of the folding wing in the time domain can be represented as follows: " #
1 2 ^f s ðθÞ ¼ q A ^ 0 ðθ Þ þ b A ^ 2 ðθÞs2 zs þ q D ^ 1 ðθÞs þ b A ^ ðθ Þ ^ ðθÞ sI V R (28) E^ ðθÞszs 1 1 V b V2 Combining the structure model of Eq. (15) and the aerodynamics model of Eq. (28), the aeroelastic response of the folding wing can be simulated by time marching technique. Starting from an initial state, the motion responses of the reduced structural nodes and the folding angle at the current time step are obtained by solving Eq. (15). Then substituting these results into Eq. (28), the unsteady aerodynamic force at the reduced structural nodes is calculated. Substituting the aerodynamic force back into Eq. (15), the motion responses of the reduced structural nodes and the folding angle at the next time step can be obtained. By repeating these steps until sufficient response time history is achieved, the aeroelastic response of the folding wing can be obtained.
4. Numerical examples The geometries of the folding wing in the present study are shown in Fig. 4. The folding wing contains three substructures, i.e. the central wing, the inboard wing and the outboard wing. In the present study, the structure of the folding wing is represented by three uniform plates and the thicknesses of them are 5 cm, 2 cm and 2 cm. These substructures are made of aluminum, and are connected by rotation hinge springs, which are on the hinge lines as shown in Fig. 4 and marked by circles. The rotation stiffness of each spring is 800 kN m/rad.
Fig. 5. The comparison of natural frequencies at different folding angles.
Fig. 6. The locations of the reduced nodes of the aerodynamics model.
224
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
^ 0 at different folding angles. Fig. 7. The comparisons of some element values of A0 and A
^ 0 when θ¼82.5°, (a) matrix A0 and (b) matrix A ^ 0. Fig. 8. The comparison of element values of A0 and A
4.1. Examination of structure model and aerodynamics model The flexible multi-body dynamics model of the folding wing is built by the procedures proposed in Sections 2.1 and 2.2, and the natural frequencies of this model at different fixed folding angles are compared with those obtained from the finite element model, which is built by MSC.Nastran. The comparison of results in Fig. 5 shows a good agreement between the natural frequencies obtained by using these two approaches. By reducing the aerodynamics model with the spline interpolation technique, the unsteady AIC matrix is reduced into matrix Ass ðωÞ with the dimension of 24 24, and the locations of these 24 structural nodes are shown in Fig. 6. Then, the minimum state approximation is adopted to calculate the coefficient matrices at some specific fixed folding angles, and the Kriging agent model technique is used to build the aerodynamics model. A series of folding angles θ ¼ 01 151 301 451 601 751 901 1051 1201 are selected as the sample points to ^ 0 at different build the agent model. Fig. 7 shows the comparisons of some element values of A0 and its agent model A folding angles. Additionally, in order to further verify the accuracy of the agent model obtained by the Kriging technique, the matrices ^ 0 at the folding angle of 82:51, which is not included in the sample points, are compared graphically in Fig. 8. The A0 and A ^ 0 are shown in Fig. 9, and the errors are within 5%. percent errors of element values of A0 and A 4.2. Flutter analysis Based on the structure model and aerodynamics model obtained in Section 4.1, the aeroelastic responses of the folding wing at a specific fixed folding angle can be simulated under a series of air speeds. The flutter speed can be obtained at the air speed when the aeroelastic response is harmonic, and the flutter frequency can be calculated from this response. By these procedures, the flutter characteristics of the folding wing at different folding angles are obtained. These flutter
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
225
^ 0 when θ¼ 82.5°. Fig. 9. The percent errors of element values of A0 and A
Fig. 10. The comparison of the flutter characteristics obtained by different methods at different folding angles.
Fig. 11. The dynamic aeroelastic response of the folding wing with folding rate of 51/s.
characteristics are compared with those obtained by the g method using the finite element models of the folding wing structure and the Doublet Lattice based aerodynamics models corresponding to the configurations of different folding angles. The comparison between these results is shown in Fig. 10, and as expected, the results agree well. Fig. 10 also shows that the flutter speed decreases as the folding angle is increased from 01 to 601, whereas the flutter speed increases as the folding angle is increased from 601 to 1201. The lowest flutter speed is about 128 m/s when the folding angle is about 601. It also indicates that one simply needs to fly below the lowest flutter speed while morphing, but then there is some additional speed margin to the flight envelope based on the folding angle, because the flutter speed is shifting according to the folding angle. What if the flight speed is larger than the lowest flutter speed while morphing? A baseline speed of 147 m/s is chosen to examine the aeroelastic response of the folding wing during the morphing process. It is indicated by the solid line in Fig. 10 and this value is chosen because it is expected the wing will flutter during the morphing process as the flutter speed between 331 and 771 drops below this value.
226
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
Fig. 12. The dynamic aeroelastic response of the folding wing with folding rate of 101/s.
Fig. 13. The dynamic aeroelastic response of the folding wing with folding rate of 201/s.
In the simulation, the angle of attack of the folding wing is 01. The aerodynamic force during the morphing process consists of two parts: the unsteady aerodynamic force induced by the vibration of the folding wing, and the unsteady aerodynamic force induced by the rigid body motion of the folding wing. In order to examine the influence of the aerodynamic force induced by the rigid body motion on the aeroelastic characteristics, two cases are studied for the aeroelastic response during the morphing process: Case 1. Aeroelastic response of the folding wing at the speed of 147 m/s considering the aerodynamic force induced by the rigid body motion. Case 2. Aeroelastic response of the folding wing at the speed of 147 m/s without considering the aerodynamic force induced by the rigid body motion. The deflection responses at the trailing edge tip of the outboard wing during the folding processes with different folding rates are shown in Figs. 11–13. The dynamic aeroelastic stability obtained by the quasi-steady aeroelastic analysis is also shown in the figures, so that it can be compared directly with those obtained by the dynamic aeroelastic analysis during the morphing process. It can be seen from Fig. 10 that when the air speed is 147 m/s, the wing loses its dynamic aeroelastic stability when the folding angles are in the range of 331 to 771. Thus the dynamic aeroelastic stability of the folding wing at different folding angles can be divided into three regions: two stable regions corresponding to the folding angle smaller than 331 (S1) and larger than 771 (S2), and one unstable region (U) corresponding to the folding angle ranging from 331 to 771. At the same folding rate, the difference of the responses between Case 1 and Case 2 indicates that the aerodynamic force induced by the rigid body motion just changes the static equilibrium position of the responses and has no influence on the dynamic aeroelastic stability, i.e. the influence of this aerodynamic force on the aeroelastic characteristics is similar to the influence of the gust aerodynamic force. And it also shows that when the folding angle is increased from 01 to 1201, the shifting value of the static equilibrium position of the responses of Case 1 changes from negative to positive, and the shifting value is zero when the folding angle is 901. This phenomenon can be explained as follows. The shifting of the static equilibrium position of the response is caused by the aerodynamic force induced by the rigid body motion, i.e., the motion velocity and acceleration normal to the wing surface caused by the rigid body motion. As shown in Fig. 14, for the outboard wing, when the folding angle is increased from 01 to 1201, the direction of the motion velocity component normal to the wing surface is changed from upward to downward, and becomes zero when the folding angle is 901. And because the direction of aerodynamic force induced by the rigid body motion is opposite to the direction of
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
227
Fig. 14. The direction of the motion velocity component during the folding process.
Fig. 15. The dynamic aeroelastic response of the folding wing from 701 to 901 with folding rate of 0:51/s.
Fig. 16. The dynamic aeroelastic response of the folding wing with unfolding rate of 51/s.
the motion velocity component, thus, when the folding angle is increased from 01 to 901, the direction of the aerodynamic force induced by the rigid body motion is downward, and it becomes upward when the folding angle is increased from 901 to 1201, so that the shifting value of the static equilibrium position changes from negative to positive, and the shifting value is zero when the folding angle is 901. Figs. 11–13 also indicate that the dynamic aeroelastic stabilities obtained by the aeroelastic responses of Case 1 are different from that obtained by the quasi-steady aeroelastic analysis. During the folding process, the dynamic aeroelastic stability changes from stable state to unstable state when the folding angle is slightly larger than 331, and with the increasing of the folding rate, the critical folding angle at which the dynamic aeroelastic stability is lost becomes much larger than 331. The similar conclusion can be found in [21], where the aeroelastic characteristics of a variable-span wing during the morphing process were investigated. Besides, according to the aeroelastic stability obtained by the quasi-steady aeroelastic analysis, the amplitude of the response should increase when the folding angle is increased from 331 to 771, and decrease when the folding angle is increased from 771 to 1201. The responses in Figs. 11–13 indicate that when the folding angle is changed into a region nearby 771, a dip of the amplitude of the response appears, and the width of the region increases with the increasing of the folding rate. The two aforementioned phenomena are caused by the motions of the inboard and the outboard wings. According to Eq. _ and (15), the quadratic velocity vector is generated because of the motions of the substructures, and can be represented by X
228
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
Fig. 17. The dynamic aeroelastic response of the folding wing with unfolding rate of 101/s.
Fig. 18. The dynamic aeroelastic response of the folding wing with unfolding rate of 201/s.
Fig. 19. The dynamic aeroelastic response of the folding wing from 401 to 251 with unfolding rate of 0:51/s.
its coefficient matrices. If we treat the coefficient matrices as damping matrices, it means that the motions of the inboard and the outboard wings will induce additional damping into the aeroelastic system of the folding wing. The additional damping will influence the aeroelastic characteristics of the folding wing, especially when the folding angle is approaching to the specific folding angle at which the flutter speed equals to the air speed. Moreover, this influence increases with the increasing of the folding rate. From the conclusion in the above paragraph, it is definite that if the folding rate is very small, the dynamic aeroelastic stability of the folding wing must be similar to that obtained by the quasi-steady aeroelastic analysis. Considering the aerodynamic force induced by the rigid body motion, when the folding wing is folding from 701 to 901 with the folding rate of 0:51/s, the deflection response is shown in Fig. 15, and the dynamic aeroelastic stability of the folding wing agrees well with that obtained by the quasi-steady aeroelastic analysis. The unfolding processes with different unfolding rates are also simulated, and the deflection responses at the trailing edge tip of the outboard wing are shown in Figs. 16–18 , respectively. Both the deflection responses of Cases 1 and 2 are
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
229
displayed in these figures, and the difference of the responses between Cases 1 and 2 show the same conclusion about the influence of the aerodynamic force induced by the rigid body motion on the aeroelastic characteristics of the folding wing. The aeroelastic responses of Case 1 as displayed in Figs. 16–18 show that during the unfolding process, the dynamic aeroelastic stability changes from unstable state to stable state when the folding angle is less than 331, and this angle becomes much less than 331 with the increasing of the unfolding rate. While as shown in Figs. 11–13, the folding angle at which the dynamic aeroelastic stability changes from stable state to unstable state is larger than 331, this phenomenon indicates that the folding and unfolding processes have opposite influences on the dynamic aeroelastic stability of the folding wing, and the influences become much more significant with the increasing of the folding and unfolding rates. The similar conclusion can also be found in [21]. The deflection response of the folding wing from 401 to 251 with the unfolding rate of 0:51/s is shown in Fig. 19, and the dynamic aeroelastic stability of the folding wing also agrees well with that obtained by the quasi-steady aeroelastic analysis.
5. Conclusions To study the aeroelastic characteristics of a folding wing during the morphing process, the structure model is built by incorporating the Craig–Bampton modal synthesis technique with the flexible multi-body dynamics approach. A method based on the Kriging agent model technique is proposed to build the aerodynamics model, which is proved to be accurate and can be used to calculate the aerodynamic force of the folding wing during the morphing process efficiently. The aeroelastic responses during the folding and unfolding processes are obtained by numerical simulation. The results show that the morphing process has an influence on the dynamic aeroelastic stability of the folding wing, which makes the dynamic aeroelastic stabilities of the folding wing different from that obtained by the quasi-steady aeroelastic analysis, and the influences of the folding and unfolding processes on the dynamic aeroelastic stability are also different. Furthermore, the influence of the morphing process on the dynamic aeroelastic stability becomes much more significant with the increasing of the morphing rate. Acknowledgments This work was founded by the National Natural Science Foundation of China (Grants no. 11472216), 111 Project of China (B07050).
References [1] M.H. Love, P.S. Zink, R.L. Stroud, et al., Demonstration of morphing technology through ground and wind tunnel tests, AIAA Paper 2007-1729, May 2007. [2] M. Snyder, G. Frank, B. Sanders, Aeroelastic analysis of a morphing Z-wing configuration, Aerospace Flutter and Dynamics Council Meeting, May 2004. [3] J. Dunn, L. Horta, T. Ivanco, et al., NASA contributions to DARPA MAS program, Aerospace Flutter and Dynamics Council Meeting, May 2004. [4] M. Snyder, B. Sanders, F. Eastep, et al., Sensitivity of flutter to fold orientation and spring stiffness of a simple folding wing, Proceedings of the International Forum on Aeroelasticity and Structural Dynamics 2005, June–July 2005. [5] D.H. Lee, T.A. Weisshaar, Aeroelastic studies on a folding wing configuration, AIAA Paper 2005-1996, April 2005. [6] Y. Zhao, H. Hu, Parameterized aeroelastic modeling and flutter analysis for a folding wing, Journal of Sound and Vibration 331 (2012) 308–324. [7] D.H. Lee, P.C. Chen, Nonlinear aeroelastic studies on a folding wing configuration with free-play hinge nonlinearity, AIAA Paper 2006-1734, May 2006. [8] D. Tang, E.H. Dowell, Theoretical and experimental aeroelastic study for folding wing structures, Journal of Aircraft 45 (2008) 1136–1147. [9] P.J. Attar, D. Tang, E.H. Dowell, Nonlinear aeroelastic study for folding wing structures, AIAA Journal 48 (2010) 2187–2195. [10] Y. Zhao, H. Hu, Prediction of transient responses of a folding wing during the morphing process, Aerospace Science and Technology 24 (2013) 89–94. [11] G.W. Reich, J.C. Bowman J., B. Sanders, et al., Development of an integrated aeroelastic multibody morphing simulation tool, AIAA Paper 2006-1892, May 2006. [12] J.C. Bowman, G.W. Reich, B. Sanders, et al., Simulation tool for analyzing complex shape-changing mechanisms in aircraft, AIAA Paper 2006-6727, August 2006. [13] J.N. Scarlett, R.A. Canfield, B. Sanders, Multibody dynamic aeroelastic simulation of a folding wing aircraft, AIAA Paper 2006-2135, May 2006. [14] R.R. Craig, Coupling of substructures for dynamic analyses, AIAA Journal 6 (1968) 1313–1319. [15] A.A. Shabana, Dynamics of Multibody System, Cambridge University Press, New York, 1998. [16] W.P. Rodden, J.P. Giesing, T.P. Kalman, New method for nonplanar configurations, AGARD Conference Proceedings, 1971. [17] R.L. Harder, R.N. Desmarais, Interpolation using surface splines, Journal of Aircraft 9 (1972) 189–191. [18] J. Duchon, Splines minimizing rotation-invariant semi-norms in sobolev spaces, Constructive Theory of Functions of Several Variables 571 (1977) 85–100. [19] M. Karpel, Procedures and models for aeroservoelastic analysis and design, Journal of Applied Mathematics and Mechanics 81 (2001) 579–592. [20] J. Sacks, W.J. Welch, T.J. Mitchell, et al., Design and analysis of computer experiments, Statistical Science 4 (1989) 409–435. [21] R. Huang, Z. Qiu, Transient aeroelastic responses and flutter analysis of a variable-span wing during the morphing process, Chinese Journal of Aeronautics 26 (6) (2013) 1430–1438.
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
Aeroelastic study for folding wing during the morphing process Wei Hu, Zhichun Yang n, Yingsong Gu School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China
a r t i c l e in f o
abstract
Article history: Received 15 April 2015 Received in revised form 14 October 2015 Accepted 29 November 2015 Handling editor: W. Lacarbonara Available online 14 December 2015
This paper focuses on the aeroelastic characteristics of a folding wing during the morphing process. The folding wing structure is modeled by using the flexible multi-body dynamics approach, and an efficient method is proposed to calculate the aerodynamic force of the folding wing during the morphing process. The aerodynamic influence coefficient (AIC) matrices at different folding angles are obtained by the Doublet Lattice based aerodynamics theory, and then the orders of these AIC matrices are reduced by the spline interpolation technique. Through the minimum state approximation, the reduced AIC matrices are described as rational functions in the Laplace domain. Then the Kriging agent model technique is used to interpolate the coefficient matrices of the rational functions obtained from several different folding angles and to build the aerodynamics model in the time domain. At some different folding angles, the element values of the coefficient matrices before and after the interpolation are compared to verify the accuracy of the aerodynamics model, and then the aeroelastic responses of the folding wing during its morphing processes are simulated. The results demonstrate that the folding and unfolding processes have opposite influences on the dynamic aeroelastic stability of the folding wing, and the influences become much more significant with the increasing of the folding and unfolding rates. When the folding wing is morphing with a very slow rate, the dynamic aeroelastic stability will be similar to that obtained by the quasi-steady aeroelastic analysis. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Folding wing Aeroelasticity Flexible multi-body Time varying system Morphing wing
1. Introduction Since the beginning of the Morphing Aircraft Structures (MAS) Program, some novel morphing aircraft concepts have been proposed. Morphing aircraft holds several advantages over conventional aircraft designs including multi-mission capability, stemming from its ability of varying wing geometries during the flight. One of the famous morphing aircraft concepts is the folding wing aircraft, proposed by Lockheed Martin Corporation [1]. In recent years, many studies have been conducted on the aeroelastic characteristics of the folding wing aircraft. Several investigations [2–4] were performed to study the influences of structural parameters, such as the folding angle and the hinge spring stiffness, on the aeroelastic characteristics. Weisshaar and Lee [5] built a high fidelity aeroelastics model of the folding wing, and the results showed that the flutter dynamic pressure increases with the increasing of the inboard wing folding angle. For the extended wing configuration, the flutter dynamic pressure is much more sensitive to the changes in n
Corresponding author. Tel./fax: þ86 29 88460461. E-mail addresses: [email protected], [email protected] (Z. Yang).
http://dx.doi.org/10.1016/j.jsv.2015.11.043 0022-460X/& 2015 Elsevier Ltd. All rights reserved.
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
217
the inboard hinge stiffness, while for the folded wing configuration, the flutter dynamic pressure is much more sensitive to the changes in the outboard hinge stiffness. Zhao and Hu [6] proposed a parameterized aeroelastic modeling procedure, and the flutter characteristics of a folding wing with different configurations were investigated. Some studies on the nonlinear aeroelastic characteristics of the folding wing were also carried out. Lee and Chen [7] investigated the nonlinear aeroelastic characteristics of a folding wing with free-play hinge nonlinearity. In their study, even with a small free-play nonlinearity, the limit cycle oscillation (LCO) was observed when the folding angle is 0–301. And the LCO phenomenon would disappear when the folding angle is larger than 301. Tang, Dowell and Attar [8,9] investigated the influence of geometric nonlinearity on aeroelastic behaviors of a folding wing through experimental and computational methods, and the responses obtained by these two methods agreed well with each other. The result showed that the LCO appeared at air speed beyond the linear flutter speed of the folding wing, and the LCO phenomenon was changed at different folding angles. It is noted that all of these studies mentioned above were analyzed in the quasi-steady condition with fixed folding angles, so the results were only suitable for the folding wing during very slow morphing process. Few studies have explored whether these results are suitable for the folding wing with a rather rapid morphing rate in normal morphing process. Zhao and Hu [10] have developed the flexible multi-body dynamics formulation by combining the Craig–Bampton synthesis technique with the flexible multi-body dynamics approach, and the accuracy of this method was verified by simulating the morphing process according to a certain morphing schedule. Reich, Bowman [11,12] and Scarlett [13] built a time-varying aeroelastic simulating system of the folding wing using the flexible multi-body dynamics method and their inhouse vortex lattice code, and the details about the development of this system were shown in [11,12]. The flight performances during the morphing process of a folding wing aircraft were investigated in [13], and the results indicated the influence of the folding angle on the structural load paths and the fold hinge moments. In this paper, the study is focused on the aeroelastic behavior of a folding wing during the morphing process, and the aerodynamics model is built based on the Doublet Lattice based aerodynamics theory, which is more accurate than the vortex lattice based aerodynamics. The structure model is built by incorporating the Craig–Bampton mode with the flexible multi-body dynamics approach, and the aerodynamics model is built by the Kriging agent model technique with the interpolation of the reduced AIC matrices in the time domain at several different folding angles. The aeroelastic responses of the folding wing during the morphing process are simulated, and then the influences of some morphing parameters, such as the morphing mode (folding or unfolding) and the morphing rate, on the dynamic aeroelastic stability of the folding wing are examined.
2. Structure model 2.1. Substructure mode orthonormalization The folding wing structure can be treated as a flexible multi-body structure, and each flexible body is a substructure, i.e. the central wing (I), the inboard wing (II) and the outboard wing (III) as shown in Fig. 1. The folding angle is defined as the angle between the inboard wing and the x–y plane, and to guarantee the lift performance of the folding wing, the outboard wing remains parallel to the x–y plane at all time. In order to build the dynamic model for the flexible multi-body structure, the Craig–Bampton synthesis technique [14] is employed. For each substructure, the equation of motion is written as follows: " ι #( ι ) " ι ( ) ι #( ι ) €i kii kij ui mii mιij 0 u þ ¼ (1) ι ι ι ι ι € mιji mιjj u f uj kji kjj j j
Fig. 1. Sketch of a folding wing.
218
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
where the subscript i and j denote the interior and interface coordinates of the substructure, respectively. The superscript ι ι ¼ I, II and III denotes the three different substructures. f j is the constraint force applied to the interface coordinates by adjacent substructure. The substructure mode calculated by Craig–Bampton synthesis technique is a combination of normal mode φιk and constraint mode ψιc as follows: " ι # h i φik ψιic ι ι ι φ ψ (2) Φ ¼ c ¼ k 0ιjk Iιjc where φιik is obtained from the eigenvalue Eq. (3). The constraint mode ψιc , which includes the rigid body mode, is calculated from the static equilibrium Eq. (4) as follows: ι ι ι (3) kii ωι2 i mii φik ¼ 0 "
ι
ι
kii
kij
ι kji
ι kjj
#"
ψιic Iιjc
#
" ¼
0 ι f jc
# (4)
Using the Craig–Bampton mode in Eq. (2), the structural displacements in Eq. (1) can be transformed into modal coordinate: uι ¼ Φι qι
(5)
ι
where q is the substructure modal coordinate. Substitution of Eq. (5) into Eq. (1) and premultiplying Eq. (1) by ΦιT yields ( ι ι ι ι m q€ þ k qι ¼ f ι
ι
ι
mι ¼ Φι T mι Φι ; k ¼ Φι T k Φι ; f ¼ Φι T f ι
(6)
ι
ι
where mι , k and f are the mass matrix, stiffness matrix and applied force vector in Eq. (1), respectively. When formulating the dynamics equation of the folding wing by using the flexible multi-body dynamics approach, the rigid body coordinates of each substructure will be defined to describe the rigid body motion of the substructure, and the substructure mode is used to describe the elastic deformation of the substructure. So the rigid body mode included in the substructure mode is redundant and must be eliminated. To eliminate the rigid body mode, the eigenvalue problem in Eq. (7) should be solved: ι
ι
k qι ¼ λ mι qι
(7) ι
ι
ι
From Eq. (7) a transform matrix S is obtained, which is orthogonal with respect to m and k as follows: ι ι SιT mι S ¼ I
(8)
ι
SιT k Sι ¼ Ωι2
where Ωι is a diagonal matrix with natural frequencies of the substructural system obtained by Craig–Bampton synthesis technique. The terms on the right-hand side of Eq. (8) are the corresponding general mass and general stiffness matrices, respectively. The first 6 columns of the transform matrix Sι are the rigid body displacements and rotations, so the first 6 columns of matrix Sι are deleted to eliminate the rigid body mode. Then, the displacement vector in Eq. (1) can be expressed as ι
uι ¼ Φι Sι qι ¼ Φ qι
(9)
ι
where Φ is referred as the orthogonal substructure mode, which is required by the flexible multi-body dynamics method, qι is the corresponding general coordinate. 2.2. Flexible multi-body dynamics model In order to build the flexible multi-body dynamics model, two types of coordinate systems are required. A local coordinate is assigned to each substructure, which translates and rotates with the substructure, and a ground coordinate is defined to describe the motion of the local coordinates. The displacement of any point p on a substructure of the flexible multi-body system can be expressed as follows: (10) rιp ¼ rι0 þ Aι sιp þuιp where rι0 is the displacement vector from the origin of the ground coordinate to the origin of the local coordinate, and Aι is the coordinate transformation matrix from the local coordinate to the ground coordinate. sιp is the position vector from the origin of the local coordinate to point p, and uιp is the deformation vector of point p in the local coordinate.
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
219
According to Eq. (9), Eq. (10) can be rewritten as ι rιp ¼ rι0 þAι sιp þ Φ qι The velocity of point p can be obtained by differentiating Eq. (11) with respect to time _ ι sι þ Φι qι þAι Φι q_ ι vιp ¼ r_ ι0 þ A p
(11)
(12)
With Eq. (12) the kinetic energy of the substructure can be obtained, and then the Lagrange equation can be applied to derive the dynamic equation of the substructure [15]: ι T d ∂T ι T ∂T þ CTXι λ ¼ Q ι (13) ι _ dt ∂X ∂Xι T where T ι is the kinetic energy of the substructure, Xι ¼ Rι ϑι qι is the generalized coordinate of the substructure, R ι and ϑι are a set of Cartesian coordinate and rotational coordinate of the local coordinate, respectively. CXι is the constraint Jacobian matrix. Q ι ¼ Kι Xι þ Q ιe is the generalized force, where Q ιe is the generalized applied force. Expanding the first two terms on the left-hand side of Eq. (13), we obtain ι T
T d ∂T ι T ∂T _ ιT Mι X € ι þM _ ι ∂ 1X _ι _ ιX ¼ Mι X ι ι ι _ dt ∂X ∂X 2 ∂X
(14)
where the last two terms on the right-hand side of Eq. (14) can be called quadratic velocity vector [15], which is resulting from the differentiation of the kinetic energy with respect to time and with respect to the generalized coordinate. Building the dynamic equations of all the substructures according to Eqs. (13) and (14), and considering the constraint conditions, the flexible multi-body dynamics equation can be expressed as follows: h iT _ X _ þKX þCX T λ ¼ Q € þM _X _ 1 ∂MX MX e 2 ∂X C¼0 (15) h iT h iT I II III where M ¼ blockdiag MI MII MIII , K ¼ blockdiag KI KII KIII , X ¼ XI XII XIII , Q e ¼ Q e Q e Q e , C ¼ 0 is the constraint equation. Compared with the dynamics equation of the time invariant system, several terms in Eq. (15) are special: M is a nonlinear mass matrix with respect to the generalized coordinate X, which indicates a typical inertia nonlinearity feature of the flexible multi-body dynamics system. The second and third terms on the left-hand side of Eq. (15) contain the gyroscopic and Coriolis force components [15], and they are also the special terms in the dynamics equation of a time variant system.
Fig. 2. The flow chart of building the aerodynamics model.
220
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
3. Aerodynamics model In order to obtain the aeroelastic response of the folding wing efficiently in the time domain, the aerodynamic force is calculated using the Doublet Lattice method. As the geometry of the folding wing changes during the morphing process, the aerodynamics model must be rebuilt rapidly and efficiently. So the Kriging agent model technique is applied to build the aerodynamics model for the folding wing during morphing process. Before applying the Kriging technique to build the agent model, the order of aerodynamics model should be reduced at first. When calculating the substructure mode by the Craig–Bampton synthesis technique, the order of the substructure mode is larger than the number of the interface coordinate. So if the aerodynamics model is reduced by the substructure modes, the reduced order of the aerodynamics model will be larger than the number of the interface coordinate. Besides, by using this reduced aerodynamics model, the physical velocity and acceleration caused by the rigid body motion must be transformed into substructure modal coordinate firstly, before applying them to calculate the aerodynamic force induced by rigid body motion. In the present study, the reduction of aerodynamics model based on the spline interpolation technique is proposed, and by using this method, the reduced aerodynamics model is independent of the number of the interface coordinate, and the velocity and acceleration can be used directly to calculate the aerodynamic force induced by rigid body motion. The flow chart of building the aerodynamics model is shown in Fig. 2. Several AIC matrices of the folding wing at different folding angles θi are obtained by the Doublet Lattice method, and reduced by the spline interpolation technique (Sections 3.1 and 3.2). Then the minimum state approximation is used to transform these reduced AIC matrices into the time domain (Section 3.3). Finally, the Kriging method is used to build the agent model of the aerodynamics model (Section 3.4). 3.1. Aerodynamics model reduction For the configuration of the folding wing with a specific fixed folding angle θ, the aerodynamics model can be formulated by the Doublet Lattice method [16] as follows: f ¼ q1 AðωÞz
(16)
where q1 is the dynamic pressure, ω is the oscillating frequency, ΑðωÞ is the unsteady AIC matrix, f and z are the force vector and the displacement vector for aerodynamic panels, respectively. The spline interpolation technique, such as the plate spline interpolation technique [17,18] and the beam spline interpolation technique, can be adopted to interpolate the displacement or force through a spline matrix G. Interpolating the
Fig. 3. The flow chart of the genetic optimization process.
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
221
aerodynamics model in Eq. (16) into a set of Ns structural nodes, the relationship in Eq. (17) can be obtained as follows: ( f s ¼ GTs f (17) z ¼ Gs zs where f s and zs are the force vector and the displacement vector for the N s structural nodes, respectively. Substitution of Eq. (17) into Eq. (16) yields the reduced aerodynamics model based on the spline interpolation technique as follows: f s ¼ q1 Ass ðωÞzs where
(18)
Ass ðωÞ ¼ GTs AðωÞGs .
3.2. Optimization of reduced aerodynamics model Although the aerodynamics model can be reduced by the method proposed in Section 3.1, the locations of the interpolated structural nodes must be selected carefully to satisfy the accuracy of the reduced aerodynamics model. In the present study, the genetic algorithm is used to optimize the locations of the interpolated structural nodes. The design variables P are the identification numbers of structural nodes, which can be used to describe the location of each structural node. The optimization problem can be defined mathematically as follows: θn X V f af t V f þ ωf af t ωf Minimize: V ω f f θ θ ¼ θ0 ( node_IDmin oP i o node_IDmax ð1 ri aj r Ns Þ Subject to P aP (19) i
j
where θ is the folding angle, V f and V f af t are the flutter speeds before and after the reduction, respectively. ωf and ωf the flutter frequencies before and after the reduction, respectively. The flow chart of the genetic optimization process can be described as shown in Fig. 3:
af t
are
3.3. Reduced aerodynamics model in the time domain After obtaining the reduced aerodynamics model in frequency domain, the minimum state approximation [19] is employed to transform the aerodynamics model into the time domain, and the unsteady aerodynamics model can be expressed as
2 b b V 1 Ass ðt Þ ¼ A0 þ A1 sþ 2 A2 s2 þ D sI R Es (20) V b V where s is the Laplace variable, V is the air speed, b is the reference length, A0 , A1 , A2 , D, E and R are the coefficient matrices obtained by the minimum state approximation. Substitution of Eq. (20) into Eq. (18) yields the reduced unsteady aerodynamics in the time domain as follows: " #
2 b b V 1 f s ¼ q1 A0 þ A1 s þ 2 A2 s2 zs þ q1 D sI R Eszs (21) V b V 3.4. Agent model for interpolation of aerodynamics Even though the reduced aerodynamics model in the time domain is obtained by the procedures described in Sections 3.1 and 3.3, it will be of poor efficiency to repeat these procedures at each different folding angle during the morphing process. As shown in Eq. (21), the reduced aerodynamics model in the time domain is represented by the coefficient matrices, therefore we can calculate these coefficient matrices at a series of different folding angles, and then apply the Kriging technique [20] to build the agent models for each of the coefficient matrices. So, the reduced aerodynamics model can be rebuilt rapidly by these agent models of the coefficient matrices. Taking the matrix A0 as an example, the procedures to create its agent model are as follows. The coefficient matrices A0θi at different folding angles θi are calculated by the procedures in Sections 3.1 and 3.3, where i ¼ 1; 2; :::; N and N is the number of samples of the folding angle. Thus, the observed data set Y ¼ h iT A0θ1 ðm; nÞ A0θ2 ðm; nÞ ::: A0θN ðm; nÞ can be obtained, where A0 ðm; nÞ is the element value in the mth row and the nth θi column of matrix A0θi . The mathematical model of the Kriging technique can be expressed as T
yðθÞ ¼ f ðθÞβ þzðθÞ
(22)
222
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
h iT h iT where f ðθÞ ¼ f 1 ðθÞ f 2 ðθÞ ::: f p ðθÞ , β ¼ β1 β2 ::: βp . f ðθÞ is the regression function, which is a polynomial of θ, β is the coefficient of f ðθÞ, zðθÞ is a stochastic error with mean zero and its covariance is Cov½zðθm Þ; zðθn Þ ¼ σ 2 Rðξ; θm ; θn Þ
ð 1 r m; n rNÞ
(23)
where σ 2 is the variance of zðθÞ, Rðξ; θm ; θn Þ is the correlation function to represent the correlation between the sample point θm and θn , ξ is the parameter of this function. And the Gaussian correlation function is used in the present study, which is expressed as follows: (24) Rðξ; θm ; θn Þ ¼ exp ξ U ðθm θn Þ2 Based on the assumption of Kriging technique, the response at point θ is evaluated by the linear combination of the observed data set Y y^ ðθÞ ¼ cT ðθÞY
(25)
E y^ ðθÞ yðθÞ ¼ E cT ðθÞY yðθÞ ¼ 0
(26)
Subject to the unbiasedness constraint
Combining Eqs. (22) and (26), cðθÞ can be obtained. Substituting cðθÞ into Eq. (25), the following response can be obtained: T y^ ðθÞ ¼ f ðθÞβ^ þrT ðθÞR 1 Y Fβ^ (27) where 1 FT R 1 Y β^ ¼ FT R 1 F rðθÞ ¼ Rðξ; θ; θ1 Þ 2 6 R¼4
Rðξ; θ; θ2 Þ
Rðξ; θ; θN Þ
:::
R ðξ; θ1 ; θ1 Þ
:::
Rðξ; θ1 ; θN Þ
⋮
⋱
⋮
R ðξ; θN ; θ1 Þ
:::
Rðξ; θN ; θN Þ
h F ¼ f ðθ 1 Þ
f ðθ2 Þ
:::
f ðθ N Þ
T
3 7 5
iT
^ 0 ðθÞ can be obtained. After calculating the response of every element of matrix A0 by Eq. (27), the estimated matrix A
Fig. 4. The geometries of the folding wing.
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
223
The agent models of the other coefficient matrices can be obtained through the same procedure. Then the estimated reduced aerodynamics model of the folding wing in the time domain can be represented as follows: " #
1 2 ^f s ðθÞ ¼ q A ^ 0 ðθ Þ þ b A ^ 2 ðθÞs2 zs þ q D ^ 1 ðθÞs þ b A ^ ðθ Þ ^ ðθÞ sI V R (28) E^ ðθÞszs 1 1 V b V2 Combining the structure model of Eq. (15) and the aerodynamics model of Eq. (28), the aeroelastic response of the folding wing can be simulated by time marching technique. Starting from an initial state, the motion responses of the reduced structural nodes and the folding angle at the current time step are obtained by solving Eq. (15). Then substituting these results into Eq. (28), the unsteady aerodynamic force at the reduced structural nodes is calculated. Substituting the aerodynamic force back into Eq. (15), the motion responses of the reduced structural nodes and the folding angle at the next time step can be obtained. By repeating these steps until sufficient response time history is achieved, the aeroelastic response of the folding wing can be obtained.
4. Numerical examples The geometries of the folding wing in the present study are shown in Fig. 4. The folding wing contains three substructures, i.e. the central wing, the inboard wing and the outboard wing. In the present study, the structure of the folding wing is represented by three uniform plates and the thicknesses of them are 5 cm, 2 cm and 2 cm. These substructures are made of aluminum, and are connected by rotation hinge springs, which are on the hinge lines as shown in Fig. 4 and marked by circles. The rotation stiffness of each spring is 800 kN m/rad.
Fig. 5. The comparison of natural frequencies at different folding angles.
Fig. 6. The locations of the reduced nodes of the aerodynamics model.
224
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
^ 0 at different folding angles. Fig. 7. The comparisons of some element values of A0 and A
^ 0 when θ¼82.5°, (a) matrix A0 and (b) matrix A ^ 0. Fig. 8. The comparison of element values of A0 and A
4.1. Examination of structure model and aerodynamics model The flexible multi-body dynamics model of the folding wing is built by the procedures proposed in Sections 2.1 and 2.2, and the natural frequencies of this model at different fixed folding angles are compared with those obtained from the finite element model, which is built by MSC.Nastran. The comparison of results in Fig. 5 shows a good agreement between the natural frequencies obtained by using these two approaches. By reducing the aerodynamics model with the spline interpolation technique, the unsteady AIC matrix is reduced into matrix Ass ðωÞ with the dimension of 24 24, and the locations of these 24 structural nodes are shown in Fig. 6. Then, the minimum state approximation is adopted to calculate the coefficient matrices at some specific fixed folding angles, and the Kriging agent model technique is used to build the aerodynamics model. A series of folding angles θ ¼ 01 151 301 451 601 751 901 1051 1201 are selected as the sample points to ^ 0 at different build the agent model. Fig. 7 shows the comparisons of some element values of A0 and its agent model A folding angles. Additionally, in order to further verify the accuracy of the agent model obtained by the Kriging technique, the matrices ^ 0 at the folding angle of 82:51, which is not included in the sample points, are compared graphically in Fig. 8. The A0 and A ^ 0 are shown in Fig. 9, and the errors are within 5%. percent errors of element values of A0 and A 4.2. Flutter analysis Based on the structure model and aerodynamics model obtained in Section 4.1, the aeroelastic responses of the folding wing at a specific fixed folding angle can be simulated under a series of air speeds. The flutter speed can be obtained at the air speed when the aeroelastic response is harmonic, and the flutter frequency can be calculated from this response. By these procedures, the flutter characteristics of the folding wing at different folding angles are obtained. These flutter
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
225
^ 0 when θ¼ 82.5°. Fig. 9. The percent errors of element values of A0 and A
Fig. 10. The comparison of the flutter characteristics obtained by different methods at different folding angles.
Fig. 11. The dynamic aeroelastic response of the folding wing with folding rate of 51/s.
characteristics are compared with those obtained by the g method using the finite element models of the folding wing structure and the Doublet Lattice based aerodynamics models corresponding to the configurations of different folding angles. The comparison between these results is shown in Fig. 10, and as expected, the results agree well. Fig. 10 also shows that the flutter speed decreases as the folding angle is increased from 01 to 601, whereas the flutter speed increases as the folding angle is increased from 601 to 1201. The lowest flutter speed is about 128 m/s when the folding angle is about 601. It also indicates that one simply needs to fly below the lowest flutter speed while morphing, but then there is some additional speed margin to the flight envelope based on the folding angle, because the flutter speed is shifting according to the folding angle. What if the flight speed is larger than the lowest flutter speed while morphing? A baseline speed of 147 m/s is chosen to examine the aeroelastic response of the folding wing during the morphing process. It is indicated by the solid line in Fig. 10 and this value is chosen because it is expected the wing will flutter during the morphing process as the flutter speed between 331 and 771 drops below this value.
226
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
Fig. 12. The dynamic aeroelastic response of the folding wing with folding rate of 101/s.
Fig. 13. The dynamic aeroelastic response of the folding wing with folding rate of 201/s.
In the simulation, the angle of attack of the folding wing is 01. The aerodynamic force during the morphing process consists of two parts: the unsteady aerodynamic force induced by the vibration of the folding wing, and the unsteady aerodynamic force induced by the rigid body motion of the folding wing. In order to examine the influence of the aerodynamic force induced by the rigid body motion on the aeroelastic characteristics, two cases are studied for the aeroelastic response during the morphing process: Case 1. Aeroelastic response of the folding wing at the speed of 147 m/s considering the aerodynamic force induced by the rigid body motion. Case 2. Aeroelastic response of the folding wing at the speed of 147 m/s without considering the aerodynamic force induced by the rigid body motion. The deflection responses at the trailing edge tip of the outboard wing during the folding processes with different folding rates are shown in Figs. 11–13. The dynamic aeroelastic stability obtained by the quasi-steady aeroelastic analysis is also shown in the figures, so that it can be compared directly with those obtained by the dynamic aeroelastic analysis during the morphing process. It can be seen from Fig. 10 that when the air speed is 147 m/s, the wing loses its dynamic aeroelastic stability when the folding angles are in the range of 331 to 771. Thus the dynamic aeroelastic stability of the folding wing at different folding angles can be divided into three regions: two stable regions corresponding to the folding angle smaller than 331 (S1) and larger than 771 (S2), and one unstable region (U) corresponding to the folding angle ranging from 331 to 771. At the same folding rate, the difference of the responses between Case 1 and Case 2 indicates that the aerodynamic force induced by the rigid body motion just changes the static equilibrium position of the responses and has no influence on the dynamic aeroelastic stability, i.e. the influence of this aerodynamic force on the aeroelastic characteristics is similar to the influence of the gust aerodynamic force. And it also shows that when the folding angle is increased from 01 to 1201, the shifting value of the static equilibrium position of the responses of Case 1 changes from negative to positive, and the shifting value is zero when the folding angle is 901. This phenomenon can be explained as follows. The shifting of the static equilibrium position of the response is caused by the aerodynamic force induced by the rigid body motion, i.e., the motion velocity and acceleration normal to the wing surface caused by the rigid body motion. As shown in Fig. 14, for the outboard wing, when the folding angle is increased from 01 to 1201, the direction of the motion velocity component normal to the wing surface is changed from upward to downward, and becomes zero when the folding angle is 901. And because the direction of aerodynamic force induced by the rigid body motion is opposite to the direction of
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
227
Fig. 14. The direction of the motion velocity component during the folding process.
Fig. 15. The dynamic aeroelastic response of the folding wing from 701 to 901 with folding rate of 0:51/s.
Fig. 16. The dynamic aeroelastic response of the folding wing with unfolding rate of 51/s.
the motion velocity component, thus, when the folding angle is increased from 01 to 901, the direction of the aerodynamic force induced by the rigid body motion is downward, and it becomes upward when the folding angle is increased from 901 to 1201, so that the shifting value of the static equilibrium position changes from negative to positive, and the shifting value is zero when the folding angle is 901. Figs. 11–13 also indicate that the dynamic aeroelastic stabilities obtained by the aeroelastic responses of Case 1 are different from that obtained by the quasi-steady aeroelastic analysis. During the folding process, the dynamic aeroelastic stability changes from stable state to unstable state when the folding angle is slightly larger than 331, and with the increasing of the folding rate, the critical folding angle at which the dynamic aeroelastic stability is lost becomes much larger than 331. The similar conclusion can be found in [21], where the aeroelastic characteristics of a variable-span wing during the morphing process were investigated. Besides, according to the aeroelastic stability obtained by the quasi-steady aeroelastic analysis, the amplitude of the response should increase when the folding angle is increased from 331 to 771, and decrease when the folding angle is increased from 771 to 1201. The responses in Figs. 11–13 indicate that when the folding angle is changed into a region nearby 771, a dip of the amplitude of the response appears, and the width of the region increases with the increasing of the folding rate. The two aforementioned phenomena are caused by the motions of the inboard and the outboard wings. According to Eq. _ and (15), the quadratic velocity vector is generated because of the motions of the substructures, and can be represented by X
228
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
Fig. 17. The dynamic aeroelastic response of the folding wing with unfolding rate of 101/s.
Fig. 18. The dynamic aeroelastic response of the folding wing with unfolding rate of 201/s.
Fig. 19. The dynamic aeroelastic response of the folding wing from 401 to 251 with unfolding rate of 0:51/s.
its coefficient matrices. If we treat the coefficient matrices as damping matrices, it means that the motions of the inboard and the outboard wings will induce additional damping into the aeroelastic system of the folding wing. The additional damping will influence the aeroelastic characteristics of the folding wing, especially when the folding angle is approaching to the specific folding angle at which the flutter speed equals to the air speed. Moreover, this influence increases with the increasing of the folding rate. From the conclusion in the above paragraph, it is definite that if the folding rate is very small, the dynamic aeroelastic stability of the folding wing must be similar to that obtained by the quasi-steady aeroelastic analysis. Considering the aerodynamic force induced by the rigid body motion, when the folding wing is folding from 701 to 901 with the folding rate of 0:51/s, the deflection response is shown in Fig. 15, and the dynamic aeroelastic stability of the folding wing agrees well with that obtained by the quasi-steady aeroelastic analysis. The unfolding processes with different unfolding rates are also simulated, and the deflection responses at the trailing edge tip of the outboard wing are shown in Figs. 16–18 , respectively. Both the deflection responses of Cases 1 and 2 are
W. Hu et al. / Journal of Sound and Vibration 365 (2016) 216–229
229
displayed in these figures, and the difference of the responses between Cases 1 and 2 show the same conclusion about the influence of the aerodynamic force induced by the rigid body motion on the aeroelastic characteristics of the folding wing. The aeroelastic responses of Case 1 as displayed in Figs. 16–18 show that during the unfolding process, the dynamic aeroelastic stability changes from unstable state to stable state when the folding angle is less than 331, and this angle becomes much less than 331 with the increasing of the unfolding rate. While as shown in Figs. 11–13, the folding angle at which the dynamic aeroelastic stability changes from stable state to unstable state is larger than 331, this phenomenon indicates that the folding and unfolding processes have opposite influences on the dynamic aeroelastic stability of the folding wing, and the influences become much more significant with the increasing of the folding and unfolding rates. The similar conclusion can also be found in [21]. The deflection response of the folding wing from 401 to 251 with the unfolding rate of 0:51/s is shown in Fig. 19, and the dynamic aeroelastic stability of the folding wing also agrees well with that obtained by the quasi-steady aeroelastic analysis.
5. Conclusions To study the aeroelastic characteristics of a folding wing during the morphing process, the structure model is built by incorporating the Craig–Bampton modal synthesis technique with the flexible multi-body dynamics approach. A method based on the Kriging agent model technique is proposed to build the aerodynamics model, which is proved to be accurate and can be used to calculate the aerodynamic force of the folding wing during the morphing process efficiently. The aeroelastic responses during the folding and unfolding processes are obtained by numerical simulation. The results show that the morphing process has an influence on the dynamic aeroelastic stability of the folding wing, which makes the dynamic aeroelastic stabilities of the folding wing different from that obtained by the quasi-steady aeroelastic analysis, and the influences of the folding and unfolding processes on the dynamic aeroelastic stability are also different. Furthermore, the influence of the morphing process on the dynamic aeroelastic stability becomes much more significant with the increasing of the morphing rate. Acknowledgments This work was founded by the National Natural Science Foundation of China (Grants no. 11472216), 111 Project of China (B07050).
References [1] M.H. Love, P.S. Zink, R.L. Stroud, et al., Demonstration of morphing technology through ground and wind tunnel tests, AIAA Paper 2007-1729, May 2007. [2] M. Snyder, G. Frank, B. Sanders, Aeroelastic analysis of a morphing Z-wing configuration, Aerospace Flutter and Dynamics Council Meeting, May 2004. [3] J. Dunn, L. Horta, T. Ivanco, et al., NASA contributions to DARPA MAS program, Aerospace Flutter and Dynamics Council Meeting, May 2004. [4] M. Snyder, B. Sanders, F. Eastep, et al., Sensitivity of flutter to fold orientation and spring stiffness of a simple folding wing, Proceedings of the International Forum on Aeroelasticity and Structural Dynamics 2005, June–July 2005. [5] D.H. Lee, T.A. Weisshaar, Aeroelastic studies on a folding wing configuration, AIAA Paper 2005-1996, April 2005. [6] Y. Zhao, H. Hu, Parameterized aeroelastic modeling and flutter analysis for a folding wing, Journal of Sound and Vibration 331 (2012) 308–324. [7] D.H. Lee, P.C. Chen, Nonlinear aeroelastic studies on a folding wing configuration with free-play hinge nonlinearity, AIAA Paper 2006-1734, May 2006. [8] D. Tang, E.H. Dowell, Theoretical and experimental aeroelastic study for folding wing structures, Journal of Aircraft 45 (2008) 1136–1147. [9] P.J. Attar, D. Tang, E.H. Dowell, Nonlinear aeroelastic study for folding wing structures, AIAA Journal 48 (2010) 2187–2195. [10] Y. Zhao, H. Hu, Prediction of transient responses of a folding wing during the morphing process, Aerospace Science and Technology 24 (2013) 89–94. [11] G.W. Reich, J.C. Bowman J., B. Sanders, et al., Development of an integrated aeroelastic multibody morphing simulation tool, AIAA Paper 2006-1892, May 2006. [12] J.C. Bowman, G.W. Reich, B. Sanders, et al., Simulation tool for analyzing complex shape-changing mechanisms in aircraft, AIAA Paper 2006-6727, August 2006. [13] J.N. Scarlett, R.A. Canfield, B. Sanders, Multibody dynamic aeroelastic simulation of a folding wing aircraft, AIAA Paper 2006-2135, May 2006. [14] R.R. Craig, Coupling of substructures for dynamic analyses, AIAA Journal 6 (1968) 1313–1319. [15] A.A. Shabana, Dynamics of Multibody System, Cambridge University Press, New York, 1998. [16] W.P. Rodden, J.P. Giesing, T.P. Kalman, New method for nonplanar configurations, AGARD Conference Proceedings, 1971. [17] R.L. Harder, R.N. Desmarais, Interpolation using surface splines, Journal of Aircraft 9 (1972) 189–191. [18] J. Duchon, Splines minimizing rotation-invariant semi-norms in sobolev spaces, Constructive Theory of Functions of Several Variables 571 (1977) 85–100. [19] M. Karpel, Procedures and models for aeroservoelastic analysis and design, Journal of Applied Mathematics and Mechanics 81 (2001) 579–592. [20] J. Sacks, W.J. Welch, T.J. Mitchell, et al., Design and analysis of computer experiments, Statistical Science 4 (1989) 409–435. [21] R. Huang, Z. Qiu, Transient aeroelastic responses and flutter analysis of a variable-span wing during the morphing process, Chinese Journal of Aeronautics 26 (6) (2013) 1430–1438.