Research on Aeroelasticity Phenomenon in Aeronautical Engineering

Abstract

Aeroelasticity phenomena arise when structural deformations induce changes on aerodynamic forces due to airplane structures that are not completely rigid. The additional aerodynamic forces cause an increase in the structural deformations, which leads to greater aerodynamic forces in a feedback process. These interactions may become smaller until reaching a condition of equilibrium or may diverge catastrophically if resonance occurs. Flutter is an instability aeroelasticity phenomenon which is the most difficult to predict. In this chapter, a numerical method and an experimental method were realized to predict aeroelastic response and characteristic parameters of a wing structure. The numerical method was firstly developed based on the interaction between computational fluid dynamic and computational structural dynamic methods using a coupling system, fluid–solid interaction (FSI), in the ANSYS software. Then, an experiment was set up in suitable conditions to study aeroelasticity characteristics with the goal of comparing the numerical results with the experimental results on the same wing structure at low speed. After that, a developed code based on immersed boundary method (IBM) was realized to predict aeroelasticity response and characteristic parameters of the wing structure. AGARD 445.6 wing model was chosen for this developed procedure at high speed. Obtained results were compared to other numerical and experimental results.

1. Introduction

Flutter is defined as the dynamic instability of aeroelasticity. Flutter is one of the most dangerous aeroelasticity phenomena as it could lead to a destroyed structure. The reason is the unsteady aerodynamic forces generated from elastic deformations of the structure that are usually involved with complicated phenomena such as the shock wave/boundary layer interaction, flow separation, nonlinear limited cycle oscillation, and more. Flutter is determined as a critical issue determining the reliability of the airplane wings or aircraft engine turbo-machine blades. Therefore, in the early phase of the structural design of the air vehicle, aircraft engine turbo-machinery, flutter problems should be calculate and predicted. However, accurate prediction of the flutter is very challenging due to the perplexing physical phenomena and the required large amount of computation [1, 2, 3, 4].

Coupled aeroelastic solution procedures use strongly coupled algorithms which contained sufficient interaction between computational fluid dynamics (CFD) and computational structural dynamics (CSD). As computer technology progresses, higher-order methods of CFD based on the Euler and the Navier-Stokes equations become more attractive due to the ability of the model and its more accurately transonic, nonlinear, and viscous effects. CFD has also advanced from two-dimensional problems to fully three-dimensional problems with or without coupled solution of the structural equations (CSD). The flow solvers used in aeroelastic analysis include 3D Euler and Navier-Stokes solvers which assumed inviscid flow.

The dynamic response of flutter characteristics of the first AGARD 445.6 wing standard aeroelastic configuration was studied using an unsteady Navier-Stokes algorithm in order to investigate a previously noted discrepancy between Euler flutter characteristics and the experimental data [5]. The 3D implicit upwind Euler/Navier-Stokes code (CFL3D Version 2.1) was previously modified for the time-marching aeroelastic analysis of wings using the unsteady Euler equations. A linear stability analysis and a time-marching aeroelastic analysis were used to determine the flutter characteristics of the isolated 45° swept-back wing. The flutter characteristics of the wing were determined using traditional V-g analysis. This stability analysis was determined at free-stream Mach numbers of 0.96 and 1.141 using the generalized aerodynamic forces calculated by solving the Euler equations and the Navier-Stokes equations.

Computational flutter required a fluid-structure interface as a common boundary to exchange the aerodynamic loads and structural displacements at the wing surface. But, the aerodynamic and structural grids were not coincident due to different systems used (fluent solver for aerodynamic grids and mechanical ADPL solver for structural grids). Therefore, the system coupling tool in ANSYS Workbench was used to transfer the aerodynamic pressure loads from the CFD grid points to the CSD grid points and vice versa, which ensured a conservative transfer of energy between the two systems [6].

