Open access peer-reviewed chapter

Nonlinear Aeroelastic Response of Highly Flexible Flying Wing Due to Different Gust Loads

By Ehsan Izadpanahi and Pezhman Mardanpour

Submitted: November 28th 2017Reviewed: February 21st 2018Published: July 18th 2018

DOI: 10.5772/intechopen.75804

Downloaded: 348

Abstract

Nonlinear aeroelastic responses of a flying wing aircraft due to different gust profiles are investigated. Three different gust profiles are obtained considering light, moderate, and severe turbulence. A flying wing configuration is designed for the purpose of this investigation. The structural properties of the wings are obtained using VABS software, and then the flying wing is simulated with Nonlinear Aeroelastic Trim and Stability of HALE Aircraft (NATASHA) computer program. The results of time domain analysis are reported for the cases when engine is placed at the root of the wing and close to the area of maximum flutter speed. It has been found that the flying wing experiences limit cycle oscillation, when the engines are mounted at the root of the aircraft, for all three gust profiles. However, when the engines are placed at the area of maximum flutter speed, the oscillations die out. In addition, the real and imaginary part of eigenvalues and the unstable mode shape of the aircraft are reported.

Keywords

  • gust response
  • flying wing
  • nonlinear time domain analysis
  • flutter analysis
  • gust suppression

1. Introduction

Very flexible high-aspect-ratio wings are widely used in the design of high altitude long endurance (HALE) aircrafts. These wings due to their characteristics may subject to large deformation, which causes geometric nonlinearities. As a result, conducting the nonlinear aeroelastic analysis is necessary when it comes to the design of very flexible configurations [1, 2, 3]. In addition, time-dependent external excitation including gust [4, 5, 6, 7] and blast [8, 9, 10, 11, 12] can lead to instability even if the aircraft is flying below the stability boundary. Therefore, the determination of nonlinear aeroelastic responses to time-dependent excitation is a crucial topic for the design of very flexible flying wings.

Gust loads can result in large deformations in the case of a highly flexible aircraft. The flight dynamic characteristics and gust response of highly flexible aircraft were investigated by Patil and Taylor [13]. It was reported that the non-uniform gust creates higher responses in a case of high-aspect-ratio flying wing compared to uniform gusts. In addition, the nonlinear gust response of a highly flexible aircraft was reported by Patil [14], which he found that the time domain response matches with frequency domain response presented in the work by Patil and Taylor [13]. Ricciardi et al. [15] investigated the accuracy of the Pratt method for unconventional HALE aircraft. The Pratt method and transient method were used to analyze the gust response on the joined-wing and flying-wing model. It was found that Pratt method is only useful for the preliminary design of the joined-wing model. However, when it comes to the design of flying-wing model, the Pratt method is inadequate. Yi et al. [16] compared a theoretical and experimental approach of a flexible high-aspect-ratio wing exposed to a harmonic gust. It was found that a very flexible wing experiences different gust response characteristics under different load conditions and the responses are difficult to evaluate using linear analysis.

On the other hand, finding ways to suppress the responses of a highly flexible configuration due to time-dependent excitations is a challenging aspect of design. Tang et al. [17] conducted an experimental and theoretical study to investigate the effect of store span location and its pitch stiffness on the flutter velocity and LCO. A delta wing for the purpose of experimentation was chosen. In addition, the von-Karman plate theory, three-dimensional vortex lattice model, and slender body aerodynamic theory were used for modeling the wing structure and determining the aerodynamic loads, respectively. It was reported that the experimental investigation and theoretical studies were in good agreement, and they showed that the structural natural frequency of the wing/store declines as the store moves from the root to the tip of the wing. They concluded that mounting the store at the leading edge of the wing tip leads to a higher critical flutter velocity. Moreover, Mardanpour et al. [18] found that the maximum flutter speed happened for engine placement at 60% of span forward the reference line. It was reported that the body-freedom flutter mode was unaffected by the engine location except for cases in which the engine was mounted at the wing tip and near the reference line.

