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

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.

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, $B\left(x_1,t\right)$ and $b\left(x_1\right)$, respectively, see Figure 1.

$$\begin{array}{c}{F}_B^{\prime }+{{\displaystyle \tilde{K}}}_B{F}_B+f_B={{\displaystyle \dot{P}}}_B+{{\displaystyle \tilde{\mathrm{\Omega }}}}_B{P}_B\\ M_B^{\prime }+{{\displaystyle \tilde{K}}}_B{M}_B+\left({{\displaystyle \tilde{e}}}_1+{\displaystyle \tilde{\gamma }}\right)F_B+m_B={{\displaystyle \dot{H}}}_B+{{\displaystyle \tilde{\mathrm{\Omega }}}}_B{H}_B+{{\displaystyle \tilde{V}}}_B{P}_B\end{array}$$E1

Figure 1.

Sketch of beam kinematics.

In this set of equations, $F_B$ and $M_B$ represent the column matrices of cross-sectional stress and moment resultant; $V_B$ and ${\mathrm{\Omega }}_B$ define column matrices of cross-sectional frame velocity and angular velocity; $P_B$ and $H_B$ indicate the column matrices of cross-sectional linear and angular momentum measures; ${\tilde{K}}_B$ is Column matrix of deformed beam’s curvature and twist. All of the abovementioned variables measure in ${\mathbf{B}}_i$ basis. The structural and the inertial constitutive equations relate the stress resultants and moments to the generalized strains and velocities as follows:

$$\left\{\begin{array}{l}\gamma \\ \kappa \end{array}\right\}=\left[\begin{array}{ll}R& S\\ S^T& T\end{array}\right]\left\{\begin{array}{l}{F}_B\\ M_B\end{array}\right\}$$E2

$$\left\{\begin{array}{l}{P}_B\\ H_B\end{array}\right\}=\left[\begin{array}{ll}\mu \mathrm{\Delta }& -\mu {\displaystyle \tilde{\xi }}\\ \mu {\displaystyle \tilde{\xi }}& I\end{array}\right]\left\{\begin{array}{l}{V}_B\\ {\mathrm{\Omega }}_B\end{array}\right\}$$E3

here, $R$, $S$, and $T$ represent 3 $\times$ 3 partitions of the cross-sectional flexibility matrix; $\mu$ is the mass per unit length; $\mathrm{\Delta }$ is the 3 $\times$ 3 identity matrix; $I$ defines the 3 $\times$ 3 cross-sectional inertia matrix; $\xi$ is ${⌊0{\xi }_2{\xi }_3⌋}^T$ in which ${\xi }_2$ and ${\xi }_3$ represent 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].

$$\begin{array}{c}{V}_B^{\prime }+{{\displaystyle \tilde{K}}}_B{V}_B+\left({{\displaystyle \tilde{e}}}_1+{\displaystyle \tilde{\gamma }}\right){\mathrm{\Omega }}_B={\displaystyle \dot{\gamma }}\\ {\mathrm{\Omega }}_B^{\prime }+{{\displaystyle \tilde{K}}}_B{\mathrm{\Omega }}_B=\dot{\kappa }\end{array}$$E4