Time-accurate aeroelastic simulations were carried out using the modal coupled aeroelastic implementation for a standard experimental test case: the AGARD 445.6 aeroelastic wind tunnel model in the subsonic and transonic regions [7]. A numerical methodology coupling Navier-Stokes equations and structural modal equations for predicting 3D transonic wing flutter was described in [8]. A modal approach is used for the structural response. The results indicate that the first five modes are sufficient to accurately model the wing structure response. In Ref. [9], an unsteady Reynolds-averaged Navier-Stokes (RANS) model was coupled with normal modes of structure to predict the flutter boundary for the AGARD 445.6 wing. A new integrated CFD-CSD simulation for flutter calculations based on a parallel, multiblock, multigrid flow solver for the Euler/Navier-Stokes equations using the ANSYS software was found in [10]. Computations were performed for a three-dimensional test case of AGARD 445.6 wing to validate and establish the usefulness of the simulation. Immersed boundary method solved Navier-Stokes equations for flow in couple with the Newton equation for structure movement under the effect of friction force exerted on the structure surface to carry out fluid–structure interaction (FSI). However, computational grids needed to be re-meshed in each time step due to changes of the structure position in time. To overcome this obstacle, immersed boundary method and finite volume methods were both invoked in solving the interaction between fluid flow and moving structure [11]. This research estimated the numerical and experimental results on the wing structure at a low speed with four different wing models such as two rectangular and two trapezoid 3D-shape wings, each 3D-shape wing had symmetric and asymmetric airfoil, respectively.

The subsonic aeroelastic stability of a two-dimensional panel resting on a continuous elastic foundation was investigated in Ref. [12]. Tests were conducted experimentally on a 104 × 24 × 0.018 in. rectangular aluminum panel in a low-speed wind tunnel. Comparison of experiment and theory showed a good agreement in flutter speed and wavelength but poor agreement in wave speed and frequency at flutter. This discrepancy was attributed to the limitations in the test setup as well as to the general difficulty of predicting the wave speed and frequency as accurately as the flutter speed. Reference [13] tested the first AGARD standard aeroelastic configuration for dynamic response, 445.6 wing, in the 16 foot transonic dynamics tunnel at the National Aeronautics and Space Administration (NASA) Langley Research Center. Several models of the wing were tested in the transonic dynamic tunnel including full-span and semispan models over a range of Mach number from 0.338 to 1.141. The NASA conducted experiments in wind tunnel to estimate the aeroelastic characteristics of new and advanced flight vehicles, including fixed-wing, rotary-wing, and space-launch configurations. Reviews and assessments were made regarding available facilities, measurement techniques, and other means and devices useful in testing. The needs and requirements for advances and improvements in testing capabilities for future experimental research and development programs are described [14].

In [15], aeroelastic concepts for increased aircraft performance were mentioned. Active aeroelastic concepts as well as robust analysis concepts aiming at efficient analysis were carried out using numerical models with uncertain or varying model parameters. A high aspect ratio wing in wind tunnel testing conditions was considered for exploitation of fluid-structure interaction of active aeroelastic structures. The structural flexibility was exploited by using multiple control surfaces such that the deformed wing shape gives minimum drag for different flight conditions. Two different drag minimization methods were carried out: one was to reduce induced drag based on numerical optimization techniques, and another was to reduce measured total drag using real-time optimization in the wind tunnel experiment.

An approach for the prediction of dynamic modal transient response and flutter characteristics of structures with unknown system parameters, such as stiffness and mass, using experimental modal parameters was remarked in [16]. A finite element model was created by using the actual material properties of the structure to study the correlation of the results. The computed transient responses and flutter velocities by the proposed method using experimental modal parameters observed that material properties were not a prerequisite.