Fazelzadeh et al. [19] investigated the effects of a nonlinear active control system on the flutter vibration of a wing/store exposed to a random gust disturbance. It was found that the control system is effective in suppressing the flutter vibration. In addition, Mardanpour et al. [7, 20] reported that the gust response of a very flexible high-aspect-ratio wing can be suppressed by changing the location of the engine. It was found that placing the engine close to 75% of the span forward of the reference line increases the flutter speed and also leads to suppression of the LCO due to gust loads.

In this chapter, the effect of engine placement on nonlinear aeroelastic gust response of a flying wing aircraft is investigated using three gust profiles with different gust intensities. The gust profiles are obtained utilizing different magnitude of turbulence at 10,000 m of altitude [21]. A flying wing aircraft is simulated for this study. The wings are designed using the structural properties which were obtained utilizing VABS software for NACA0012 airfoil. The computer program Nonlinear Aeroelastic Trim and Stability of High Altitude Long Endurance Aircraft (NATASHA) [3, 22] is used to simulate the nonlinear behavior of the flying wing aircraft. NATASHA is a powerful tool for the simulation of nonlinear behavior of HALE aircraft. It uses the nonlinear composite beam theory [23] that accommodates the modeling of high-aspect-ratio wings and the aerodynamic theory of Peters et al. [24] to model the aerodynamic forces and the p method to evaluate the aeroelastic stability. NATASHA has been verified and validated against experimental and theoretical studies many times [25, 26]. The nonlinear responses of the aircraft are obtained for the cases when the engines are mounted at the root of wings and at the area of maximum flutter speed (i.e., 60% of span forward of reference line).

2. Theory

2.1. Nonlinear composite beam theory

The equations of motion, which are presented in Eq. (1), are based on force, moment, angular velocity, and velocity with nonlinearities of second order. These variables can be expressed in the bases of the deformed and undeformed frames, Bx1tand bx1, respectively, see Figure 1.

FB+K˜BFB+fB=ṖB+Ω˜BPBMB+K˜BMB+e˜1+γ˜FB+mB=ḢB+Ω˜BHB+V˜BPBE1

Figure 1.

Sketch of beam kinematics.

In this set of equations, FBand MBrepresent the column matrices of cross-sectional stress and moment resultant; VBand ΩBdefine column matrices of cross-sectional frame velocity and angular velocity; PBand HBindicate the column matrices of cross-sectional linear and angular momentum measures; K˜Bis Column matrix of deformed beam’s curvature and twist. All of the abovementioned variables measure in Bibasis. The structural and the inertial constitutive equations relate the stress resultants and moments to the generalized strains and velocities as follows:

γκ=RSSTTFBMBE2
PBHB=μΔμξ˜μξ˜IVBΩBE3

here, R, S, and Trepresent 3 ×3 partitions of the cross-sectional flexibility matrix; μis the mass per unit length; Δis the 3 ×3 identity matrix; Idefines the 3 ×3 cross-sectional inertia matrix; ξis 0ξ2ξ3Tin which ξ2and ξ3represent the position coordinates of the cross-sectional mass center with respect to the reference line. Finally, strain- and velocity-displacement equations are utilized to derive the intrinsic kinematical partial differential Equations [23].

VB+K˜BVB+e˜1+γ˜ΩB=γ̇ΩB+K˜BΩB=κ̇E4

In these equations, the tilde ˜represents the antisymmetric 3 ×3 matrix associated with the column matrix over which the tilde is placed, (̇)defines the partial derivative with respect to time, and is the partial derivative with respect to the axial coordinate, x1. More details about these equations can be found in Ref. [27]. In order to solve these first-order, partial differential equations, one may eliminate γand κusing Eq. (2) and PBand HBusing Eq. (3), and also 12 boundary conditions are required, in terms of force (FB), moment (MB), velocity (VB), and angular velocity (ΩB). Displacement and rotation variables do not appear in this formulation, and singularities due to finite rotations are avoided. The position and the orientation can be obtained as postprocessing operations by integrating

ri=Cibe1ri+ui'=CiBe1+γE5

and

Cbi=k˜CbiCBi=k˜+κ˜CBiE6

2.2. Finite state-induced model of Peters et al.

