**4. IBM method**

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,

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

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

 24.3 23.5 67.0 63.0 103.6 115.3 132.0 139.0

**Non-structure wing Structure wing**

**Mode Frequency (Hz)**

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

*Amplitude of oscillation of the structure wing—Mode 1.*

and the minimum force was 4.32 N.

estimated as shown in **Table 6**.

*Aerodynamics*

**Figure 12.**

**Table 6.**

**Figure 13.**

**66**

## **4.1 Methodology**

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

$$\frac{\partial \boldsymbol{u}}{\partial t} = \frac{1}{\text{Re}} \nabla^2 \boldsymbol{u} - \boldsymbol{u} \nabla \boldsymbol{u} - \nabla p + \boldsymbol{f} \tag{5}$$

$$
\nabla \iota = 0 \tag{6}
$$

$$\frac{d\left(m\_p U\_c\right)}{dt} = F + \mathbf{g} \tag{7}$$

$$\frac{d\mathbf{x}\_c}{dt} = U\_c \tag{8}$$

$$\frac{d\left(I\_p \cdot \mathcal{O}\_p\right)}{dt} = T\tag{9}$$

$$\frac{d\theta\_p}{dt} = \phi\_p \tag{10}$$

where: u: fluid velocity vector p: fluid pressure f: force that affected on wing Re: Reynolds number Uc: displacement velocity ωp: angular velocity xc: 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

*φ*ð Þ¼ *r*

and Newton equations as follows:

where:

of the wing:

the wing:

Runge-Kutta method:

**69**

1

8 >>>>>><

*DOI: http://dx.doi.org/10.5772/intechopen.91748*

>>>>>>:

immersed boundary surface, i.e., f = 0:

αk: coefficient of kth step calculation γk: coefficient of kth step calculation ζk: coefficient of kth step calculation

to calculate force of Eulerian points (f).

υ: kinematic viscosity

1 3 q

0 j j*r* ≥

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

<sup>6</sup> <sup>5</sup> � <sup>3</sup>j j*<sup>r</sup>* �

*Research on Aeroelasticity Phenomenon in Aeronautical Engineering*

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 � 3j j*r* 2

> 3 2

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 � 3 1ð Þ � j j*r*

� � q 1

j j*r* ≤ 1 2

Three-step Runge–Kutta method was applied to solve Navier-Stokes equations

Apply this instantaneous velocity to calculate Lagrangian velocity on the surface

This Lagrangian velocity was combined with wing velocity, U(b)x*l*, 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)

• 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

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

• Step 3: Solve Eq. (17), and calculate velocity and pressure at kth step of the

• Step 1: Calculate the instantaneous velocity at Eulerian points with no

2

<sup>2</sup> <sup>≤</sup>j j*<sup>r</sup>* <sup>≤</sup>

3 2

(15)

ð16Þ

ð17Þ

ð18Þ

ð19Þ

ð20Þ

**Figure 14.** *Cartesian grid and immersed boundary.*

Lagrangian points (rounded points) and fixed points in fluid were called Eulerian points. Parameters of Lagrangian points were noted as capitalization.

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:

$$F = \frac{U^{k+1} - \hat{U}^k}{\Delta t} - \left(\frac{1}{\text{Re}}\nabla^2 U - U\nabla U - \nabla P\right) \tag{11}$$

where:

k: time step.

Uk+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:

$$f = \sum\_{l=1}^{N} F\left(\mathbf{x}\_{l}\right) \delta\_{h}\left(\mathbf{x} - \mathbf{x}\_{l}\right) \Delta V\_{l} \tag{12}$$

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:

$$\delta\_h \left( \mathbf{x} - \mathbf{x}\_l \right) = \delta\_h^{l,D} \left( \mathbf{x} - \mathbf{x}\_l \right) \delta\_h^{l,D} \left( \mathbf{y} - \mathbf{y}\_l \right) \delta\_h^{l,D} \left( z - z\_l \right) \tag{13}$$

$$\delta\_h^{1D}\left(r\right) = \frac{1}{h}\varphi\left(r\right) \quad r = \frac{\varkappa - \varkappa\_l}{h} \tag{14}$$

*Research on Aeroelasticity Phenomenon in Aeronautical Engineering DOI: http://dx.doi.org/10.5772/intechopen.91748*

$$\rho(r) = \begin{cases} \frac{1}{6} \left( 5 - 3|r| - \sqrt{1 - 3\left(1 - |r|\right)^2} \right) & \frac{1}{2} \le |r| \le \frac{3}{2} \\\frac{1}{3} \sqrt{1 - 3|r|^2} & |r| \le \frac{1}{2} \\\ 0 & |r| \ge \frac{3}{2} \end{cases} \tag{15}$$

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:

$$\tilde{u} = u^{k \ast 1} + \Delta t \left( 2\alpha\_k \nu \nabla^2 u^{k \cdot 1} - 2\alpha\_k \nabla p^{k \cdot 1} - \gamma\_k \left( (u \cdot \nabla) u \right)^{k \cdot 1} - \zeta\_k \left( (u \cdot \nabla) u \right)^{k \cdot 2} \right) \tag{16}$$

where:

Lagrangian points (rounded points) and fixed points in fluid were called Eulerian

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

Uk+1: identified from the moving surface of wing, so this velocity was known

Force was created from displacement and affected on the fluid element. Force

ð11Þ

ð12Þ

ð13Þ

ð14Þ

points. Parameters of Lagrangian points were noted as capitalization.

was calculated using the following interpolation formula:

N: set of Lagrangian points round Eulerian point *l.*

ΔU*l*: volume of effect corresponding to Lagrangian point *l.*

were estimated as follows:

*Cartesian grid and immersed boundary.*

where: k: time step.

where:

**68**

x: coordinate of Eulerian point.

x*l*: coordinate of Lagrangian point.

δh: 3D delta function identified as follows:

as U(b).

**Figure 14.**

*Aerodynamics*

k: step calculation of Runge-Kutta method (k = 1, 2, 3) αk: coefficient of kth step calculation γk: coefficient of kth step calculation ζk: coefficient of kth step calculation

υ: kinematic viscosity

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

$$\tilde{U}\_{\mathbf{x}\_{\parallel}} = \sum\_{k=1}^{N} \tilde{\mu}\_{\mathbf{x}} \delta\_{h} \left( \mathbf{x} - \mathbf{x}\_{k} \right) h^{3} \tag{17}$$

This Lagrangian velocity was combined with wing velocity, U(b)x*l*, 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:

$$\nabla^2 u^\ast - \frac{u^\ast}{\alpha\_k \nu \Delta t} = -\frac{1}{\nu \alpha\_k} \left( \frac{\tilde{u}}{\Delta t} + f\_\times \right) + \nabla^2 u^{k-1} \tag{18}$$

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

$$
\nabla^2 \Phi^k = \frac{\nabla \cdot \boldsymbol{u}^\*}{2\alpha\_k \Delta t} \tag{19}
$$

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

$$
u^{\!\!\!\!\/} = \boldsymbol{u}^{\!\!\!\!\/} - 2\alpha\_{\boldsymbol{k}} \Delta \boldsymbol{\!\!\!\/} \nabla \boldsymbol{\Phi}^{\!\!\!\/} \tag{20}$$

$$p^k = p^{k-1} + \Phi^k - \alpha\_k \Delta t \nabla^2 \Phi^k \tag{21}$$

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

*Research on Aeroelasticity Phenomenon in Aeronautical Engineering*

*DOI: http://dx.doi.org/10.5772/intechopen.91748*

*Instant displacement of wings. (a) NACA65A004 – Rectangular wing. (b) NACA65A004 – Trapezoidal*

*wing. (c) Supercritical – Rectangular wing. (d) Supercritical – Trapezoidal wing.*

**Figure 16.**

**71**

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

$$U\_c^k = U\_c^{k-1} + 2\alpha\_k \frac{\Delta t}{M} (F + G) \tag{22}$$

$$
\alpha \rho\_p^k = \alpha \rho\_p^{k-1} + 2\alpha\_k \, \frac{\Delta t}{I\_c} T \tag{23}
$$

After calculating the velocity of center of gravity (Uc k ) and angular velocity of wing (ω<sup>p</sup> 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<sup>2</sup> 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**).

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<sup>6</sup> . 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°.

**Figure 15.** *Wing models. (a) Rectangular wing. (b) Trapezoidal wing.*

ð21Þ

ð22Þ

ð23Þ

.

) and angular velocity of

From the estimated forces at Lagrangian points, translational and angular

), coordinates of Lagrangian points were estimated by the same

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<sup>2</sup> 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

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<sup>6</sup>

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

k

movements of the wing were carried out by solving Eqs. (7) and (9):

After calculating the velocity of center of gravity (Uc

performed at air velocity of 20 m/s and attack angle of 5°.

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

wing (ω<sup>p</sup>

**Figure 15.**

**70**

expressions.

*Aerodynamics*

k

**4.2 Wing model**

of aluminum (**Figure 15**).

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

#### **Figure 16.**

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