In Ref. [17], transonic flutter characteristics of AGARD 445.6 wing between the numerical method [5] and the experimental data of NASA’s experiment were compared [13]. The comparison provided the basis for developing the numerical setup and the experimental setup to check out subsonic flutter characteristics of some simple wing structures (e.g., a thin plate). Aeroelasticity on airplane wing, which had supercritical airfoil, was carried out in [18]. The model wings were made from different materials and dimensions. Hence, varied wing structures were created to accomplish a comprehensive analysis, and the flutter velocity was also restricted to appropriate values within the working range of experiment devices. At the same time, infinite element method with the help of the ANSYS software was also conducted to simulate the phenomena on the same model wings as a verification for the precision of the experimental models.

2. Two-way FSI method

2.1 FSI method

The generalized equations of aeroelasticity motion:

E1

E2

where:

w: structural displacement at any time instant and position.

q: generalized displacement vector.

[M]: generalized mass matrice.

[C]: damping matrice.

[K]: stiffness matrices.

ϕ: normal modes of the structure.

N: total number of modes of the structure.

F: generalized force vector, which is responsible for linking the unsteady aerodynamics and inertial loads with the structural dynamics.

Eq. (1) shows that there are distinct terms representing the structure, aerodynamic, and dynamic disciplines. This equation was solved numerically by integrated CFD-CSD tool on the ANSYS software. This method was also called the two-way fluid-solid interaction (FSI). Aerodynamic loads were first calculated by CFD solver. Then, these loads were used to calculate the structural response of the wing structure through the fluid-solid interface. By using CSD solver, the structural deflection was estimated, and the mesh in each time step was deformed. The simulation of the two-way FSI is presented in Figure 1 [10].

Figure 1.

Two-way FSI algorithm.

Following [8], the first four modes of vibration were sufficient to accurately model the wing structure response. So, in order to determine the time step of the unsteady problem, a modal method was applied to estimate first the natural frequency of the first four modes and then calculate the time step using the following formula:

E3

where:

Δt: time step

f: natural frequency

The flutter velocity was estimated from the vibration of the wing tip position, the most dangerous position of the wing [5, 8, 9, 13]. From the variation of this position, the damping coefficient except the influence of structural damping was measured. It means that there was an effect of aerodynamic damping coefficient on the vibration of the wing structure. The aerodynamic damping coefficient was calculated as follows:

E4

where:

X_i: i^th peak of vibration.

n: number of periods.

ζ: aerodynamic damping coefficient.

The damping coefficient increased/decreased when the air velocity decreased/increased. The type of vibration was distinguished from the value of damping coefficient:

• ζ > 0: The vibration was convergent. The wing structure was stable.

• ζ < 0: The vibration was divergent. The wing structure was unstable.

• ζ = 0: The vibration was harmonic oscillation. The wing was in critical state. The air velocity was in the flutter velocity.

2.2 Wing model

AGARD 445.6 wing with aspect ratio of 1.65, taper ratio of 0.66, and 45^o sweep angle at quarter chord line was studied as seen in Figure 2. The cross section of the wing was NACA65A004 airfoil in the stream-wise direction. This NACA65A004 airfoil was a symmetric airfoil with a maximum thickness of 4% of the local chord. The dimensions of the wing were root chord of 0.558 m, tip chord of 0.368 m, and semispan of 0.762 m. The wing model used in aeroelastic experiments [13] was constructed by laminated mahogany which was modeled as an orthotropic material with different material properties in different directions. The properties of the laminated mahogany are given in Table 1. The modal analysis was performed using mechanical APDLs solver to evaluate the accuracy of the constructed model.

Figure 2.

Semispan AGARD 445.6 wing model.

Material property	E11	E22	E33	G	υ	ρ
Value	3.151	0.416	0.416	0.439	0.310	397.5
Unit	GPa	GPa	GPa	GPa	N/m^2	kg/m^3

Table 1.

Mechanical properties for the weakened AGARD 445.6 wing.

The AGARD 445.6 wing was modeled at a zero attack angle and at altitudes of 9.65 and 14 km, the same conditions of experimental study in [13]. The wing was meshed in 9257 nodes and 1350 elements (Figure 3a). The fluid domain in CFD problem was meshed in 67,949 nodes and 279,535 elements (Figure 3b).