The aerodynamic model of Peters et al. [24] is utilized in this study. This finite state model is a state-space, thin-airfoil, inviscid, incompressible approximation of an infinite-state representation of the aerodynamic loads. By using known airfoil parameters, it can consider induced flow in the wake and apparent mass effects. In addition, it can accommodate large motion of the airfoil as well as deflection of a small trailing-edge flap. Available studies in literature [24, 25, 26] indicate that although this model cannot simulate the three-dimensional effects associated with the wing tip, it can accurately approximate the aerodynamic loads acting on high-aspect-ratio wings. The lift, drag, and pitching moment at the quarter-chord are given by

Laero=ρbcl0+clββVTVa2clαV̇a3b/2clαVa2Va3+λ0Ωa1b/2cd0VTVa3E7
Daero=ρbcl0+clββVTVa3+clαVa3+λ02cd0VTVa2E8
Maero=2ρbcm0+cmββVTcmαVTVa3bclα/8Va2Ωa1b2clαΩ̇a1/32+bclαV̇a3/8E9

Where,

VT=Va22+Va32.E10
sinα=Va3VTE11
αrot=Ωa1b/2VTE12

and βis the angle of flap deflection, Va2and Va3denote the measure numbers of Va. The effect of unsteady wake (induced flow) and apparent mass included as λ0and acceleration terms in the force and moment equation, which λ0can be calculated using the induced flow model of Peters et al. [24]:

Ainduced flowλ̇+VTbλ=V̇a3+b2Ω̇a1cinduced flowE13
λ0=12binduced flowTλE14

here, λdefines the column matrix of induced flow states, and Ainduced flow, cinduced flow, binduced flowrepresent constant matrices, which are derived in Ref. [24].

2.3. Gust airloads model

The gust airloads are taken into account separately from the aerodynamic forces of the flight dynamic velocities. The unsteady gust model measures the chordwise variation of the gust field on the deformed state of the wing. Here, an interpretation of the Peters and Johnson [28] theory that considers these effects is provided. The total induced flow is ωB, defining the vertical gust velocity in the deformed beam frame

L¯=ω0+12ω1+12ω̇0+12ω̇1bVTE15

here, L¯denotes the velocity-normalized lift coefficient presented by Peters and Johnson [28]; ωnis the coefficient of the nth Chebychev polynomial mode shape. ωBcan approximated as

ωB=0NωnTnE16

where Tnis the nth order Chebyshev polynomial. The gust force can be provided as

fgust=0ρbCV3+ω0L¯ρbCV2L¯E17

and the gust contribution to the induced flow can be presented as

λ0gust=ω̇0+12ω̇1E18

2.4. Aeroelastic system

By unifying the aerodynamic equations with the structural equations, the aeroelastic system is constructed

Aẋ+Bx=fcont+fgustE19

here, x, fcont, and fgustdefine the vectors of all of the aeroelastic variables, the flight controls, and gust loads, respectively. The resulting nonlinear ordinary differential equations are then linearized about a static equilibrium state, which is obtained by nonlinear algebraic equations. Utilizing the Newton-Raphson procedure, NATASHA solves these equations to obtain the steady-state trim solution [3]. The stability of the structure can be analyzed by linearizing this system of nonlinear aeroelastic equations about the resulting trim state, which leads to a standard eigenvalue problem. The linearized system is represented as

Ax̂̇+Bx̂=f̂cont+fgustE20

where (̂)is the perturbation about the steady-state values.

2.5. Transient gust response

The dynamic aeroelastic equations are solved in time to obtain the transient gust response. A central difference scheme in time-marching algorithm is used with a high-frequency damping. The linearized system in time can be written as follows:

1δtAx̂t+δtx̂t+12B1+ςx̂t+δt+1ςx̂t=f̂cont+fgustE21

here, δtand ςare the time step and the high-frequency-damping parameter, respectively. Utilizing ςapproximately equal to 0.01 provides a good time-marching algorithm, which the results are close to the central difference method.

The gust profiles are presented in Figure 2. These profiles presented in Figure 2 are generated by passing the Gaussian white noise through the Dryden spectrum model.

Figure 2.

Gust velocity profile versus time.

3. Case study

A very flexible high-aspect-ratio flying wing (see Figure 3) is designed in order to investigate the effects of different gust loads. The properties of the flying wing are presented in Table 1. The wings are aft swept 15 , and each wing has 20 elements. The fuselage is considered as a rigid body which contains four elements. The weight of each element of fuselage is five times of the weight of the elements of the wings. The aircraft has two engines with the mass of 10 kg.