In these equations, the tilde $\widetilde{\left(\right)}$ represents the antisymmetric 3 $\times$ 3 matrix associated with the column matrix over which the tilde is placed, $\dot{(})$ defines the partial derivative with respect to time, and ${\left(\right)}^{\prime }$ is the partial derivative with respect to the axial coordinate, $x_1$. More details about these equations can be found in Ref. [27]. In order to solve these first-order, partial differential equations, one may eliminate $\gamma$ and $\kappa$ using Eq. (2) and $P_B$ and $H_B$ using Eq. (3), and also 12 boundary conditions are required, in terms of force ($F_B$), moment ($M_B$), velocity ($V_B$), and angular velocity (${\mathrm{\Omega }}_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

$$\begin{array}{c}{r}_i^{\prime }=C^{\mathit{ib}}{e}_1\\ r_i+{u_i}^{\text{'}}=C^{\mathit{iB}}\left(e_1+\gamma \right)\end{array}$$E5

and

$$\begin{array}{c}{\left(C^{\mathit{bi}}\right)}^{\prime }=-{\displaystyle \tilde{k}}{C}^{\mathit{bi}}\\ {\left(C^{\mathit{Bi}}\right)}^{\prime }=-\left(\tilde{k}+\tilde{\kappa }\right)C^{\mathit{Bi}}\end{array}$$E6

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

$$L_{\text{aero}}=\mathit{\rho b}\left[\left(c_{l_0}+c_{l_{\beta }}\beta \right)V_T{V}_{a_2}-c_{l_{\alpha }}{\dot{V}}_{a_3}b/2-c_{l_{\alpha }}{V}_{a_2}\left(V_{a_3}+{\lambda }_0-{\mathrm{\Omega }}_{a_1}b/2\right)-c_{d_0}{V}_T{V}_{a_3}\right]$$E7

$$\begin{array}{l}{D}_{\text{aero}}=\mathit{\rho b}\left[-\left(c_{l_0}+c_{l_{\beta }}\beta \right)V_T{V}_{a_3}+c_{l_{\alpha }}{\left(V_{a_3}+{\lambda }_0\right)}^2-c_{d_0}{V}_T{V}_{a_2}\right]\end{array}$$E8

$$\begin{array}{l}{M}_{\text{aero}}=2\mathit{\rho b}\left[\left(c_{m_0}+c_{m_{\beta }}\beta \right)V_T-c_{m_{\alpha }}{V}_T{V}_{a_3}-{\mathit{bc}}_{l_{\alpha }}/8V_{a_2}{\mathrm{\Omega }}_{a_1}-b^2{c}_{l_{\alpha }}{\dot{\mathrm{\Omega }}}_{a_1}/32+{\mathit{bc}}_{l_{\alpha }}{\dot{V}}_{a_3}/8\right]\end{array}$$E9

Where,

$$V_T=\sqrt{V_{a_2}^2+V_{a_3}^2}.$$E10

$$\text{sin}\alpha =\frac{-V_{a_3}}{V_T}$$E11

$${\alpha }_{\mathrm{rot}}=\frac{{\mathrm{\Omega }}_{a_1}b/2}{V_T}$$E12

and $\beta$ is the angle of flap deflection, $V_{a_2}$ and $V_{a_3}$ denote the measure numbers of $V_a$. The effect of unsteady wake (induced flow) and apparent mass included as ${\lambda }_0$ and acceleration terms in the force and moment equation, which ${\lambda }_0$ can be calculated using the induced flow model of Peters et al. [24]:

$$\begin{array}{l}\left[A_{\text{induced} \text{flow}}\right]\left\{\dot{\lambda }\right\}+\left(\frac{V_T}{b}\right)\left\{\mathrm{\lambda }\right\}=\left(-{\dot{V}}_{a_3}+\frac{b}{2}{\dot{\mathrm{\Omega }}}_{a_1}\right)\left\{c_{\text{induced} \text{flow}}\right\}\end{array}$$E13

$$\begin{array}{l}{\lambda }_0=\frac{1}{2}{\left\{b_{\text{induced} \text{flow}}\right\}}^T\left\{\lambda \right\}\end{array}$$E14

here, $\lambda$ defines the column matrix of induced flow states, and $\left[A_{\text{induced} \text{flow}}\right]$, $\left\{c_{\text{induced} \text{flow}}\right\}$, $\left\{b_{\text{induced} \text{flow}}\right\}$ represent 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 ${\omega }^B$, defining the vertical gust velocity in the deformed beam frame

$$\begin{array}{l}\overline{L}={\omega }_0+\frac{1}{2}{\omega }_1+\frac{1}{2}\left({\dot{\omega }}_0+\frac{1}{2}{\dot{\omega }}_1\right)\frac{b}{V_T}\end{array}$$E15

here, $\overline{L}$ denotes the velocity-normalized lift coefficient presented by Peters and Johnson [28]; ${\omega }_n$ is the coefficient of the nth Chebychev polynomial mode shape. ${\omega }^B$ can approximated as

$$\begin{array}{l}{\omega }^B=\sum _0^N{\omega }_n{T}_n\end{array}$$E16

where $T_n$ is the nth order Chebyshev polynomial. The gust force can be provided as

$$f_{\mathit{\text{gust}}}=\left\{\begin{array}{l}0\\ -\rho {\mathit{bC}}_{\mathit{l\alpha }}\left(V_3+{\omega }_0\right)\overline{L}\\ \rho {\mathit{bC}}_{\mathit{l\alpha }}{V}_2\overline{L}\end{array}\right\}$$E17

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

$$\begin{array}{l}{\lambda }_{0_{\mathit{\text{gust}}}}={\dot{\omega }}_0+\frac{1}{2}{\dot{\omega }}_1\end{array}$$E18

2.4. Aeroelastic system

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

$$\begin{array}{l}\left[A\right]\left\{\dot{x}\right\}+\left\{B\left(x\right)\right\}=\left\{f_{\mathit{\text{cont}}}\right\}+\left\{f_{\mathit{\text{gust}}}\right\}\end{array}$$E19

here, $\left\{x\right\}$, $\left\{f_{\mathit{\text{cont}}}\right\}$, and $\left\{f_{\mathit{\text{gust}}}\right\}$ define 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

$$\begin{array}{l}\left[A\right]\left\{\dot{\widehat{x}}\right\}+\left[B\right]\left\{\widehat{x}\right\}=\left\{{\widehat{f}}_{\mathit{\text{cont}}}\right\}+\left\{f_{\mathit{\text{gust}}}\right\}\end{array}$$E20

where $\widehat{(})$ 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:

$$\begin{array}{l}\frac{1}{\mathit{\delta t}}\left[A\right]\left(\left\{{\widehat{x}}^{t+\mathit{\delta t}}\right\}-\left\{{\widehat{x}}^t\right\}\right)+\frac{1}{2}\left[B\right]\left(\left(1+\varsigma \right)\left\{{\widehat{x}}^{t+\mathit{\delta t}}\right\}+\left(1-\varsigma \right)\left\{{\widehat{x}}^t\right\}\right)=\left\{{\widehat{f}}_{\mathit{\text{cont}}}\right\}+\left\{f_{\mathit{\text{gust}}}\right\}\end{array}$$E21

here, $\mathit{\delta t}$ and $\varsigma$ are the time step and the high-frequency-damping parameter, respectively. Utilizing $\varsigma$ 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 ${}^{\circ }$, 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.

Property	Value
Span	16
Number of elements	20
Sweep angle	$15$
R	$\left[\begin{array}{lll}9.06\times {10}^{-9}& 0& 0\\ 0& 3.50\times {10}^{-8}& 7.22\times {10}^{-13}\\ 0& 7.22\times {10}^{-13}& 1.18\times {10}^{-6}\end{array}\right]$
S	$\left[\begin{array}{lll}0& 2.63\times {10}^{-12}& 7.57\times {10}^{-11}\\ -3.01\times {10}^{-12}& 0& 0\\ -1.02\times {10}^{-6}& 0& 0\end{array}\right]$
T	$\left[\begin{array}{lll}4.33\times {10}^{-6}& 0& 0\\ 0& 5.53\times {10}^{-6}& 2.42\times {10}^{-14}\\ 0& 2.42\times {10}^{-14}& 8.43\times {10}^{-8}\end{array}\right]$
I	$\left[\begin{array}{lll}4.78\times {10}^{-1}& 0& 0\\ 0& 7.2\times {10}^{-3}& -1.04\times {10}^{-10}\\ 0& -1.04\times {10}^{-10}& 4.71\times {10}^{-1}\end{array}\right]$
$\xi$	$\left[\begin{array}{l}0\\ 8.98\times {10}^{-4}\\ -4.76\times {10}^{-7}\end{array}\right]$
Mass per unit length	4.38
Chord, c	1
Offset of aerodynamic center from reference line, e	0.125
$c_{l_{\alpha }}$	$2\pi$
$c_{l_{\delta }}$	$1$
$c_{d_0}$	$0.01$
$c_{m_0}$	$0.0$
$c_{m_{\alpha }}$	$-0.08$
Gravity, $g$	$9.8$
Air Density, $\rho$	$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 6–11 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 $\frac{r_1+u_1}{l}$ versus normalized time $\frac{t}{t_N}$. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 7.

Normalized velocity vector of wing tip $\frac{V_1}{V_N}$ versus normalized wing tip position $\frac{r_1+u_1}{l}$. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 8.

Normalized wing tip position $\frac{r_2+u_2}{l}$ versus normalized time $\frac{t}{t_N}$. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 9.

Normalized velocity vector of wing tip $\frac{V_2}{V_N}$ versus normalized wing tip position $\frac{r_2+u_2}{l}$. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 10.

Normalized wing tip position $\frac{r_3+u_3}{l}$ versus normalized time $\frac{t}{t_N}$. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 11.

Normalized velocity vector of wing tip $\frac{V_3}{V_N}$ versus normalized wing tip position $\frac{r_3+u_3}{l}$. (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 14–16 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 $\frac{r_1+u_1}{l}$ versus normalized time $\frac{t}{t_N}$. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 15.

Normalized wing tip position $\frac{r_2+u_2}{l}$ versus normalized time $\frac{t}{t_N}$. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 16.

Normalized wing tip position $\frac{r_3+u_3}{l}$ versus normalized time $\frac{t}{t_N}$. (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

a	deformed beam aerodynamic frame of reference
b	undeformed beam cross-sectional frame of reference
B	deformed beam cross-sectional frame of reference
bi	unit vectors in undeformed beam cross-sectional frame of reference (i=1,2,3)
Bi	unit vectors of deformed beam cross-sectional frame of reference (i=1,2,3)
c	chord
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 (β)
e1	column matrix 100T
e	offset of aerodynamic center from the origin of frame of reference along b2
f	column matrix of distributed applied force measures in Bi basis
F	column matrix of internal force measures in Bi basis
g	gravitational vector in Bi basis
H	column matrix of cross-sectional angular momentum measures in Bi basis
i	inertial frame of reference
ii	unit vectors for inertial frame of reference (i=1,2,3)
I	cross-sectional inertia matrix
k	column matrix of undeformed beam initial curvature and twist measures in bi basis
K	column matrix of deformed beam curvature and twist measures in Bi basis
l	wing length
L¯	velocity-normalized lift coefficient
m	column matrix of distributed applied moment measures in Bi basis
M	column matrix of internal moment measures in Bi basis
P	column matrix of cross-sectional linear momentum measures in Bi basis
r	column matrix of position vector measures in bi basis
u	column matrix of displacement vector measures in bi basis
U∞	free stream velocity
V	column matrix of velocity measures in Bi basis
x1	axial 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