Figure 3.

Computational grids. (a) AGARD 445.6 wing grids. (b) Fluid computational grids.

The free-stream air velocity was from 0.29 to 0.59 M at an altitude of 9.65 km and from 0.47 M to 0.73 M at an altitude of 14 km.

2.3 Results

2.3.1 Modal analysis

The mode shapes are obtained from the finite element analysis of the modeled wing. The deflection contours between the modal analysis and experiment [13] were compared as shown in Figure 4. The natural frequencies between the developed solution, experiment [13], and other researches were also compared as shown in Table 2. It could be concluded that the obtained results were in good agreement with the experimental results in [4, 6, 13] within an error relative of 8%. The frequency of the first mode was around 9.96 Hz. Following Eq. (3), the time step of the unsteady problem was estimated about 0.005 s.

Figure 4.

Comparison of the mode shapes. (a) Mode 1. (b) Mode 2. (c) Mode 3. (d) Mode 4.

	Mode 1	Mode 2	Mode 3	Mode 4
Modal analysis	9.96	41.5	52.95	99.99
Kolonay [4]	9.63	37.12	50.50	89.94
Erkut and Ali [6]	10.85	44.57	56.88	109.10
Yates [13]	9.6	38.1	50.7	98.5

Table 2.

Natural frequencies of modal analysis.

2.3.2 Damping analysis

At altitude 9.65 km, the Mach number of air velocity was varied from 0.29 to 0.59 (M = 0.29; 0.35; 0.41; 0.47; 0.53; 0.59). From Eq. (4), the aerodynamic damping coefficient was calculated and presented in Figure 5a. The zero damping coefficients were interpolated at Mach number 0.46. While the experimental zero damping was 0.499 in [13], the difference between simulation results and results of [13] was about 8%.

Figure 5.

Aerodynamic damping coefficients. (a) At altitude 9.65 km. (b) At altitude 14 km.

At altitude 14 km, the Mach numbers of air velocity were varied from 0.47 to 0.73 (M = 0.47; 0.53; 0.59; 0.65; 0.67; 0.73). From Eq. (4), the aerodynamic damping coefficient was calculated and presented in Figure 5b. The zero damping coefficients were interpolated at Mach number 0.58. While the experimental zero damping was 0.678 in [13], the difference within simulation results was about 14%.

In both the two considered altitudes (9.65 and 14 km), the numerical results agreed well with the experimental results of [13] with a relative error about 14%. This difference would be from the computation such as the quality of mesh and order of model in CFD and CSD.

For more details of the stability of the wing structure, the vibrated value of the wing tip position at three Mach numbers were plotted as shown in Figure 5.

At a Mach number smaller than the flutter value, the damping coefficient was positive, and the wing was stable (Figure 6a).

Figure 6.

Vibration of wing at altitude H = 9.65 km. (a) M = 0.29 – ζ = 0.0071. (b) M = 0.47 – ζ = −0.000125. (c) M = 0.53 – ζ = −0.006.

At a Mach number near the flutter value, the damping coefficient was zero, and the vibration was harmonic oscillation (Figure 6b).

At a Mach number greater than the flutter value, the damping coefficient was negative, and the vibration was divergent (Figure 6c). This divergent vibration would create the damage of the wing such as loss of control for the flap, aileron, fracture of wing, etc.

2.3.3 Simulation analysis

Dynamic aeroelastic analysis was a problem related to fluid-structure interaction over a period of time. Therefore, the quality of aerodynamic grid and the time step strongly influenced the results of aeroelastic analysis. These parameters were also two of the most important problems in the dynamic aeroelastic analysis.

In order to evaluate the quality of aerodynamic grid, the coefficient of pressure of AGARD 445.6 wing was first estimated at 26% semispan and at 75.5% semispan and then was compared with Ref. [6] as shown in Figure 7. The presented results were in good agreement with the results in Ref. [6]. It could be concluded that these simulation settings were appropriate for solving the transonic flow.

