Approximation of Distributed Aerodynamic Force to a Few Concentrated Forces for Studying Supersonic Panel Flutter

초고속 패널 플러터 연구를 위한 분포 공기력의 집중하중 근사화

Abstract

The present study considers the usage of concentrated forces to simulate real panel flutter. The concept of using concentrated forces have been validated for studying the flutter of wing structure in subsonic flow, yet its application in the supersonic region remained to be explored. Hence, a simply supported panel subjected to forces, equivalent to aerodynamic force is considered for studying supersonic panel flutter. The distributed aerodynamic forces are approximated to few concentrated forces by taking numerical integration. The aeroelastic equation is formulated using the classical small-deflection theory and the piston theory for linear panel flutter whereas for emulated panel flutter the flutter equation is derived by replacing the pressure due to aerodynamic loading with pressure from concentrated loading. Finally, flutter frequency, flutter dynamic pressure, and corresponding mode shape are found for emulated panel flutter and compared with linear panel flutter. Two important parameters, the number of concentrated forces and their location are discussed through numerical examples and optimization process respectively. So far, the flutter results acquired in this study are reasonable to suggest the feasibility of reproducing panel flutter using concentrated forces.

이 연구는 패널 플러터 시뮬레이션을 위한 집중 하중의 사용을 연구한다. 이러한 구상은 날개 구조의 아음속 플러터 연구에 대해서는 검증된 바 있으나 초음속 영역에서는 그렇지 못하다. 따라서, 4면 단순 지지 경계 조건의 패널에 공기력과 등가의 집중하중을 가하여 초음속 패널 플러터를 연구한다. 분포된 공기력은 수치 적분 계산을 통해 집중 하중들로 근사된다. 선형 패널 플러터에 대한 공탄성 방정식은 고전적인 small-deflection theory와 piston theory를 이용하여 세워지는 반면, 모방된 패널 플러터에서 플러터 방정식은 분포 공기력에 의한 압력을 집중 하중에 의한 압력으로 대체함으로써 유도된다. 최종적으로 플러터 주파수, 임계 동압, 그리고 그에 상응하는 모드형상이 모방된 패널 플러터에 대해 구해지고, 그 결과를 선형 패널 플러터로부터 얻은 결과와 비교하여 검증하였다. 또한 두 가지 중요한 파라미터인 집중 하중의 개수와 위치는 수치적 예제들과 최적화 과정을 통해 각각 논의되었다. 이 연구에서 얻어진 플러터 결과는 집중하중들을 이용하여 패널 플러터를 재현하는 가능성을 논의하는데 타당한 것으로 생각된다.

1. Introduction

The panel under air loads can cause self-excited vibration due to the mutual interaction among elastic, inertial, and aerodynamic forces. For increasing dynamic pressure, at flutter onset, a catastrophic panel failure may occur. Therefore, to avoid such an undesirable phenomenon, flutter analyses in supersonic flow have been performed using analytical, finite element, and experimental methods to set the flutter boundaries. A comprehensive review of linear and nonlinear panel flutter have been published by Dowell(1) and Mei et al.(2). All the literature on panel flutter has considered distributed aerodynamic loading. Concerning experimental studies, so far, flutter wind tunnel testing has been carried out to investigate the aeroelastic instabilities of an aircraft. However, in general, the testing is expensive and complex, and unaffordable for general research propose. Similarly, the structure discrepancy from the scaled down model could result in flutter boundaries error.

An alternate way of treating flutter without the need for wind tunnel, known as "Dry Wind Tunnel"(DWT)(3) based on ground vibration test (GVT) has been proposed. The DWT system utilizes GVT hardware for analyzing flutter boundaries of a full-scale structure subjected to concentrated forces. The DWT system has been demonstrated on a rectangular plate, clamped at one edge for a subsonic Mach number and a very good agreement with their finite element method counterparts. However, the flutter characteristics in supersonic Mach number such as panel flutter using DWT concept has not been investigated. And, it is to understand that the experimental procedure is far more complicated for supersonic flow than subsonic flow.