Figure 3.

A schematic 3D view of a very flexible high-aspect-ratio wing.

PropertyValue
Span16
Number of elements20
Sweep angle15
R9.06×1090003.50×1087.22×101307.22×10131.18×106
S02.63×10127.57×10113.01×1012001.02×10600
T4.33×1060005.53×1062.42×101402.42×10148.43×108
I4.78×1010007.2×1031.04×101001.04×10104.71×101
ξ08.98×1044.76×107
Mass per unit length4.38
Chord, c1
Offset of aerodynamic center from reference line, e0.125
clα2π
clδ1
cd00.01
cm00.0
cmα0.08
Gravity, g9.8
Air Density, ρ0.4135

Table 1.

Properties of wing in SI units.

4. Results and discussion

In this section, two cases are considered. First, when the engine mounted at the root of the wing and the second case when the engines are located at 60% of the span forward of the reference line. For each case, the eigenvalues, the unstable mode shape of the aircraft, and the nonlinear time domain responses to the gust profiles are reported. The velocity results are normalized with the aircraft cruise speed of 50 m/s. The wing tip deflections also normalized with the length of the entire flying wing (i.e., 35.2 m), and the time is normalized with the period of oscillation of the flying wing at the flutter boundary when the engines are located at the root (i.e., 0.129 s).

4.1. Engine at the root

When the engines are located at the root of the flying wing, the wings experience a flutter at the speed of 48.9 m/s with a frequency of 7.7 rad/s. The real and imaginary parts of the eigenvalues are shown in Figure 4. In addition, the mode shape of the unstable mode is shown in Figure 5. The mode shape seems to contain first and second free-free bending mode.

Figure 4.

(a) Real part of eigenvalues and (b) imaginary part of eigenvalues.

Figure 5.

Unstable mode of the flying wing.

Figures 611 illustrate the results of time domain analysis when the engine is mounted at the root of the flying wing for different gust profiles in which Case 1, Case 2, and Case 3 indicate the results when the flying wing is exposed to light, moderate, and severe turbulence, respectively. It is found that the tip deflection increases in all directions when the gust load changes from light to severe turbulence. The same also happens for velocities. The velocity of the wing tip in different directions increases.

Figure 6.

Normalized wing tip position r1+u1l versus normalized time ttN. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 7.

Normalized velocity vector of wing tip V1VN versus normalized wing tip position r1+u1l. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 8.

Normalized wing tip position r2+u2l versus normalized time ttN. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 9.

Normalized velocity vector of wing tip V2VN versus normalized wing tip position r2+u2l. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 10.

Normalized wing tip position r3+u3l versus normalized time ttN. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 11.

Normalized velocity vector of wing tip V3VN versus normalized wing tip position r3+u3l. (a) Case 1, (b) Case 2, and (c) Case 3.

4.2. Engine at 60% of span forward of reference line

In another case, the engines are mounted at 60% of span forward of reference line. Mardanpour et al. [18] reported that this area coincides with the area of maximum flutter speed. It is found that the flying wing becomes unstable at the speed of 75.6 m/s. The real and imaginary parts of the eigenvalues are shown in Figure 12, and the mode shape of the unstable mode is displayed in Figure 13. Apparently, the mode shape only contains the first symmetric free-free bending mode.

Figure 12.

(a) Real part of eigenvalues and (b) imaginary part of eigenvalues.

Figure 13.

Unstable mode of the flying wing.

Figures 1416 show the results of time domain analysis when the engine is located at the area of maximum flutter speed (i.e., 60% of span forward of reference line). The results are reported for three different gust profiles. The results for this arrangement indicate that all the excitations from gust loads with different strength ranges from light to severe loads die out and the wing remains stable.

Figure 14.

Normalized wing tip position r1+u1l versus normalized time ttN. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 15.

Normalized wing tip position r2+u2l versus normalized time ttN. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 16.

Normalized wing tip position r3+u3l versus normalized time ttN. (a) Case 1, (b) Case 2, and (c) Case 3.

5. Conclusion