Figure 7.

Comparison of Cp distribution at M = 1.141, α = 0°. (a) 26 % semispan. (b) 75.5 % semispan.

To evaluate the time step size, three different time sizes were examined such as 0.001, 0.002, and 0.005 s. As it could be seen in Figure 8, the displacement of the wing tip was reduced with the reduction of time step size up to 0.002 s until the aeroelastic simulation did not change [6]. Therefore, the value 0.002 s of time step size in the numerical solution was chosen for both aerodynamic and structural analysis.

Figure 8.

Wing-tip oscillation depends on the selected time step.

The limit of flutter was identified by using damping estimations for a large test point at each Mach number. At M = 0.499, the oscillation of the displacement of the wing tip was harmonic (Figure 9), and it was considered as a flutter point. At this limit of flutter, the air speed was calculated as 174.26 m/s, and the density of air was calculated as 0.432 kg/m³. These values were very close to the experimental values: 172.46 m/s for flutter speed and 0.428 kg/m³ for density of air (Table 3). This remark illustrated that the developed solution could be used to specify the transonic flutter characteristics with errors less than 10%.

Figure 9.

Wing-tip oscillation of flutter point at M = 0.499.

	Velocity (m/s)	Density (kg/m^3)	Frequency at flutter (Hz)
Baskut [6]	171.84	0.3987	21.67
Yates [13]	172.5	0.42770	20.39
Simulation	174.257	0.43164	18.87
Relative error with [13]	1%	1%	7%

Table 3.

Flutter characteristics at M = 0.499.

3. Experimental method

3.1 Experimental setup

The test model was set in AF6116 (M = 0.1) subsonic wind tunnel located at the Hanoi University of Science and Technology, which was of a blowdown type with a closed test section (0.4 × 0.5 × 1.0 m³). The wind speed could be arbitrarily varied up to 30 m/s, where the Reynolds number based on wing root chord was 10⁶, which was driven by an 8 kW electric motor.

Flutter characteristics were determined with the help of the frequency meter and load cell, which allowed to specify the flutter frequency and root wing force, respectively. The oscillated frequency was measured by the DT-2234C+ frequency meter. The signal of measured frequency was averaged by five measurements. The force applied to the wing was measured by load cell system. In the experimental aeroelastic analysis, the flutter frequency and the flutter amplitude were measured at different velocities ranging from 10 to 30 m/s using an oscillator generator system.

3.2 Wing model

Two wing models with the parameter and dimension are shown in Figure 10 and Table 4. The non-structure wing had only balsa wood, while the structure wing had balsa wood for skin and carbon rod and hard wood for the inner parts.

Figure 10.

Experimental models. (a) Wing model. (b) Wing with support. (c). Broken wing.

	Wing 1 (non-structure)	Wing 2 (with structure)
Chord length	300 mm	300 mm
Root chord length	500 mm	500 mm
Tip chord length	100 mm	100 mm
Profile	NACA65A004	NACA65A004
Taper ratio	0.5	0.5
Material	Balsa wood	Balsa wood, carbon rod, hard wood
Mass	5.1 g	11.6 g

Table 4.

Wing models.

3.3 Results

Experimental results showed that a flutter phenomenon appeared with the non-structure wing (broken wing in Figure 10c), but this phenomenon did not happen with the structure wing model. It could be explained by the more durability of structure wing than that of the non-structure wing with the testing range of velocity. Experiments also demonstrated that the combination of multiple materials to more durability of structure of wing could be highly effective in preventing flutter phenomenon [18].

The measurement results of the non-structure wing were shown in Table 5. When the attack angle increased, the velocity of first oscillation, flutter velocity, and frequency decreased.

Attack angle (^0)	Velocity of first oscillation (m/s)	Flutter velocity (m/s)	Frequency (Hz)
0	20	22	37.38
5	17.8	20.5	32.22
10	15.2	19.5	30.31

Table 5.