In the present paper, the flutter boundaries of a simply supported panel in supersonic Mach number have been analyzed using distributed aerodynamic forces and a few concentrated forces respectively. The key challenge of reducing distributed aerodynamic loads to a few concentrated forces is carried out through the numerical integration of pressure distribution based on the piston theory and finding the best-concentrated forces location through an optimization approach.

 

2. Linear Panel Flutter

2.1 Mathematical Formulation

A flat panel with one side exposed to supersonic airflow is idealized to be uniform, thin, and isotropic as shown in Fig. 1. The structural formulation of the panel is carried out using the classical small-deflection theory and the aerodynamic formulation is based on the piston theory. Thus, the governing equation of panel motion is written as(4)

Fig. 1A flat panel

where, ΔP is an aerodynamic pressure due to panel motion and free stream velocity and in a direction opposite to transverse deflection, expressed as(2)

The solution to linear panel flutter problem is based on modal expansion. The panel deflection (w) is expanded in terms of its structural normal modes, Wmn, for example, Wmn is an eigenfunction and ωmn is a natural frequency from free vibration analysis of Eq. (1) that satisfy same boundary conditions as w. The modal expansion of plate is written as

Then, substituting Eq. (3) and Eq. (2) into Eq. (1), the system of continuous equations to the discrete equations for qmn is obtained by introducing orthogonality condition and Galerkin's procedure of multiplying resultant equation by Wuv and integrating over the plate surface.

Where,

and,

where, Quv is a generalized force. The Eq. (4) is mass normalized and expressed in matrix notation.

where, is a modal mass matrix, is a modal aerodynamic damping matrix, is a modal stiffness matrix, is a modal aerodynamic stiffness matrix, , and g = λ(M2 - 2) /U∞(M2-1). Apparently, Eq. (6) is a second-order linear equation which does not allow direct eigenvalue analysis, so the transformation of variables defined by q1 = q and = are used to simplify Eq. (6) in a standard first-order form.

If one considers, qmn(t) = Amnept, a form of solution as simple harmonic time dependence, the resulting eigenvalues ( p = pR + ipI) from Eq. (7), in general, are complex that determines the stability of the system. So, while the real part of the eigenvalue is negative, the system is stable and if the real part is positive the system is unstable. Similarly, as the real part becomes zero, the system is considered neutral, a point, generally referred to flutter onset at which critical flutter velocity and critical flutter frequency are acquired. Non-dimensional parameters are introduced to generalize the results. They are non-dimensional dynamic pressure , non-dimensional frequency , and damping .

2.2 Convergence Test

The flutter characteristics of typical aluminum rectangular simply supported panel (a/b = 1, a/h = 100) is presented in Table 1. The following properties are used for the convergence test: E = 70 GPa, ρ = 2700 kg/m3, and μ = 0.3. A simply supported panel has a deformation shape in a form of traveling wave with exponential spatial growth(4). In such case, a large number of modes are required to obtain accurate results. Therefore, to specify the necessary number of modes, a convergence test is preferable. The test is carried out for various combination of modes given by the number of half sine waves in the flow direction (m) and perpendicular to the flow direction (n) as shown in Fig. 4. It is observed that for m > 5 the flutter boundary values are converging well irrespective of a number of half sine wave perpendicular to the flow direction. Therefore, to obtain quantitatively accurate results, 24 modes for m = 8 and n = 3 are selected.

Table 1Comparison of flutter boundary for simply supported square panel

Fig. 4Flutter convergence test of simply supported panel

Figure 2 illustrate eigenvalue plot considering 24 modes in which the non-dimensional pressure parameter, λ* = 0, represent free vibration solution. While for λ* > 0, there exist an infinite number of coalescence of two structural modes. However, the lowest value of λ* at which the first flutter occur is physically accepted. Thus, the first flutter is observed at and . Similarly, Figure 3 illustrates the evolution of mode shape of first two dominant modes when the dynamic pressure is increased until mode 11 and mode 21 become single flutter mode at . This shows that the lateral deflection is concentrated at the rare end of the panel. The flutter bounds obtained from the present study are compared with the results from the available literature shows a very good agreement as listed in Table 1.