The nonlinear aeroelastic responses of a flying wing aircraft are investigated when the aircraft is exposed to different gust profiles with different intensities. The aircraft is designed with two aluminum wings with NACA 0012 airfoil. The properties of the wings are obtained using VABS software. The properties are then used in geometrically exact beam formulation, which is coupled with two-dimensional finite state aerodynamic model of Peters. The flutter characteristics for two configurations of the aircraft (i.e., engines at the root of the wings and engines at 60% of span forward of the reference line) as well as the eigenvalues and mode shape of the unstable modes for each configuration are studied. The flutter results are in agreement with the previous conclusion by Mardanpour et al. [18], which shows a higher flutter speed when the engines are mounted at 60% of span forward of reference line.

Three different gust profiles are then produced by passing white noise through Dryden gust model. The gust loads with light, moderate, and severe intensities are applied to the aircraft in time domain when the aircraft is cruising at 50 m/s. The results indicate that when the engines are mounted at the root of the wings, large oscillations exist, which their amplitude increases as the intensity of the gust loads increases. On the contrary, for all of the gust loads, when the engines are located at 60% of span forward of the reference line, the oscillations suppress. Previous study on gust alleviation by Mardanpour et al. [7, 20] for a cantilever wing also showed the suppression of gust responses when the engines are mounted at the area of maximum flutter speed.

Nomenclature

adeformed beam aerodynamic frame of reference
bundeformed beam cross-sectional frame of reference
Bdeformed beam cross-sectional frame of reference
biunit vectors in undeformed beam cross-sectional frame of reference (i=1,2,3)
Biunit vectors of deformed beam cross-sectional frame of reference (i=1,2,3)
cchord
cmβpitch moment coefficient w.r.t. flap deflection (β)
clαlift coefficient w.r.t. angle of attack (α)
clβlift coefficient w.r.t. flap deflection (β)
e1column matrix 100T
eoffset of aerodynamic center from the origin of frame of reference along b2
fcolumn matrix of distributed applied force measures in Bi basis
Fcolumn matrix of internal force measures in Bi basis
ggravitational vector in Bi basis
Hcolumn matrix of cross-sectional angular momentum measures in Bi basis
iinertial frame of reference
iiunit vectors for inertial frame of reference (i=1,2,3)
Icross-sectional inertia matrix
kcolumn matrix of undeformed beam initial curvature and twist measures in bi basis
Kcolumn matrix of deformed beam curvature and twist measures in Bi basis
lwing length
velocity-normalized lift coefficient
mcolumn matrix of distributed applied moment measures in Bi basis
Mcolumn matrix of internal moment measures in Bi basis
Pcolumn matrix of cross-sectional linear momentum measures in Bi basis
rcolumn matrix of position vector measures in bi basis
ucolumn matrix of displacement vector measures in bi basis
U∞free stream velocity
Vcolumn matrix of velocity measures in Bi basis
x1axial coordinate of beam
βtrailing edge flap angle
Δidentity matrix
γcolumn matrix of 1D-generalized force strain measures
κcolumn matrix of elastic twist and curvature measures (1D-generalized moment strain measures)
ηdimensionless position of the engine along the span
λcolumn matrix of induced flow states
Λsweep angle
μmass per unit length
ξcolumn matrix of center of mass offset from the frame of reference origin in bi basis
ψcolumn matrix of small incremental rotations
ωinduced flow velocity
Ωcolumn matrix of cross-sectional angular velocity measures in Bi basis
partial derivative of with respect to x1
(̇)partial derivative of with respect to time
(̂)nodal variable

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Ehsan Izadpanahi and Pezhman Mardanpour (July 18th 2018). Nonlinear Aeroelastic Response of Highly Flexible Flying Wing Due to Different Gust Loads, Nonlinear Systems - Modeling, Estimation, and Stability, Mahmut Reyhanoglu, IntechOpen, DOI: 10.5772/intechopen.75804. Available from:

chapter statistics

348total chapter downloads

1Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

A Reduced-Order Gauss-Newton Method for Nonlinear Problems Based on Compressed Sensing for PDE Applications

By Horacio Florez and Miguel Argáez

Related Book

First chapter

On Nonoscillatory Solutions of Two-Dimensional Nonlinear Dynamical Systems

By Elvan Akın and Özkan Öztürk

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us