Flutter characteristics at different attack angles—Non-structure wing.

Figure 11 resumed the measurement of the force at the wing root in varying velocities from zero to flutter velocity and more by using the load cell system. After increasing the air velocity from zero to the limit of non-structure wing, the limit velocity of non-structure wing was found at 19.5 m/s. The load at the wing root of this wing at flutter velocity is shown in Figure 12. The maximum force was 5.44 N, and the minimum force was 4.32 N.

Figure 11.

Force at the wing root of the non-structure wing—Attack angle 10°.

Figure 12.

Force at the wing root of the non-structure wing—Attack angle 10° and velocity 19.5 m/s.

With structure wing, different modes of vibration appeared depending on the characteristics of the structure as remarked in [8]. With the help of the oscillator generator system, the specific oscillation frequencies of the first four modes were estimated as shown in Table 6.

Mode	Frequency (Hz)
	Non-structure wing	Structure wing
1	24.3	23.5
2	67.0	63.0
3	103.6	115.3
4	132.0	139.0

Table 6.

Specific oscillation frequency of the non-structure and structure wings.

The frequencies of the first and second modes of the non-structure wing were higher than those of the structure wing, while the frequencies of third and fourth mode of structure wing were higher than those of non-structure wing (Table 6). Considering both wings in the first mode of oscillation, when the force was applied to the wing, the amplitude of the non-structure wing was higher than that of the structure wing (Figure 13). In conclusion, the non-structure wing was easier to resonate than the structure wing.

Figure 13.

Amplitude of oscillation of the structure wing—Mode 1.

4. IBM method

4.1 Methodology

Fluid flow and deformation of structure were governed by the following equations with assumption of linear elastic structure [11]:

E5

E6

E7

E8

E9

E10

where:

u: fluid velocity vector

p: fluid pressure

f: force that affected on wing

Re: Reynolds number

U_c: displacement velocity

ω_p: angular velocity

x_c: center of gravity

θ_p: rotation of wing r

mp: mass of wing

Ip: inertial moment of wing

F: force created by fluid passing through the wing

T: moment created by fluid passing through the wing

To solve out these equations using IBM method, the most important was that the velocity of fluid at fluid-solid interface was equal to the velocity of the wing. It means that the interaction force (f) between the wing and fluid was calculated such that the boundary condition of fluid was satisfied on the surface of the wing.

IBM method used the Cartesian grid and immersed boundary that were illustrated in Figure 14, in which the moving surface of wing was described by Lagrangian points (rounded points) and fixed points in fluid were called Eulerian points. Parameters of Lagrangian points were noted as capitalization.

Figure 14.

Cartesian grid and immersed boundary.

Discrete partial derivative of velocity over time in Eq. (5) with denote intermediate velocity at zero force of Lagrangian point, Û_k, force F of Lagrangian points were estimated as follows:

E11

where:

k: time step.

U_k+1: identified from the moving surface of wing, so this velocity was known as U^(b).

Force was created from displacement and affected on the fluid element. Force was calculated using the following interpolation formula:

E12

where:

x: coordinate of Eulerian point.

N: set of Lagrangian points round Eulerian point l.

x_l: coordinate of Lagrangian point.

∆U_l: volume of effect corresponding to Lagrangian point l.

δ_h: 3D delta function identified as follows:

E13

E14