Fig. 2Eigenvalues plot of 1st and 2nd mode

Fig. 3The evolution of mode 11 and mode 21and the flutter mode shape

 

3. Emulated Panel Flutter using Concentrated forces

The emulating of real panel flutter by replacing distributed aerodynamic forces with concentrated forces has not been researched previously. Practically, the main application of concentrated forces is to simulate aerodynamics. The dynamic shakers can be controlled to provide necessary concentrated force to mimic aerodynamics in order to determine flutter speed and flutter frequency. This study could be significant to circumvent complicated wind tunnel testing as well as applicable for general research propose. Meanwhile, the viability of approximating distributed aerodynamic loads to few concentrated forces through numerical integration method is considered.

3.1 Mathematical Formulation

An attempt is made to replace aerodynamic forces with a few concentrated forces using numerical integration method. We assumed that the total force contributed by a number of concentrated forces is same as that of the distributed aerodynamic forces respectively. In other words, the resultant force from numerical integration of distributed aerodynamic pressure for entire panel and the sum of k number of concentrated force is equivalent. And, the distributed aerodynamic pressure is evaluated using the piston theory. Thus, the total resultant force is expressed as

Similarly, the work done by each of the concentrated forces are proportional to displacement it creates at unknown coordinates (xk, yk) in the panel. The displacement at unknown coordinates are calculated using normal modes (Wuv) with Dirac delta function. Thus the mathematical expression of work done by the concentrated forces is the product of Eq. (8) and displacement given by normal modes at the location of concentrated forces given as

where, Dirac delta function: δ ( x - xk, y - yk ) = 0, x ≠ xk, y ≠ yk

Furthermore, for a system of equations, the total work done by the concentrated forces can be expressed in matrix notation given by Eq. (10) where the first term is approximated aerodynamic stiffness matrix and the second term is approximated aerodynamic damping matrix due to concentrated forces.

Depending upon the selection of m and n, the mode indicators, the matrices in Eq. (10) are expressed as follows for each concentrated forces.

Additionally, the elements in above matrices for a fully simply supported panel are expressed in terms of normal modes.

Thus, the governing equation of plate motion for concentrated loading is written as

The structural formulation is similar to linear panel flutter whereas formulation of pressure due to concentrated forces is deduced from Eq. (10) and the mass normalized flutter equation in modal coordinate is constructed.

The following equation is in second-order and does not allow direct eigenvalue solution so it is reduced to standard first-order form by considering the transformation of variables, q1 = q and .

Further analysis in obtaining eigenvalues require assuming a form of solution as simple harmonic. Then, the flutter characteristics are determined by the real part of the eigenvalues as explained in linear panel flutter section.

There are two important parameters that require special attention in order to make current work significant. Firstly, it is choosing the number of concentrated forces (k) and secondly, their locations (xk, yk). A number of concentrated forces need to be few for practical application. Meanwhile, they should be sufficient to perform the flutter test. Therefore, the best location of concentrated forces is determined by a sequence of location optimization.

3.2 Optimization of Locations

The accuracy of flutter boundary depends on the locations of concentrated forces on the panel. Therefore, an optimization is carried out to determine the best locations of concentrated forces that would result in a close solution with linear panel flutter. The optimization is based on the inbuilt optimization algorithm of MATLAB program to minimize the given objective(3).

where, and are referential non-dimensional flutter dynamic pressure and frequency, and are non-dimensional flutter dynamic pressure and frequency for optimized force location.

The location optimization process uses gradient-based methods for minimizing an objective. The optimization process is a function of f, initial guess coordinates (xi, yi), the lower bound (LB), upper bound (UB), and constraints (C) as expressed in Eq. (15). At first, a set of initial guess coordinates of concentrated forces is selected. Each of these coordinates has a freedom of movement within the panel defined by a lower bound and upper bound. Besides that, the constraints between concentrated forces locations are useful in reducing the computation time. Then, the optimization process examines for the best locations within the sub-panels defined by the lower and upper bounds by performing flutter analysis using Eq. (13) to minimize the objective function.

As referred to DWT, the flutter test was accomplished on a rectangular plate with one edge clamped and remaining edge free using four concentrated forces and considering first four modes. But in the present study, for simply supported boundary condition, the structural modes take a form of traveling waves and the comparatively large number of modes are required to obtain quantitatively accurate results. Hence, compromising between the concentrated forces number and the necessity to obtain meaningful results, more than four concentrated forces are desirable.

 

4. Results and Discussions

An aluminum square panel with simply supported boundary is selected for numerical analysis. The studies are conducted with three different panel configuration as shown in Fig. 5 containing four, six, and nine concentrated forces respectively. These concentrate forces are marked as k = 1, 2, … 9. Similarly, the panel is divided into k equal sub-panels for a respective number of concentrated forces. Dividing panel into a number of sub-panel is an efficient way to obtain concentrated forces as well as to perform location optimization. The non-dimensional dynamic pressure and frequency are obtained by solving eigenvalue problem given by Eq. (13) for 24 modes (m = 8 and n = 3) using MATLAB program.

Fig. 52 × 2, 3 × 2, and 3 × 3 panel configurations

Figure 6 shows the eigenvalue plot and a flutter mode shape for three different panel configurations where the concentrated forces are placed at the center of sub-panel. Using four concentrated forces, no flutter behavior is observed. So, the number of concentrated forces are increased along the flow direction. With six concentrated forces, a frequency coalescence between first two natural modes is observed. Again, more concentrated forces are added in perpendicular to flow direction. With the use of nine concentrated forces, the coalescence of frequencies is also observed.

Fig. 6Emulated panel flutter analysis of panel configuration 2 × 2, 3 × 2, and 3 × 3 (m = 8, n = 3)

In order to determine emulated panel flutter bounds that are closer to reference linear panel flutter bounds, an optimization is carried out to find the best locations of six and nine concentrated forces as shown in Fig. 7. Figures 8 and 9 shows an eigenvalues plot for both distributed loading and concentrated loading using six and nine concentrated forces respectively. It can be observed that there is a good convergence of emulated flutter results for 3 × 2 panel configuration using six concentrated forces. However, using nine concentrated forces with 3 × 3 panel configuration a noticeable difference between linear and emulated panel flutter bounds is observed. Similarly, the objective function from both configurations, taking forces at the center of sub-panels and with optimized locations are presented in Fig. 10. Such that, a minimum objective function is observed for the six concentrated forces.

Fig. 7Location optimization of concentrated forces

Fig. 8Emulate flutter analysis with optimized location using 6 concentrated forces (3 × 2 configuration)

Fig. 9Emulate flutter analysis with optimized location using 9 concentrated forces (3 × 3 configuration)

Fig. 10Objective function

 

5. Conclusion

The numerical integration method to emulate real panel flutter is performed for a simply supported square panel. The aerodynamic loading was based on the piston theory and structural formulation was based on the classical small-deflection theory. The results propose that emulating real panel flutter using concentrated forces is feasible. Unlike wing flutter, panel flutter requires a large number of modes for calculation. Hence, different panel configuration can be approached to acquire concentrated forces in order to perform parametric studies.

 

Nomenclature

a, b: Panel length and width h: Panel thickness E: Panel young's modulus μ: Poisson's ratio D: Flexural rigidity, Eh3/12(1-μ2) ρ: Panel density w(x,y,t): Panel transverse deflection m, n, u, v: Mode indicators Wmn, Wuv: Normal modes ΔP: Aerodynamic pressure ΔPc: Pressure due to concentrated loading Quv: Generalized aerodynamic force p: Complex eigenvalue, pR+ipI k: Number of concentrated forces (xk,yk): Location coordinate of concentrated forces Smnuv: Aerodynamic stiffness coefficient due to concentrated forces Dmnuv: Aerodynamic damping coefficient due to concentrated forces f: Objective function