$$\phi \left(r\right)=\left\{\begin{array}{l}\frac{1}{6}\left(5-3\left|r\right|-\sqrt{1-3{\left(1-\left|r\right|\right)}^2}\right)\frac{1}{2}\le \left|r\right|\le \frac{3}{2}\\ \frac{1}{3}\sqrt{1-3{\left|r\right|}^2}\left|r\right|\le \frac{1}{2}\\ 0\left|r\right|\ge \frac{3}{2}\end{array}\right.$$E15

Three-step Runge–Kutta method was applied to solve Navier-Stokes equations and Newton equations as follows:

• Step 1: Calculate the instantaneous velocity at Eulerian points with no immersed boundary surface, i.e., f = 0:

E16

where:

k: step calculation of Runge-Kutta method (k = 1, 2, 3)

α_k: coefficient of k^th step calculation

γ_k: coefficient of k^th step calculation

ζ_k: coefficient of k^th step calculation

υ: kinematic viscosity

Apply this instantaneous velocity to calculate Lagrangian velocity on the surface of the wing:

E17

This Lagrangian velocity was combined with wing velocity, U^(b)_xl, which was calculated from the dynamic equation of wing, to calculate forces of Lagrangian points (F) following Eq. (11). Then, apply this force of Lagrangian points to Eq. (12) to calculate force of Eulerian points (f).

• Step 2: Solve out Navier-Stokes Eqs. (5) using calculated forces of Eulerian points to estimate the effect of flutter of the wing into the velocity field around the wing:

E18

To satisfy the continuity equation, a temporary pressure was described:

E19

• Step 3: Solve Eq. (17), and calculate velocity and pressure at k^th step of the Runge-Kutta method:

E20

E21

From the estimated forces at Lagrangian points, translational and angular movements of the wing were carried out by solving Eqs. (7) and (9):

E22

E23

After calculating the velocity of center of gravity (U_c^k) and angular velocity of wing (ω_p^k), coordinates of Lagrangian points were estimated by the same expressions.

4.2 Wing model

Four different wing models were carried out in order to analyze the effect of the wing structure. There were two rectangular and two trapezoid 3D-shape wings, and each 3D-shape wing had symmetric (NACA65A004) and asymmetric (supercritical) airfoil, respectively. The wings had the same area of 450 cm² and the same semispan-wise length of 30 cm. Therefore, the rectangular wing had a chord length of 7.5 cm, while the trapezoidal wing has no sweep angle, and the leading edge line had a tip chord length of 5 cm and root chord length of 10 cm. The wings were made of aluminum (Figure 15).

Figure 15.

Wing models. (a) Rectangular wing. (b) Trapezoidal wing.

Experiments were performed using a low-speed blowdown wind tunnel, which belongs to the Department of Aerospace Engineering at the Hanoi University of Science and Technology, Vietnam. This wind tunnel had a maximum free-stream velocity in empty test section of 30 m/s that corresponded to Reynolds number 10⁶. The wind tunnel was operated continuously by an 8 kW electric motor. The turbulence level in test section was slightly less than 1%. Free-stream velocity was kept constant in test section within ±2%. Total pressure of free-stream and dynamic pressures were measured by pitot tube within ±2%. Ambient temperature was measured within ±1%. Both experimental and numerical researches were performed at air velocity of 20 m/s and attack angle of 5°.

For the experimental study, 160 pressure taps were applied on the wing model (Figure 15). These pressure taps were connected to an external digital manometer via stainless and silicon tubes. Each pressure tap was measured one time with waiting time of 5 s (average of about 1000 instant values) using the Keyence pressure measurement. The standard deviation of the Keyence pressure measurement errors was within ±0.001 Pa. Moreover, flutter of wing was captured with help of high performance camera.

4.3 Results

The results of IBM method were analyzed at three different instants (Figure 16):

• Time T₀: initial time when distortion did not occurred

• Time T₁: time between maximum deformation and non-deformation

• Time T₂: time of maximum distortion

Figure 16.

Instant displacement of wings. (a) NACA65A004 – Rectangular wing. (b) NACA65A004 – Trapezoidal wing. (c) Supercritical – Rectangular wing. (d) Supercritical – Trapezoidal wing.

The wing deformation was maximum at the tip of the wing and decreased gradually into the root of the wing over time. However, the normal stress was found to have an opposite tendency in comparison with deformation. The maximum normal stress was observed at the root of the wing, while the minimum normal stress was found at the tip of the wing (Figures 17 and 18). It could be explained by the fixed support with fuselage at the root of the wing and free support at the tip of the wing [13]. These remarks were in well agreement with the experimental results within a relative error less than 10% (Table 7).

Figure 17.

Instant normal stress—Rectangular wing.

Figure 18.

Instant normal stress—Trapezoidal wing.

Wing	IBM method (mm)	Experiment method (mm)	Relative error (%)
NACA65A004-rectangular	0.119	0.131	9.7
NACA65A004-trapezoidal	0.030	0.033	8.6
Supercritical-rectangular	0.035	0.039	9.9
Supercritical-trapezoidal	0.034	0.037	9.2

Table 7.

Maximum deformation of the wing tip.

At T0 instant, the normal stress had important value near the wing tip. During flutter behaviors of wing, this important normal stress propagated from the tip of the wing to the root of the wing. The maximum value of the normal stress was found out at the root of the wing and at T₂ instant.

With the same airfoil, the rectangular wing was found to be more distorted and have higher maximum deformation and higher maximum normal stress than the trapezoid wing. Thus, 3D-shape wing contributed significantly to the deformation of wing when aeroelasticity phenomenon occurred (Table 7).

With the same 3D-shape wing, the maximum deformation and maximum normal stress of NACA65A004 rectangular wing were higher than those of the rectangle supercritical wing. Meanwhile, the maximum deformation and maximum normal stress of NACA65A004 trapezoidal wing were less than those of the supercritical trapezoidal wing. It could be concluded that the 3D shape of wing played an important role in the durability of the structure (Table 7).

5. Conclusions

The flutter phenomenon of AGARD 445.6 wing was determined by (a) a modal approach for a structural response; (b) an aerodynamic damping coefficient to predict the appearance of flutter phenomenon; (c) a strongly coupled FSI method to predict the aeroelastic response for subsonic and transonic flutter characteristics; (d) an experiment method to predict the aeroelastic response for subsonic flutter characteristics with wing structure; (e) and an IBM method to improve the interface between the fluid and solid of aircraft wing.

In brief, the major results could be summarized as follows:

• Experimental results were in good agreement with numerical results within a relative error less than 10%.

• During aeroelasticity phenomenon, deformation of the wing tip was maximum while it was minimum at the wing root. The tendency of normal stress was in contrast with deformation. The minimum normal stress was observed at the wing tip, while the maximum normal stress was observed at the wing root;

• Geometry of wing (3D shape, airfoil) had a significantly contribution to the deformation of wing when aeroelasticity phenomenon occurred.

For further research of aeroelasticity in the future, both experimental and numerical researches at low and high speed should be performed.

Acknowledgments

A part of this work was supported by the Ministry of Science and Technology in Vietnam through the bilateral and multilateral research project HNQT/SPĐP/12.19.

Appendices and nomenclature

w

structural displacement at any time instant and position

q

generalized displacement vector

[M]

generalized mass matrice

[C]

damping matrice

[K]

stiffness matrices

ϕ

normal modes of the structure

N

total number of modes of the structure

F

generalized force vector

X_i

i^th peak of vibration

n

number of period

ζ

aerodynamic damping coefficient

M

Mach number

u

fluid velocity vector

p

fluid pressure

f

force that affected on wing

Re

Reynolds number

U_c

displacement velocity

ω_p

angular velocity

x_c

center of gravity

θ_p

rotation of wing

m_p

mass of wing

I_p

inertial moment of wing

F

force that created by fluid go pass through the wing

T

moment that created by fluid go pass through the wing

X

coordinate of Eulerian point

N

set of Lagrangian points round Eulerian point l

x_l

coordinate of Lagrangian point

∆U_l

volume of effect corresponded to Lagrangian point l

δ_h

3D delta function

α γ ζ

coefficient of step Runge-Kutta calculation

υ

kinematic viscosity

E₁₁

normal stress following x-coordinates

E₂₂

normal stress following y-coordinates

E₃₃

normal stress following z-coordinates

G

shear stress

ρ

density

Cp

coefficient of pressure

α

attack angle