- Email: [email protected]
Projection-free proper orthogonal decomposition method for a cantilever plate in supersonic flow
Journal of Sound and Vibration 333 (2014) 6190–6208
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
Projection-free proper orthogonal decomposition method for a cantilever plate in supersonic flow Dan Xie a,n, Min Xu a, Earl H. Dowell b a b
College of Astronautics, Northwestern Polytechnical University, Xi'an 710072, China Duke University, Durham, NC 27708-0300, United States
a r t i c l e in f o
abstract
Article history: Received 5 May 2014 Received in revised form 25 June 2014 Accepted 25 June 2014 Handling Editor: L.G. Tham Available online 19 July 2014
This paper explores the use of the proper orthogonal decomposition (POD) method for supersonic nonlinear flutter of a cantilever plate or wing. The aeroelastic equations are constructed using von Karman plate theory and first-order piston theory. The twodimensional POD modes (POMs) in xy plane are determined from the chaotic results given by the traditional Rayleigh–Ritz (RR) approach. For a specific structure, the POMs need to be calculated once and then can be used for various parameters of interest. The derivatives of the POMs are calculated numerically to avoid the complex projection from the POMs to the Rayleigh–Ritz modes (RRMs). Numerical examples demonstrate that the POD method using 4 POMs can obtain accurate limit cycle oscillation (LCO) results with substantial computational cost savings, compared with 12 RRMs by the Rayleigh–Ritz method. The POD method is employed for the analysis of the chaotic oscillations. It is also demonstrated that the POD modes are robust over a range of flight parameters. & 2014 Elsevier Ltd. All rights reserved.
1. Introduction The aeroelastic system of a fluttering panel undergoing supersonic flow, using von Karman theory accounting for the structural nonlinearity and piston theory accounting for the aerodynamics, is a classical self-excited autonomous system, and has been intensively investigated by many researchers [1–7]. The mathematical model for this problem is a set of partial differential equations (PDEs), with both spatial and temporal variables. Analytical solutions to this set of PDEs are not available due to the existence of the structural nonlinearity which nullifies the additivity pertaining to the linear systems [8,9]. Nevertheless, considering the simple spatial domain (due to the simple rectangular geometry of the plate), several semi-analytical methods can be applied to solve this problem. For a simply supported plate [1,10,11], the Galerkin method has been used in the spatial domain to transform the governing PDEs into a set of coupled ordinary differential equations (ODEs) in time. Subsequently, the resulting ODEs were integrated by the numerical integration method. In addition, the Rayleigh–Ritz method was applied by Ye and Dowell [5] to analyze the LCO of a cantilever plate. For the same model, the chaotic oscillations are explored using the Rayleigh–Ritz method by Xie et al. [12]. However, when the plate has a complex geometry, the semi-analytical methods (the Galerkin method and the Rayleigh–Ritz method) become invalid. The finite element method (FEM) is more flexible for handling complex geometries. It has been employed extensively to study the nonlinear panel flutter [4,13,14]. However, the number of degrees-of-freedom (DOFs) required in the finite element analysis
n
Corresponding author. Tel.: þ86 18291964199. E-mail address: [email protected] (D. Xie).
http://dx.doi.org/10.1016/j.jsv.2014.06.039 0022-460X/& 2014 Elsevier Ltd. All rights reserved.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Nomenclature
R, S
a, b aij ; brs D, E h I, J
r, s t u, v uiðrÞ ; vjðsÞ V, vj W, w x, y, z
i, j K k Ma M, N m, n Δp Qk , Q q, q qk
plate length, plate width mode coordinate for u, v plate stiffness, Young's modulus plate thickness total mode number retained in x; y direction for u mode number in x; y direction for u total mode number retained for w (POD) POD mode number Mach number total mode number retained in x; y direction for w (RR) mode number in x; y direction for w aerodynamic pressure, positive in direction opposite to w generalized aerodynamic force, snapshots matrix dynamic pressure, POD mode coordinate in vector POD mode coordinate for w
β λ λpk μ
ν ξ; η ρ; ρm τ Φ Ψ, φ k φ~ k , φ k
6191
total mode number retained in x; y direction for v mode number in x; y direction for v time in-plane displacement in length and width mode in x; y direction for u(v) POD eigenvectors matrix, POD eigenvector w/h, panel transverse deflection streamwise, spanwise, normal coordinates ðMa2 1Þ1=2 , compressibility correction factor 2qa3 =βD, nondimensional dynamic pressure POD eigenvalue ρa=ρm h, nondimensional fluid/structure mass ratio Poisson ratio x=a; y=b air density, plate density 4 tðD=ρm ha Þ1=2 , nondimensional time correlation matrix POD modes matrix, POD mode in xy plane POD mode in vector form, POD mode in matrix form
is very large. Therefore, the semi-analytical methods are employed for the cantilever plate in this paper. Of course, one can use finite element analysis to find the eigenmodes of the structure and this is often done in aeroelastic analysis. All these traditional approaches use the assumed modes. Consequently, the coupled nonlinear modal equations with large number of modes are computationally costly for flutter analysis. For solving this problem, the reduced-order model (ROM) is proposed by many researchers. For the systems with original full-order equations of motion unknown, Attar [15] constructed a ROM using a system identification technique to investigate the LCO of a 451 delta wing. Guo and Mei [16,17] first used the aeroelastic modes for nonlinear panel flutter to reduce the number of ODEs and save computational costs greatly. Similarly, a set of POMs for the aeroelastic panels has been used to construct a ROM [18]. The POD method was also used to identify the coherent structures of the panel dynamics [19]. Furthermore, Amabili et al. [20,21] applied the POD method with a Galerkin projection (projection from the POMs to the Galerkin modes) for a water-filled circular cylindrical shell. More recently, a POD reduced-order model (POD/ROM) with a Galerkin projection has been constructed by Xie et al. [7] to solve the supersonic nonlinear oscillations of a simply supported plate with remarkable computational cost savings. However, seldom has a reduced-order model been used for a fluttering cantilever plate. In the present study, we create a reduced-order model using a computationally fast POD method to investigate the nonlinear flutter of a cantilever plate undergoing supersonic flow. The newly proposed POD method is different from the conventional POD method because it does not require the complex and computationally expensive projection from the POMs to the RRMs. The projection-free POD method has been applied for LCO analysis of a two-dimensional panel with all edges simply supported [22]. In this paper, LCO and chaotic oscillations are both investigated with the POD method. Section 2 provides the details for the construction of the POD/ROM. The convergence and efficiency studies are performed and the numerical examples are presented in Section 3. Moreover, a set of general POD modes is used for a range of physical parameters. Finally, some conclusions are drawn in Section 4.
2. Theoretical analysis A wing-like cantilever plate shown in Fig. 1 is constrained along the root edge with supersonic flow passing over the upper surface parallel to the fixed edge. The plate is taken as homogeneous isotropic aluminum with the material constants E ¼7.10 GPa, ρ ¼ 2:8 103 kg=m3 and ν ¼0.33. The von Karman plate theory and the first-order piston aerodynamic theory are used to account for the structural nonlinearity and the aerodynamic force, respectively. The modal equations being studied are derived directly from Hamilton's Principle and Lagrange's equations and thus the von Karman plate equations per se are never explicitly use. Of course the von Karman equations can be derived from Hamilton's Principle. The POD method using the beam modes for the in-plane displacements and the POD modes for the transverse displacement is employed to discretize the governing equations.
6192
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Y
Y
h
a
Flow
b
Z
X Top view
Side view
Fig. 1. Geometry of a cantilever panel or wing in top view and side view.
2.1. Numerical POD modes in xy plane It is the authors' experience in the present study as well as earlier studies [7] that using a single chaotic time history with broad band frequency response to construct the snapshots for determining the POD modes is much better than using a time history for a single frequency response. This is because the frequency spectrum of a single chaotic time series is broad and excites many POD or eigenmodes of the system. One could use many time histories at single frequencies for constructing a collection of snapshots as well, but a single chaotic time history is usually more computationally efficient. In this paper, the snapshots are gathered from the chaotic results by the Rayleigh–Ritz method. A similar procedure for the POMs' determination as previous work [7] is performed. First, the chaotic time response of the transverse deflection is assembled in a matrix as snapshots. Different columns correspond to the response at different times τ. Note that every column consists of WðτÞ of all the nodes in both x and y directions. We sort the order by fixing an individual η with ξ changed from ξ1 to ξN and then, fixing a subsequent η with ξ changed in the same way. This is different from the POMs in one-dimension [7]. Here the snapshots matrix takes the form as follows: 0
Wðξ1 ; η1 ; t 1 Þ B ⋮ B B B WðξN ; η1 ; t 1 Þ B B B Wðξ1 ; η2 ; t 1 Þ B B ⋮ B Q ¼B B WðξN ; η2 ; t 1 Þ B B ⋮ B B B Wðξ1 ; ηN ; t 1 Þ B B ⋮ @ WðξN ; ηN ; t 1 Þ
Wðξ1 ; η1 ; t 2 Þ
⋯
⋮
⋱
WðξN ; η1 ; t 2 Þ
⋯
Wðξ1 ; η1 ; t 2 Þ ⋮
⋯ ⋱
WðξN ; η1 ; t 2 Þ
⋯
⋮
⋱
Wðξ1 ; ηN ; t 2 Þ
⋯
⋮
⋱
WðξN ; ηN ; t 2 Þ
⋯
Wðξ1 ; η1 ; t J Þ
1
C C C WðξN ; η1 ; t J Þ C C C Wðξ1 ; η2 ; t J Þ C C C ⋮ C C; WðξN ; η2 ; t J Þ C C C ⋮ C C Wðξ1 ; ηN ; t J Þ C C C ⋮ A WðξN ; ηN ; t J Þ ⋮
(1)
where Wðξi ; ηj ; τk Þ denotes the non-dimensional transverse deflection of i-th point in the x direction and j-th point in the y direction at k-th time point ði ¼ 1; 2; …; N; j ¼ 1; 2; …; N; k ¼ 1; 2; …; JÞ. Note that J and N are the number of time points and nodes in the x or y direction, respectively. Then, we can get the correlation matrix: T
Φ¼Q Q:
(2)
Φvk ¼ λpk vk ;
(3)
And an eigenvalue problem is yielded:
where λk ; vk ðk ¼ 1; 2; …; JÞ are the eigenvalues and eigenvectors of Φ, respectively. Sorting the eigenvalues in a decreasing order, the eigenvectors vk are sorted correspondingly: p
λp1 ≧λp2 ≧⋯≧λpJ :
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
6193
A commonly used criterion for choosing K, the optimal number of basis retained in the POD subspace, is the so-called energy contribution criterion. For a predefined “energy percentage”, here ϵ ¼0.9999, K is chosen by ∑Kk ¼ 1 λk
p
F ðkÞ ¼
p ≧0:9999; ∑Jj ¼ 1 j
λ
and, then, the modal vectors of the correlation matrix is V ¼ ½v1 ; v2 ; …; vK . With the matrix V, the snapshots are linearly combined to generate a smaller number of basis vectors assembled in a matrix: 0 1 φ1 ðξ1 ; η1 Þ φ2 ðξ1 ; η1 Þ ⋯ φK ðξ1 ; η1 Þ B C ⋮ ⋮ ⋱ ⋮ B C B C B φ1 ðξN ; η1 Þ φ2 ðξN ; η1 Þ ⋯ φK ðξN ; η1 Þ C B C B C B φ1 ðξ1 ; η2 Þ φ2 ðξ1 ; η2 Þ ⋯ φK ðξ1 ; η2 Þ C B C B C ⋮ ⋮ ⋱ ⋮ C Ψ¼QV¼B (4) B C; B φ1 ðξN ; η2 Þ φ2 ðξN ; η2 Þ ⋯ φK ðξN ; η2 Þ C B C B C ⋮ ⋮ ⋱ ⋮ B C B C B φ1 ðξ1 ; ηN Þ φ2 ðξ1 ; ηN Þ ⋯ φK ðξ1 ; ηN Þ C B C B C ⋮ ⋮ ⋱ ⋮ @ A φ1 ðξN ; ηN Þ φ2 ðξN ; ηN Þ ⋯ φK ðξN ; ηN Þ ~ k ðk ¼ 1; 2; …; KÞ. where the k-th column of matrix (4) is corresponding to the POD mode vector φ 2.2. Aeroelastic equations with the POD modes Based on the mode functions for a free-free beam and a cantilever beam [23], we assume the mode functions for the inplane displacements of u and v: ui ξ ¼ cos iπξ;
2j 1 vj η ¼ sin πη; 2
ur ξ ¼ cos r πξ;
2s 1 vs η ¼ sin πη: 2
(5)
For the transverse deflection, the POD mode φk ðξ; ηÞ in xy plane is used. Using the panel thickness, a set of appropriate nondimensionalizations is employed: u ¼ u=h; v ¼ v=h and w ¼ w=h. Then the streamwise, spanwise and transverse displacements in the nondimensional form can be expressed in terms of the mode functions as follows: J
I
u ¼ ∑ ∑ aij ðτÞui ðξÞvj ðηÞ; i¼1j¼1 R
S
v ¼ ∑ ∑ brs ðτÞur ðξÞvs ðηÞ; r ¼1s¼1 K
w ¼ ∑ qk ðτÞφk ðξ; ηÞ:
(6)
k¼1
Based on the von Karman nonlinear strain–displacement relation, the normal strains ϵx ; ϵy in x; y directions and the shear strain γxy in xy plane take the forms as ∂u 1 ∂w 2 ∂v 1 ∂w 2 ∂v ∂u ∂w ∂w : (7) ϵx ¼ þ ; ϵy ¼ þ ; γ xy ¼ þ þ ∂x 2 ∂x ∂y 2 ∂y ∂x ∂y ∂x ∂y Then the elastic energies US and UB induced by the stretching and bending are Z Z Eh 1ν 2 US ¼ ϵ2x þ ϵ2y þ2νϵx ϵy þ γ xy dx dy; 2 2 2ð1 ν Þ UB ¼
D 2
" 2 2 #) 2 Z Z ( 2 ∂ w ∂2 w ∂2 w ∂2 w ∂ w þ 2 ð 1 ν Þ dx dy; ∂x∂y ∂x2 ∂y2 ∂x2 ∂y2
(8)
and T is the kinetic energy: T¼
1 2
Z Z
ρm h
∂w 2 dx dy: ∂t
(9)
6194
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Then L is the Lagrangian taking the form L ¼ T ðU S þ U B Þ, and Lagrange's equations are 8 ∂L > > ¼ 0; > > > ∂aij > > > < ∂L ¼ 0; ∂brs > > > > > d ∂L ∂L > > > : dt ∂ðdq =dtÞ ∂q þ Q k ¼ 0: k k Substituting L ¼ T ðU S þ U B Þ into Eq. (10) yields 8 ∂U S > > ¼ 0; > > > ∂aij > > > < ∂U S ¼ 0; ∂brs > > > h i ∂U > > d ∂U S B > ∂T > > : dt ∂ðdqk =dtÞ þ ∂q þ ∂q þQ k ¼ 0: k k Aerodynamic force is generalized as
Z Z Qk ¼
Δpðx; y; t Þ
∂w dx dy; ∂qk
where the aerodynamic pressure Δp is modeled using first-order piston or quasi-steady aerodynamic theory ! " # 2q ∂w Ma2 2 1 ∂w þ Δp ¼ : β ∂x Ma2 1 V 1 ∂t
(10)
(11)
(12)
(13)
Performing a subsequent substitution of equations as Eqs. (5) (- 6) (- 7) - (8) - (9), along with Eq. (13) - (12), and finally substituting them into Eq. (11) with appropriate nondimensionalization, the following governing equations are yielded: Ca a þCb b ¼ C; Da a þDb b ¼ D;
(14)
2
A
d q dt
2
þ Bqþ F þ Q ¼ 0:
(15)
Eq. (14) is a set of algebraic equations which can be solved for aij and brs in terms of qk, and Eq. (15) is a set of ordinary 2 differential equations in terms of aij ; brs ; qk ; dqk =dt and d qk =dt 2 . Pre-multiplying Eq. (15) by the inverse of matrix A yields 2
d q dt
2
þ A 1 Bq þ A 1 F þ A 1 Q ¼ 0:
(16)
This set of ODEs can be solved numerically by the fourth-order Runge–Kutta method (RK4). From Eq. (6), the panel displacements in time can be calculated. The matrices Ca ; Cb ; Da and Db are all the same as Ref. [5], hence they are omitted here for avoidance of repetition. However, using the numerical derivatives of the POMs, the other coefficients of C; D; A; B; Q and F take different expressions. Based on the mode expansion for the transverse deflection in Eq. (6) compared to Eq. (7) in Ref. [5], we can replace the terms in the expressions of C; D; A; B; Q and F in Ref. [5] as follows: Z Z Z Z M N M N K K ∑ ∑ ∑ ∑ qmn qpl ϕm ϕp dx ψ n ψ l dy- ∑ ∑ qk1 qk2 φk1 φk2 dx dy (17) m n
p
l
k1 k2
And then we can obtain the expressions of C; D; A; B; Q and F given in Appendix A for the present study. 2.3. Numerical derivatives of the POD modes Ye and Dowell [5] used the analytical beam modes for the transverse deflection. The derivatives of the beam modes can be calculated analytically and exactly. By contrast, the numerical POMs are used for the transverse deflection in this study. Hence, the differential manipulation on the POMs is performed numerically. Consequently, the differential and integral operations on the numerical POMs in the transformation of PDEs into ODEs may propagate the numerical errors. To circumvent this problem, a POD method projecting the POMs onto the assumed modes (RRMs or Galerkin modes) semi-analytically is employed [7,24]. In that way, any linear or nonlinear manipulation can be performed analytically to avoid any spatial errors. However, the authors’ experience has shown that the mode projection will lead to a complex math derivation of the governing equations and subsequently increase the computational costs vastly for the cantilever plate.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
6195
Table 1 Chaotic responses for the POMs determination. a=b
λ
Number of RRMs
It
Ts
Tc
Interval
0.5 1 2
90 500 1200
M ¼ 8, N ¼2, I ¼ R ¼ 14 M ¼ 8, N ¼2, I ¼ R ¼ 14 M ¼ 8, N ¼2, I ¼ R ¼ 16
t ¼ 500 t ¼ 200 t ¼ 200
0.005 0.005 0.002
0.4 0.3 0.3
5 2 5
Hence, the derivatives of the POMs are directly calculated by a three-point difference technique to avoid the mode projection in this paper. ~ k of k-th column in Eq. (4) into a matrix form, by sorting φk ðξ; ηÞ with different η First, we rewrite the POD mode vector φ into different columns, and different ξ into different rows: 0 1 φk ðξ1 ; η1 Þ φk ðξ1 ; η2 Þ ⋯ φk ðξ1 ; ηN Þ B C B φk ðξ2 ; η1 Þ φk ðξ2 ; η2 Þ ⋯ φk ðξ2 ; ηN Þ C C: φk ¼ B (18) B C ⋮ ⋮ ⋱ ⋮ @ A φk ðξN ; η1 Þ φk ðξN ; η2 Þ ⋯ φk ðξN ; ηN Þ And then the derivatives with respect to ξ and η can be calculated numerically by a three-point difference technique. The first and second derivatives with respect to ξ are performed in every column of matrix (18) as follows: ∂φk ðξ; ηÞ 3φk ðξi ; ηÞ þ 4φk ðξi þ 1 ; ηÞ φk ðξi þ 2 ; ηÞ ¼ ∂ξi 2dξ ∂φk ðξ; ηÞ φk ðξi þ 2 ; ηÞ φk ðξi ; ηÞ ¼ ∂ξi þ 1 2dξ
∂φk ðξ; ηÞ φk ðξi ; ηÞ 4φk ðξi þ 1 ; ηÞ þ3φk ðξi þ 2 ; ηÞ ¼ ∂ξi þ 2 2dξ ∂2 φk ðξ; ηÞ φk ðξi ; ηÞ 2φk ðξi þ 1 ; ηÞ þ φk ðξi þ 2 ; ηÞ ¼ ∂ξi ðdξÞ2 ∂2 φk ðξ; ηÞ φk ðξi ; ηÞ 2φk ðξi þ 1 ; ηÞ þ φk ðξi þ 2 ; ηÞ ¼ ∂ξi þ 1 ðdξÞ2
∂2 φk ðξ; ηÞ φk ðξi ; ηÞ 2φk ðξi þ 1 ; ηÞ þ φk ðξi þ 2 ; ηÞ ¼ ; ∂ξi þ 2 ðdξÞ2
i ¼ 1; 2; …; N 1:
(19)
In the same way, the derivatives with respect to η are performed in every line of matrix (18). Finally, substitute the numerical derivatives of the POMs into Eqs. (22), (23), (24), (25), (26) and (27) to calculate the coefficients of C; D; A; B; Q and F respectively. 2.4. Numerical integration for the aeroelastic equations Because the PDEs are transformed into low-order ODEs with the numerical POMs, the differential and integral operations on the numerical POMs result in a set of numerical ODEs. Hence, the integrations in the ODEs have to be solved numerically. A one-dimensional trapezoidal integration formula has been used extensively and takes the form " # Z b N1 h f ðxÞ dx ¼ f ðaÞ þ2 ∑ f ðxi Þ þ f ðbÞ ; 2 a i¼1 h¼
ba ; xi ¼ a þ i h; i ¼ 1; 2; …; N 1: N
(20)
In the present study, due to the two-dimensional POMs in xy plane, a corresponding two-dimensional trapezoidal integration formula is proposed and derived as follows: # Z bZ a Z b " N1 h f ð0; yÞ þ 2 ∑ f ðxi ; yÞ þ f ða; yÞ dy f ðx; yÞ dx dy ¼ 0 0 0 2 i¼1 # 2 (" N 1 h ¼ f ð0; 0Þ þ 2 ∑ f ðxi ; 0Þ þ f ða; 0Þ 2 i¼1 " # N1
þ2 ∑
j¼1
"
N1
f ð0; yj Þ þ 2 ∑ f ðxi ; yj Þ þ f ða; yj Þ i¼1
N1
#)
þ f ð0; bÞ þ 2 ∑ f ðxi ; bÞ þ f ða; bÞ i¼1
6196
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
0
4
λ=90
λ=90
3
−1 −2
1
dW/dτ
W(ξ,η,τ)
2
0 −1
−3 −4
−2
−5
−3 −4
0
100
200
300
400
−6 −0.25 −0.2 −0.15 −0.1 −0.05
500
τ
0
0.05
0.1
0.15
W
Fig. 2. Time history and Poincare plot for a panel of a/b ¼0.5 at λ ¼90: (a) time history; (b) Poincare map.
3
10
λ=500
λ=500 5
2
0
dW/dτ
W(ξ,η,τ)
1 0 −1
−10 −15
−2 −3
−5
−20
0
50
100
150
−25
200
−0.5
0
τ
0.5
W
Fig. 3. Time history and Poincare plot for a panel of a/b¼ 1 at λ¼ 500: (a) time history; (b) Poincare map.
60
2
λ=1200
1.5
λ=1200
40
dW/dτ
W(ξ,η,τ)
1 0.5 0 −0.5 −1
20 0 −20
−1.5 −2
0
50
100
150
200
−40
τ
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
W
Fig. 4. Time history and Poincare plot for a panel of a/b ¼2 at λ ¼ 1200: (a) time history; (f) Poincare map.
¼
) 2 ( N1 N1 h f ð0; 0Þ þ f ða; 0Þ þ f ð0; bÞ þ f ða; bÞ þ 4 ∑ ∑ f ðxi ; yj Þ 2 i¼1 j¼1
i N1 N 1h
þ 2 ∑ f ðxi ; 0Þ þ f ðxi ; bÞ þ 2 ∑ f 0; yj þ f a; yj ; i¼1
h¼
j¼1
a b or ; xi ¼ i h; yj ¼ j h; i; j ¼ 1; 2; …; N 1: N N
Here we let N be 100 in the present paper.
(21)
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
6197
Accumulative energy contribution %
100 99 98 97 96 a/b=0.5 a/b=1 a/b=2
95 94
0
2
4 6 POD mode number
8
10
Fig. 5. Accumulative energy of POD eigenvalues versus POD mode number for panels of a/b¼ 0.5, a/b¼ 1 and a/b¼ 2.
0.02
0.05
0.01
0
0 1 0.5 ξ
0 0
0.5 η
1
−0.05 1 0.5 ξ
0.05
0.05
0
0
−0.05 1
−0.05 1
0.5 ξ
0 0
0.5 η
1
0 0
0.5 ξ
0 0
0.5 η
0.5
1
1
η
a/b=0.5, λ =90 Fig. 6. The first 4 POD modes for a panel of a/b ¼ 0.5 determined at λ¼ 90 using 16 RRMs.
3. Results and discussions In this section, the nonlinear dynamics of a cantilever plate in supersonic flow shown in Fig. 1 are investigated, using the POD method free of the mode projection from the POMs to RRMs. The results are compared with those using the Rayleigh– Ritz method. For all the calculations here, h ¼0.01 m, Ma¼2, μ ¼ 0.1 and J¼ S¼2 unless otherwise noted. All the figures are plots with respect to a typical point ξ ¼0.75, η ¼1.0. The initial conditions are q1ðτ ¼ 0Þ ¼ 0:1 with others being zero. Here q1 is q11 in Eq. (6) and deflection Wp is the peak of w. Three panels of a/b¼0.5, 1 and 2 are considered and the POMs are determined from the Rayleigh–Ritz chaotic solutions. The property of the POMs is explored by comparing the POM shapes with the panel deflections. The convergence study is performed based on the modal participation of the POMs, and the efficiency is evaluated in terms of the computational costs. Finally, the chaotic oscillations are analyzed and the robustness of the POD modes to the system parameter variations is investigated. 3.1. Property of the POD modes The results in Ref. [7] suggest that the chaotic oscillations are the best choice to determine the POD modes, which contain the richest information of the dynamic systems. Thus, the general POD modes only need to be calculated once and can be
6198
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
0.02
0.05
0.01
0
0 1
−0.05 1
0.5 ξ
0 0
1
0.5
0.5
0.05
0.05
0
0
−0.05 1
−0.05 1
0.5 ξ
0 0
0 0
ξ
η
1
0.5
0.5 η
0.5 0 0
ξ
η
0.5
1
1
η
a/b=1, λ =500 Fig. 7. The first 4 POD modes for a panel of a/b ¼1 determined at λ ¼500 using 16 RRMs.
0.04
0.05
0.05
0.02
0
0
−0.05 1
−0.05 1
0 1
1
0.5 ξ
0 0
1
0.5
0.5 η
0 0
ξ
ξ
0.05
0.05
0.05
0
0
0
−0.05 1
1
0.5 ξ
0 0
0.5 η
−0.05 1
1
0.5 ξ
0 0
0.5 η
1
0.5
0.5 η
0 0
0.5 η
0 0
0.5 η
−0.05 1
1
0.5 ξ
a/b=2, λ=1200 Fig. 8. The first 6 POD modes for a panel of a/b ¼ 2 determined at λ ¼1200 using 16 RRMs.
2
2
1
1
0 1
(w/h)p
(w/h)p
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
0.5
0 0
0.5 η
0 1
1
1
0.5
0.5
0 1
0.5 ξ
0 0
0.5 ξ
(w/h)p
(w/h)p
ξ
1
1
0.5
0 1
0.5 ξ
η
0 0
0 0
6199
0.5
1
η
0.5 η
1
Fig. 9. Deflection shapes for panels of a/b¼ 0.5, 1, 2; λ¼ 50, 200, 950 are corresponding to LCO deflection shapes and λ ¼ 1200 is about chaotic deflection shape.
Table 2 Modal participation for a/b ¼ 0.5. λ
40 50 55
Mode
RRM POM RRM POM RRM POM
ðw=hÞp
MP (%) q1
q2
q3
q4
q5
q6
0.7825 0.7508 1.0632 1.0253 1.1682 1.1509
80.5228 85.4608 82.1996 89.0744 82.8710 90.8220
8.6812 8.3048 9.6046 5.9960 9.8578 5.1710
8.2011 6.0025 6.2903 4.7139 5.5860 3.7881
1.6080 0.1591 1.2246 0.1494 1.0922 0.1528
0.5959 0.0334 0.4136 0.0356 0.3581 0.0341
0.3910 0.0395 0.2673 0.0307 0.2348 0.0321
ðw=hÞp
MP (%) q1
q2
q3
q4
q5
q6
75.4875 79.8955 75.8395 86.3509 75.6475 86.7527
1.4745 8.9275 8.4764 6.1160 9.1233 6.1182
18.7835 7.5841 11.1256 5.3574 10.7588 5.1964
2.5438 3.1181 2.3365 1.6009 1.9005 1.4971
1.1994 0.3939 1.7402 0.5033 2.0711 0.3619
0.5112 0.0809 0.4819 0.0715 0.4988 0.0737
Table 3 Modal participation for a/b ¼ 1. λ
100 200 250
Mode
RRM POM RRM POM RRM POM
0.3053 0.2963 1.0392 1.0344 1.1689 1.1631
used for the various cases. Specifically, λ ¼90, 500 and 1200 are selected for the panels of a/b¼0.5, 1, and 2, respectively to generate the chaotic responses by the Rayleigh–Ritz method [12]. Herein we choose the parameters as shown in Table 1, where “It” (Integration time of RK4), “Ts” (Time step of RK4), “Tc” (Time chosen for POMs) and “Interval” are all the related time parameters in RK4 using the Rayleigh–Ritz method. a/b¼1, λ ¼500 is taken as an example, for which the time marching (RK4) is carried out from τ ¼ 0 to τ ¼200. Only the last 200 0.3 ¼60 is chosen for the POMs determination, and one solution is selected as a snapshot every 0.005 2 ¼0.01. Figs. 2–4 show the time responses and Poincare plots for the chaotic oscillations aforementioned. Based on Section 2.1, Fig. 5 shows the accumulative energy of the POD eigenvalues. Obviously, 4, 4 and 6 POMs are required for the panels of a/b¼0.5, 1 and 2, respectively. Subsequently, Figs. 6, 7 and 8 show the first four POM shapes for the panels of a/b¼0.5 and a/b¼1, and the first six POMs for the panel of a/b¼ 2, respectively. To demonstrate the POM property, Fig. 9 presents the
6200
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
1.4
12 RRMs 4 POMs
1
1.2
0.8
Wp
Wp
a/b=1
1
0.8 0.6
0.6
0.4
0.4
0.2 0
12 RRMs 4 POMs
a/b=0.5
0.2 20
25
30
35
40
45
50
55
0
60
50
100
150
200
250
300
λ
λ
Fig. 10. Comparison of LCO deflection peak versus dynamic pressure using 12 RRMs and 4 POMs: (a) a panel of a/b ¼0.5; (b) a panel of a/b ¼ 1.
1.5
2.5
12 RRMs 4 POMs
a/b=0.5, λ=50
2
1
1.5
0.5
0.5
12 RRMs 4 POMs
a/b=0.5, λ=50
dW/dτ
W(ξ ,η ,τ)
1
0
0 −0.5 −1
−0.5
−1.5 −2
−1 20
25
30
35
40
45
50
55
−2.5 −1.5
60
−1
−0.5
0
τ
1
1.5
12 RRMs 4 POMs
a/b=0.5,λ=50 1
Amplitude
0.5
W
0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
f/Hz Fig. 11. Comparison of a LCO for a panel of a/b ¼ 0.5 at λ ¼ 50 using 12 RRMs and 4 POMs: (a) time history; (b) phase portrait; (c) frequency spectrum.
panel deflections, including LCO deflections for a/b¼0.5 at λ ¼50, a/b¼1 at λ ¼200, a/b¼2 at λ ¼950 and chaotic deflection for a/b ¼2 at λ ¼1200. Compared with the POM shapes in Figs. 6–8, it can be seen clearly that the first dominant POM holds a great similarity with the panel deflection either for the LCO or the chaotic oscillation. This indicates that the POMs are the inherent modes of the fluttering panel. 3.2. Convergence studies The aeroelastic equations (AEs) using the RRMs or the POMs can be transformed into a smaller number of ODEs. Specifically, how much the AEs can be reduced will be explored with the convergence studies. A modal participation (MP) value that evaluates the contribution of the k-th RRMs or POMs to the total deflection is defined as MP ¼ maxjqk j=∑Kj¼ 1 maxjqj j. Tables 2 and 3 show the LCO deflections and MP values of the first six RRMs and POMs for the panels of a/b¼0.5 and a/b¼1 under three different dynamic pressures. The two methods have a good agreement on the LCO deflections. The MP
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
1.5
12 RRMs 4 POMs
a/b=1, λ=250
10
1
12 RRMs 4 POMs
a/b=1, λ=250
5 dW/d τ
0.5
W(ξ,η ,τ)
6201
0 −0.5
0
−5
−1 10
12
14
τ
16
18
−10 −1.5
20
−1
0
0.5
1
1.5
W
12 RRMs 4 POMs
a/b=1,λ=250 1
Amplitude
−0.5
0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
2.5
3
3.5
4
f/Hz Fig. 12. Comparison of a LCO for a panel of a/b ¼ 1 at λ¼ 250 using 12 RRMs and 4 POMs: (a) time history; (b) phase portrait; (c) frequency spectrum.
values are decreasing rapidly for both the RRMs and the POMs. Specifically, the first two POMs are the most dominant and they contribute more than 90 percent to the deflection, with the last two contributing less than 0.1 percent and 1 percent for a/b¼0.5 and a/b ¼1, respectively. Overall, the POMs converge faster than the RRMs. Based on the convergence studies previously discussed, Fig. 10(a) and (b) plots the LCO deflection versus dynamic pressure using POMs and RRMs for comparison. Obviously, only 4 POMs are necessary to obtain a good agreement with that using 12 ðM ¼ 6; N ¼ 2Þ RRMs. Specially, in Fig. 10(b) there is a departure between the POD and RR results from λ ¼285 to λ ¼300, which is due to the duality bifurcation depending on the initial conditions rather than the inaccurate discrepancy between the two methods. Ref. [12] provides the detail analysis of the duality bifurcation. Figs. 11 and 12 show the time histories, phase portraits and frequency spectra for a/b¼0.5 at λ ¼50 and a/b ¼1 at λ ¼250, respectively. It turns out that the POD method can give accurate LCO solutions with fewer modes than the RR method. 3.3. Efficiency studies In addition to the requirement of fewer modes, a computational cost saving is also our concern using the POD method. For a quantitative evaluation of the cost savings, the efficiency of the POD method is explored in this section. Three panels of a/b¼0.5, 1 and 2 are considered. Both the RR and POD methods are employed to determine the LCO responses. Table 4 provides the comparison between the two methods in terms of the Mode number, Integration time (It), Time step (Ts), and CPU time (Ct) (Intel Corei7-2720QM at 2.20 GHz) in RK4. For the panels of a/b¼0.5 and 1, the CPU time using the RR method with 12 RRMs is approximately 30 times that using the POD method with 4 POMs. For the panel of a/b¼2, the POD method with 6 POMs is approximately 60 times faster than the RR method with 16 RRMs. Hence, the computational cost is vastly reduced with fewer modes and a larger time step using the POD method. By contrast, when the POD method with a mode projection (from the POMs to RRMs) is used, there will be a huge difference. For the sake of brevity, the results are not reported here. However, the numerical simulations indicated that the projection-free POD method is approximately two orders of magnitude faster than the POD method with the mode projection. 3.4. Chaos studies Generally, a linear stability analysis of a nonlinear system is necessary. However, such studies for this nonlinear system are voluminous in the literature and the linear stability results for the presently studied plate configurations are well known.
6202
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Table 4 Efficiency comparison between the Rayleigh–Ritz and the POD methods. Method
a=b
RR POD RR POD RR POD
0.5
λ 50
1
200
2
950
Mode number
It
Ts
Ct (s)
12 4 12 4 16 6
100
0.005 0.01 0.005 0.01 0.002 0.005
715.1 20.8 375.0 11.8 1473.5 24.8
60 30
3 2
5 a/b=0.5,μ/Ma=0.05 I=R=12,J=S=2 M=6,N=2
Local amplitude extrema
Local amplitude extrema
4
1 0 −1 −2 −3 −4 −5 20
40
60
80
100
120
140
a/b=0.5,μ/Ma=0.05 I=R=14,J=S=2 K=4
0
−5 20
160
40
60
80
λ
100
120
140
160
λ
Fig. 13. Comparison of bifurcation diagram for a panel of a/b ¼ 0.5 using 12 RRMs and 4 POMs: (a) bifurcation diagram (12 RRMs); (b) bifurcation diagram (4 POMs).
4 3
3
RR POD
λ=64
RR POD
λ=64
2 1
1
dW/dτ
dW/dτ
2
0 −1
0 −1
−2 −2
−3 −4 −2
−1.5
−1
−0.5
0
0.5
1
1.5
2
−3 −0.1
−0.05
0
0.05
0.1
W
W 1
RR POD
λ=64
0.8 0.6 0.4 0.2 0
0
0.2
0.4
0.6
0.8
1
Fig. 14. Comparison of a quasi-periodic oscillation for the panel of a/b ¼ 0.5 at λ¼ 64 using 12 RRMs and 4 POMs: (a) phase portrait; (b) Poincare plot; (c) frequency spectrum.
The linear stability results agree with the present nonlinear results in the limit as the amplitude of the limit cycle oscillations approaches zero. Therefore, instead of the linear stability study, the bifurcation diagrams are presented to observe the chaotic motions.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
10
6
RR POD
λ=90
λ=90
4
5
2
dW/dτ
dW/dτ
6203
0
0 −2
−5 −4 −10 −4
−3
−2
−1
0
1
2
3
RR POD
−6 −0.25 −0.2 −0.15 −0.1 −0.05
4
0
0.05
0.1
0.15
W
W
λ=90
RR POD
0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
Fig. 15. Comparison of a chaotic oscillation for the panel of a/b ¼ 0.5 at λ¼ 90 using 12 RRMs and 4 POMs: (a) phase portrait; (b) Poincare plot; (c) frequency spectrum.
15
RR POD
λ=110
10
10
5
5
dW/dτ
dW/dτ
15
0 −5
RR S−POD
λ=110
0 −5
−10 −10 −15
−4
−3
−2
−1
0
1
2
3
4
−0.25 −0.2 −0.15 −0.1 −0.05
W
0
0.05 0.1 0.15
W 0.5
RR POD
λ=110
Amplitude
0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
f/Hz Fig. 16. Comparison of a chaotic oscillation for the panel of a/b ¼ 0.5 at λ ¼110 using 12 RRMs and 4 POMs: (a) phase portrait; (b) Poincare plot; (c) frequency spectrum.
6204
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
3
15
λ=132
2 1
5
0
dW/dτ
W(ξ,η,τ)
λ=132
10
−1
0
−2
−5
−3
−10
−4
440
450
460
470
480
490
−15 −5
500
−4
−3
−2
−1
τ
0
1
2
3
W 1
λ=132
Amplitude
0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
f/Hz Fig. 17. A quasi-periodic oscillation for the panel of a/b ¼ 0.5 at λ¼ 132 using 4 POMs: (a) time response; (b) phase portrait; (c) frequency spectrum.
5
λ=132
15
4
10
2
5
dW/dτ
W(ξ,η,τ)
3 1 0 −1
0 −5
−2 −3 −4
λ=132
−10 420
430
440
450
460
470
480
490
−4
−3
−2
−1
0
τ
2
3
4
λ=132
1
Normalized amplitude
1
W
0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
f/Hz Fig. 18. An almost chaotic oscillation for the panel of a/b ¼ 0.5 at λ ¼ 132 using 12 RRMs: (a) time response; (b) phase portrait; (c) frequency spectrum.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
3 a/b=1,μ/Ma=0.05 I=R=14,J=S=2 M=8,N=2
2
Local amplitude extrema
Local amplitude extreme
3
6205
1 0 −1 −2 −3 300
350
400
450
a/b=1,μ/Ma=0.05 I=R=14,J=S=2 K=4
2 1 0 −1 −2 −3 300
500
350
400
450
500
λ
λ
Fig. 19. Comparison of bifurcation diagram for a panel of a/b ¼ 1 using 16 RRMs and 4 POMs: (a) bifurcation diagram (16 RRMs); (b) bifurcation diagram (4 POMs).
2
a/b=2, μ/Ma=0.05 I=R=16,J=S=2 M=8,N=2
1.5 1
Local amplitude extrema
Local amplitude extrema
2
0.5 0 −0.5 −1 −1.5 −2 950
1000
1050
1100
1150
1.5 1 0.5 0 −0.5 −1 −1.5 −2 950
1200
a/b=2,μ/Ma=0.05 I=R=16,J=S=2 K=6
1000
1050
λ
1100
1150
1200
λ
Fig. 20. Comparison of bifurcation diagram for a panel of a/b ¼ 2 using 16 RRMs and 6 POMs: (a) bifurcation diagram (16 RRMs); (b) bifurcation diagram (6 POMs).
1.5
12 RRMs 4 POMs
a/b=0.5
μ/Ma=0.05
1
μ/Ma=0.25
0.8
μ/Ma=0.5
12 RRMs 4 POMs
a/b=1 μ/Ma=0.05 μ/Ma=0.25
1
μ/Ma=0.5
Wp
Wp
1.2
0.6 0.5
0.4 0.2 0
20
25
30
35
40
λ
45
50
55
60
0
100
150
200
250
300
λ
Fig. 21. Comparison of LCO deflection peak versus dynamic pressure at varying μ=Ma ¼ 0:05, 0.25, 0.5 using 12 RRMs and 4 POMs: (a) a panel of a/b ¼ 0.5; (b) a panel of a/b ¼ 1.
Fig. 13 shows the bifurcation diagrams for the panel of a/b¼0.5 using 4 POMs and 12 RRMs. They are in a good agreement with a symmetry property in terms of the amplitude of response. The reason for the symmetry lies in the fact that the aeroelastic system of this cantilever plate is symmetric, and it leads to the symmetric solutions depending on the initial conditions. Specifically, Figs. 14, 15 and 16 show the phase portraits, Poincare plots and frequency spectra for dynamic pressures of λ ¼64, 90 and 110 respectively. Note that the phase portraits using the POD method shown in Figs. 14(a), 15(a) and 16(a) are plotted with the symmetric solutions of the original ones for clear comparison with the RR results. Along with the frequency spectra shown in Figs. 14(c), 15(c) and 16(c), the two methods are in a very good agreement. The Poincare maps using the two methods are approximately symmetric in Figs. 14(b), 15(b) and 16(b). All these indicate that the POD method can be successfully used for the analysis of chaotic oscillations.
6206
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Exceptionally at λ ¼ 132, Fig. 17 shows a quasi-periodic motion with the POD method. By contrast, Fig. 18 indicates a nearly chaotic oscillation using the Rayleigh–Ritz method. For the panels of a/b¼1 and 2 shown respectively in Figs. 19 and 20, the bifurcation diagrams using the two methods have great similarities but with differences. The possible reasons are twofold. First, the chaotic oscillation is extremely sensitive to the dynamic system, so the approximate POD/ROM may differ from the original system, which would lead to distinct chaotic motions. Another possible reason is that the POMs extracted from the chaotic results are not unique. Maybe the POMs we have determined are not the optimal ones, which could affect the resulting chaotic motions. How to obtain the optimal POMs will be a subject for the future work.
3.5. Robustness of the POD modes for various μ=Ma Both λ and μ=Ma are the typical and important flight parameters. If the POMs can be applied to the variations of the flight parameters in a range of interest, the systems are easy to be analyzed without remodeling by the new POMs determined. The numerical simulations performed previously have shown that the POMs can be applied for different λ. Hence, a question arises here: whether a set of general POMs can be used for various μ=Ma? This is explored next. The general POMs are determined at a fixed μ=Ma ¼ 0:05 with the panel undergoing chaotic oscillations. Fig. 21(a) and (b) shows the LCO deflection versus dynamic pressure using the general POMs at μ=Ma ¼ 0:05, 0.25 and 0.5 for the panels of a/b¼0.5 and 1, respectively. They are all in an excellent agreement with the Rayleigh–Ritz results. Therefore, the question arisen previously can be answered as the general POMs determined at a fixed μ=Ma can be applied for different μ=Ma. Consequently, a conclusion may be drawn that as long as the POMs are determined from the chaotic responses, the POMs are robust for the same panel at varying flight parameters. The same panel is defined as a panel with the same dimensions and support conditions. Said another way, the POMs are the inherent characteristics of the fluttering cantilever panel independent of the flight conditions in this paper.
4. Concluding remarks A computationally fast POD method has been proposed and employed to study the aeroelastic response of a cantilever wing-like plate in supersonic flow. The presently proposed POD method has been verified to be about two orders of magnitude faster than the conventional POD method, due to the fact that the present method avoids the expensive projection procedure from the POD modes to the Rayleigh–Ritz modes by performing direct numerical derivatives to the POD modes. In addition to the computational efficiency improvement, the present projection-free POD method has also been compared with the Rayleigh–Ritz method for investigating the fluttering plate with three length–width ratios, namely a/b¼0.5, 1 and 2. Drawn from numerical simulations are the following conclusions: 1. For a good convergence, the required number of POD modes is much less than that of the Rayleigh–Ritz modes. Concretely, 4, 4 and 6 POD modes in contrast to 12, 12 and 16 Rayleigh–Ritz modes are required in analyzing the panels of a/b¼0.5, 1 and 2, respectively. Hence, the computational effort has been remarkably reduced using the POD method (1/ 60 1/30 that by the Rayleigh Ritz method). Also, numerical simulation indicates that although more POD modes should be included as a=b is increased for both methods, the POD method converges faster. 2. It was shown that the dominant (or first) POD mode extracted from the Rayleigh–Ritz chaotic results resembles the physical deflection of the panel, which accounts for the fact that fewer POD modes are required for capturing the panel deflection. 3. It has been demonstrated that the POD modes determined from the chaotic results for a specified set of parameters are universal modes, and are applicable to a variety of parametric variations (e.g. varying λ, and μ=Ma). This suggests that the POD modes are inherent characteristics of a plate with a given geometry and supporting pattern. 4. Previous studies believed that the POD method only works for capturing the periodic responses of nonlinear dynamical systems. In this study, the POD method has been shown to be able to obtain the periodic solutions as well as the chaotic solutions.
Acknowledgments The first author gratefully acknowledges the helpful discussions with the aeroelasticity research group in the Department of Mechanical Engineering and Materials Science, School of Engineering at Duke University, and colleagues of the College of Astronautics in Northwestern Polytechnical University.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
6207
Appendix A. Coefficients of C, D, A, B, Q and F
C ij ¼
Drs ¼
3 K K Z Z 2 K K Z Z h ∂φk1 ∂φk2 dui h h ∂φk1 ∂φk2 dui ∑∑ vj dx dy ν v dx dy ∑∑ a k1 k2 b a k1 k2 ∂x ∂x dx ∂y ∂y dx j 2 K K Z Z h h ∂φk1 ∂φk2 dvj ∑∑ u dx dy; ð1 νÞ b a k1 k2 ∂x ∂y i dy 3 K K Z Z 2 K K Z Z h ∂φk1 ∂φk2 dvs h h ∂φk1 ∂φk2 dvs ∑∑ ur dx dy ν ur dx dy ∑∑ b a b ∂y ∂y dy ∂x ∂x dy k1 k2 k1 k2 2 K K Z Z h h ∂φk1 ∂φk2 dur ∑∑ vs dx dy; ð1 νÞ a b k1 k2 ∂x ∂y dx Akk1 ¼
Bkk1 ¼
1 6
Z Z
1 6
Z Z
φk1 φk dx dy;
"
λ K ∑q 6β k1 k1
Z Z
∂φk1 Ma2 2 φk dx dy þ 2 ∂x Ma 1
rffiffiffiffi
μ K dqk1 ∑ λ k1 dt
(23)
(24)
a 4 Z Z ∂ 2 φ ∂ 2 φ
a 2 Z Z ∂2 φ ∂2 φ ∂2 φk1 ∂2 φk ∂2 φk1 ∂2 φk k1 k k1 k dx dy þ dx dyþ ν dx dy þ dx dy b b ∂x2 ∂x2 ∂y2 ∂y2 ∂y2 ∂x2 ∂x2 ∂y2
a 2 Z Z ∂2 φ ∂2 φ k1 k dx dy ; þ 2ð1 νÞ b ∂x∂y ∂x∂y Qk ¼
(22)
Z Z
(25)
#
φk1 φk dx dy ;
(26)
Z Z K K K dui ∂φk1 ∂φk ∂φk1 ∂φk2 ∂φk3 ∂φk vj dx dy þ ∑ ∑ ∑ qk1 qk2 qk3 dx dy h i j k1 dx ∂x ∂x ∂x ∂x ∂x ∂x k1 k2 k3 Z Z
a 3 a R S K dvs ∂φk1 ∂φk þ2 ∑ ∑ ∑ brs qk1 dx dy ur b h r s k1 dy ∂y ∂y Z Z Z Z
a 4 K K K
a 2 a I J K ∂φk1 ∂φk2 ∂φk3 ∂φk dui ∂φk1 ∂φk þ ∑ ∑ ∑ aij qk1 dx dy þ 2ν v dx dy ∑ ∑ ∑ qk1 qk2 qk3 b k1 k2 k3 b h i j k1 ∂y ∂y ∂y ∂y dx j ∂y ∂y Z Z
a 2 K K K ∂φk1 ∂φk2 ∂φk3 ∂φk dx dy ∑ ∑ ∑ qk1 qk2 qk3 þν b k1 k2 k3 ∂y ∂x ∂x ∂y Z Z Z Z
a a R S K
a 2 K K K dvs ∂φk1 ∂φk ∂φk1 ∂φk2 ∂φk3 ∂φk þ 2ν ∑ ∑ ∑ brs qk1 dx dyþ ν dx dy ∑ ∑ ∑ qk1 qk2 qk3 ur b h r s k1 b k1 k2 k3 dy ∂x ∂x ∂x ∂y ∂y ∂x ( Z Z
a 2 a I J K dvj ∂φk1 ∂φk ∂φk1 ∂φk þ ð1 νÞ þ dx dy ui ∑∑∑a q b h i j k1 ij k1 dy ∂y ∂x ∂x ∂y Z Z
a a R S K dur ∂φk1 ∂φk ∂φk1 ∂φk þ ∑ ∑ ∑ brs qk1 vs þ dx dy b h r s k1 dx ∂y ∂x ∂x ∂y ) Z Z
a 2 K K K ∂φk1 ∂φk2 ∂φk3 ∂φk ∂φk1 ∂φk2 ∂φk3 ∂φk þ dx dy : (27) ∑ ∑ ∑ qk1 qk2 qk3 þ b k1 k2 k3 ∂x ∂y ∂y ∂x ∂x ∂y ∂x ∂y
Fk ¼ 2
a
I
J
K
Z Z
∑ ∑ ∑ aij qk1
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
E.H. Dowell, Nonlinear oscillations of a fluttering plate, AIAA Journal 4 (1966) 1267–1275. E.H. Dowell, Nonlinear oscillations of a fluttering plate. II, AIAA Journal 5 (1967) 1856–1862. E.H. Dowell, Aeroelasticity of Plates and Shells, Noordhoff International Publishing, Groningen, Netherlands, 1975. C. Mei, A finite-element approach for nonlinear panel flutter, AIAA Journal 15 (1977) 1107–1110. W.L. Ye, E.H. Dowell, Limit cycle oscillation of a fluttering cantilever plate, AIAA Journal 29 (1991) 1929–1936. J. Zhou, Z.C. Yang, Y.S. Gu, Aeroelastic stability analysis of heated panel with aerodynamic loading on both surfaces, Science China: Technological Sciences 55 (2012) 2720–2726. D. Xie, M. Xu, E.H. Dowell, Proper orthogonal decomposition reduced-order model for nonlinear aeroelastic oscillations, AIAA Journal 52 (2014) 1–13. H.H. Dai, M. Schnoor, S.N. Atluri, A simple collocation scheme for obtaining the periodic solutions of the duffing equation, and its equivalence to the high dimensional harmonic balance method: subharmonic oscillations, Computer Modeling in Engineering and Sciences 84 (2012) 459. H.H. Dai, X.K. Yue, J.P. Yuan, S.N. Atluri, A time domain collocation method for studying the aeroelasticity of a two dimensional airfoil with a structural nonlinearity, Journal of Computational Physics 270 (2014) 214–237. H.H. Dai, J.K. Paik, S.N. Atluri, The global nonlinear Galerkin method for the analysis of elastic large deflections of plates under combined loads: a scalar homotopy method for the direct solution of nonlinear algebraic equations, Computers Materials and Continua 23 (2011) 69–99.
6208
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
[11] H.H. Dai, J.K. Paik, S.N. Atluri, The global nonlinear Galerkin method for the solution of von Karman nonlinear plate equations: an optimal & faster iterative method for the direct solution of nonlinear algebraic equations FðxÞ ¼ 0, using x ¼ λ½αF þ ð1 αÞBT B, Computers Materials and Continua 23 (2011) 155–185. [12] D. Xie, M. Xu, H. Dai, E.H. Dowell, Observation and evolution of chaos for a cantilever plate in supersonic flow, Journal of Fluids and Structures (2014), http://dx.doi.org/10.1016/j.jfluidstructs.2014.05.015, in press. [13] I.R. Dixon, C. Mei, Finite element analysis of large-amplitude panel flutter of thin laminates, AIAA Journal 31 (1993) 701–707. [14] C.E. Gray, C. Mei, Large-amplitude finite element flutter analysis of composite panels in hypersonic flow, AIAA Journal 31 (1993) 1090–1099. [15] P.J. Attar, E.H. Dowell, A reduced order system id approach to the modelling of nonlinear structural behavior in aeroelasticity, Journal of Fluids and Structures 21 (2005) 531–542. [16] X. Guo, C. Mei, Using aeroelastic modes for nonlinear panel flutter at arbitrary supersonic yawed angle, AIAA Journal 41 (2003) 272–279. [17] X. Guo, C. Mei, Application of aeroelastic modes on nonlinear supersonic panel flutter at elevated temperatures, Computers & Structures 84 (2006) 1619–1628. [18] B.I. Epureanu, L.S. Tang, M.P. Paidoussis, Coherent structures and their influence on the dynamics of aeroelastic panels, International Journal of NonLinear Mechanics 39 (2004) 977–991. [19] B.I. Epureanu, L.S. Tang, M.P. Paidoussis, Exploiting chaotic dynamics for detecting parametric variations in aeroelastic systems, AIAA Journal 42 (2004) 728–735. [20] M. Amabili, A. Sarkar, M.P. Paidoussis, Reduced-order models for nonlinear vibrations of cylindrical shells via the proper orthogonal decomposition method, Journal of Fluids and Structures 18 (2003) 227–250. [21] M. Amabili, A. Sarkar, M.P. Paidoussis, Chaotic vibrations of circular cylindrical shells: Galerkin versus reduced-order models via the proper orthogonal decomposition method, Journal of Sound and Vibration 290 (2006) 736–762. [22] D. Xie, M. Xu, A simple proper orthogonal decomposition method for von Karman plate undergoing supersonic flow, Computer Modeling in Engineering and Sciences 93 (2013) 377–409. [23] E.H. Dowell, On asymptotic approximations to beam modal functions, Journal of Applied Mechanics 51 (1984) 439. [24] A. Sarkar, M.P. Païdoussis, A cantilever conveying fluid: coherent modes versus beam modes, International Journal of Non-Linear Mechanics 39 (2004) 467–481.
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
Projection-free proper orthogonal decomposition method for a cantilever plate in supersonic flow Dan Xie a,n, Min Xu a, Earl H. Dowell b a b
College of Astronautics, Northwestern Polytechnical University, Xi'an 710072, China Duke University, Durham, NC 27708-0300, United States
a r t i c l e in f o
abstract
Article history: Received 5 May 2014 Received in revised form 25 June 2014 Accepted 25 June 2014 Handling Editor: L.G. Tham Available online 19 July 2014
This paper explores the use of the proper orthogonal decomposition (POD) method for supersonic nonlinear flutter of a cantilever plate or wing. The aeroelastic equations are constructed using von Karman plate theory and first-order piston theory. The twodimensional POD modes (POMs) in xy plane are determined from the chaotic results given by the traditional Rayleigh–Ritz (RR) approach. For a specific structure, the POMs need to be calculated once and then can be used for various parameters of interest. The derivatives of the POMs are calculated numerically to avoid the complex projection from the POMs to the Rayleigh–Ritz modes (RRMs). Numerical examples demonstrate that the POD method using 4 POMs can obtain accurate limit cycle oscillation (LCO) results with substantial computational cost savings, compared with 12 RRMs by the Rayleigh–Ritz method. The POD method is employed for the analysis of the chaotic oscillations. It is also demonstrated that the POD modes are robust over a range of flight parameters. & 2014 Elsevier Ltd. All rights reserved.
1. Introduction The aeroelastic system of a fluttering panel undergoing supersonic flow, using von Karman theory accounting for the structural nonlinearity and piston theory accounting for the aerodynamics, is a classical self-excited autonomous system, and has been intensively investigated by many researchers [1–7]. The mathematical model for this problem is a set of partial differential equations (PDEs), with both spatial and temporal variables. Analytical solutions to this set of PDEs are not available due to the existence of the structural nonlinearity which nullifies the additivity pertaining to the linear systems [8,9]. Nevertheless, considering the simple spatial domain (due to the simple rectangular geometry of the plate), several semi-analytical methods can be applied to solve this problem. For a simply supported plate [1,10,11], the Galerkin method has been used in the spatial domain to transform the governing PDEs into a set of coupled ordinary differential equations (ODEs) in time. Subsequently, the resulting ODEs were integrated by the numerical integration method. In addition, the Rayleigh–Ritz method was applied by Ye and Dowell [5] to analyze the LCO of a cantilever plate. For the same model, the chaotic oscillations are explored using the Rayleigh–Ritz method by Xie et al. [12]. However, when the plate has a complex geometry, the semi-analytical methods (the Galerkin method and the Rayleigh–Ritz method) become invalid. The finite element method (FEM) is more flexible for handling complex geometries. It has been employed extensively to study the nonlinear panel flutter [4,13,14]. However, the number of degrees-of-freedom (DOFs) required in the finite element analysis
n
Corresponding author. Tel.: þ86 18291964199. E-mail address: [email protected] (D. Xie).
http://dx.doi.org/10.1016/j.jsv.2014.06.039 0022-460X/& 2014 Elsevier Ltd. All rights reserved.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Nomenclature
R, S
a, b aij ; brs D, E h I, J
r, s t u, v uiðrÞ ; vjðsÞ V, vj W, w x, y, z
i, j K k Ma M, N m, n Δp Qk , Q q, q qk
plate length, plate width mode coordinate for u, v plate stiffness, Young's modulus plate thickness total mode number retained in x; y direction for u mode number in x; y direction for u total mode number retained for w (POD) POD mode number Mach number total mode number retained in x; y direction for w (RR) mode number in x; y direction for w aerodynamic pressure, positive in direction opposite to w generalized aerodynamic force, snapshots matrix dynamic pressure, POD mode coordinate in vector POD mode coordinate for w
β λ λpk μ
ν ξ; η ρ; ρm τ Φ Ψ, φ k φ~ k , φ k
6191
total mode number retained in x; y direction for v mode number in x; y direction for v time in-plane displacement in length and width mode in x; y direction for u(v) POD eigenvectors matrix, POD eigenvector w/h, panel transverse deflection streamwise, spanwise, normal coordinates ðMa2 1Þ1=2 , compressibility correction factor 2qa3 =βD, nondimensional dynamic pressure POD eigenvalue ρa=ρm h, nondimensional fluid/structure mass ratio Poisson ratio x=a; y=b air density, plate density 4 tðD=ρm ha Þ1=2 , nondimensional time correlation matrix POD modes matrix, POD mode in xy plane POD mode in vector form, POD mode in matrix form
is very large. Therefore, the semi-analytical methods are employed for the cantilever plate in this paper. Of course, one can use finite element analysis to find the eigenmodes of the structure and this is often done in aeroelastic analysis. All these traditional approaches use the assumed modes. Consequently, the coupled nonlinear modal equations with large number of modes are computationally costly for flutter analysis. For solving this problem, the reduced-order model (ROM) is proposed by many researchers. For the systems with original full-order equations of motion unknown, Attar [15] constructed a ROM using a system identification technique to investigate the LCO of a 451 delta wing. Guo and Mei [16,17] first used the aeroelastic modes for nonlinear panel flutter to reduce the number of ODEs and save computational costs greatly. Similarly, a set of POMs for the aeroelastic panels has been used to construct a ROM [18]. The POD method was also used to identify the coherent structures of the panel dynamics [19]. Furthermore, Amabili et al. [20,21] applied the POD method with a Galerkin projection (projection from the POMs to the Galerkin modes) for a water-filled circular cylindrical shell. More recently, a POD reduced-order model (POD/ROM) with a Galerkin projection has been constructed by Xie et al. [7] to solve the supersonic nonlinear oscillations of a simply supported plate with remarkable computational cost savings. However, seldom has a reduced-order model been used for a fluttering cantilever plate. In the present study, we create a reduced-order model using a computationally fast POD method to investigate the nonlinear flutter of a cantilever plate undergoing supersonic flow. The newly proposed POD method is different from the conventional POD method because it does not require the complex and computationally expensive projection from the POMs to the RRMs. The projection-free POD method has been applied for LCO analysis of a two-dimensional panel with all edges simply supported [22]. In this paper, LCO and chaotic oscillations are both investigated with the POD method. Section 2 provides the details for the construction of the POD/ROM. The convergence and efficiency studies are performed and the numerical examples are presented in Section 3. Moreover, a set of general POD modes is used for a range of physical parameters. Finally, some conclusions are drawn in Section 4.
2. Theoretical analysis A wing-like cantilever plate shown in Fig. 1 is constrained along the root edge with supersonic flow passing over the upper surface parallel to the fixed edge. The plate is taken as homogeneous isotropic aluminum with the material constants E ¼7.10 GPa, ρ ¼ 2:8 103 kg=m3 and ν ¼0.33. The von Karman plate theory and the first-order piston aerodynamic theory are used to account for the structural nonlinearity and the aerodynamic force, respectively. The modal equations being studied are derived directly from Hamilton's Principle and Lagrange's equations and thus the von Karman plate equations per se are never explicitly use. Of course the von Karman equations can be derived from Hamilton's Principle. The POD method using the beam modes for the in-plane displacements and the POD modes for the transverse displacement is employed to discretize the governing equations.
6192
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Y
Y
h
a
Flow
b
Z
X Top view
Side view
Fig. 1. Geometry of a cantilever panel or wing in top view and side view.
2.1. Numerical POD modes in xy plane It is the authors' experience in the present study as well as earlier studies [7] that using a single chaotic time history with broad band frequency response to construct the snapshots for determining the POD modes is much better than using a time history for a single frequency response. This is because the frequency spectrum of a single chaotic time series is broad and excites many POD or eigenmodes of the system. One could use many time histories at single frequencies for constructing a collection of snapshots as well, but a single chaotic time history is usually more computationally efficient. In this paper, the snapshots are gathered from the chaotic results by the Rayleigh–Ritz method. A similar procedure for the POMs' determination as previous work [7] is performed. First, the chaotic time response of the transverse deflection is assembled in a matrix as snapshots. Different columns correspond to the response at different times τ. Note that every column consists of WðτÞ of all the nodes in both x and y directions. We sort the order by fixing an individual η with ξ changed from ξ1 to ξN and then, fixing a subsequent η with ξ changed in the same way. This is different from the POMs in one-dimension [7]. Here the snapshots matrix takes the form as follows: 0
Wðξ1 ; η1 ; t 1 Þ B ⋮ B B B WðξN ; η1 ; t 1 Þ B B B Wðξ1 ; η2 ; t 1 Þ B B ⋮ B Q ¼B B WðξN ; η2 ; t 1 Þ B B ⋮ B B B Wðξ1 ; ηN ; t 1 Þ B B ⋮ @ WðξN ; ηN ; t 1 Þ
Wðξ1 ; η1 ; t 2 Þ
⋯
⋮
⋱
WðξN ; η1 ; t 2 Þ
⋯
Wðξ1 ; η1 ; t 2 Þ ⋮
⋯ ⋱
WðξN ; η1 ; t 2 Þ
⋯
⋮
⋱
Wðξ1 ; ηN ; t 2 Þ
⋯
⋮
⋱
WðξN ; ηN ; t 2 Þ
⋯
Wðξ1 ; η1 ; t J Þ
1
C C C WðξN ; η1 ; t J Þ C C C Wðξ1 ; η2 ; t J Þ C C C ⋮ C C; WðξN ; η2 ; t J Þ C C C ⋮ C C Wðξ1 ; ηN ; t J Þ C C C ⋮ A WðξN ; ηN ; t J Þ ⋮
(1)
where Wðξi ; ηj ; τk Þ denotes the non-dimensional transverse deflection of i-th point in the x direction and j-th point in the y direction at k-th time point ði ¼ 1; 2; …; N; j ¼ 1; 2; …; N; k ¼ 1; 2; …; JÞ. Note that J and N are the number of time points and nodes in the x or y direction, respectively. Then, we can get the correlation matrix: T
Φ¼Q Q:
(2)
Φvk ¼ λpk vk ;
(3)
And an eigenvalue problem is yielded:
where λk ; vk ðk ¼ 1; 2; …; JÞ are the eigenvalues and eigenvectors of Φ, respectively. Sorting the eigenvalues in a decreasing order, the eigenvectors vk are sorted correspondingly: p
λp1 ≧λp2 ≧⋯≧λpJ :
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
6193
A commonly used criterion for choosing K, the optimal number of basis retained in the POD subspace, is the so-called energy contribution criterion. For a predefined “energy percentage”, here ϵ ¼0.9999, K is chosen by ∑Kk ¼ 1 λk
p
F ðkÞ ¼
p ≧0:9999; ∑Jj ¼ 1 j
λ
and, then, the modal vectors of the correlation matrix is V ¼ ½v1 ; v2 ; …; vK . With the matrix V, the snapshots are linearly combined to generate a smaller number of basis vectors assembled in a matrix: 0 1 φ1 ðξ1 ; η1 Þ φ2 ðξ1 ; η1 Þ ⋯ φK ðξ1 ; η1 Þ B C ⋮ ⋮ ⋱ ⋮ B C B C B φ1 ðξN ; η1 Þ φ2 ðξN ; η1 Þ ⋯ φK ðξN ; η1 Þ C B C B C B φ1 ðξ1 ; η2 Þ φ2 ðξ1 ; η2 Þ ⋯ φK ðξ1 ; η2 Þ C B C B C ⋮ ⋮ ⋱ ⋮ C Ψ¼QV¼B (4) B C; B φ1 ðξN ; η2 Þ φ2 ðξN ; η2 Þ ⋯ φK ðξN ; η2 Þ C B C B C ⋮ ⋮ ⋱ ⋮ B C B C B φ1 ðξ1 ; ηN Þ φ2 ðξ1 ; ηN Þ ⋯ φK ðξ1 ; ηN Þ C B C B C ⋮ ⋮ ⋱ ⋮ @ A φ1 ðξN ; ηN Þ φ2 ðξN ; ηN Þ ⋯ φK ðξN ; ηN Þ ~ k ðk ¼ 1; 2; …; KÞ. where the k-th column of matrix (4) is corresponding to the POD mode vector φ 2.2. Aeroelastic equations with the POD modes Based on the mode functions for a free-free beam and a cantilever beam [23], we assume the mode functions for the inplane displacements of u and v: ui ξ ¼ cos iπξ;
2j 1 vj η ¼ sin πη; 2
ur ξ ¼ cos r πξ;
2s 1 vs η ¼ sin πη: 2
(5)
For the transverse deflection, the POD mode φk ðξ; ηÞ in xy plane is used. Using the panel thickness, a set of appropriate nondimensionalizations is employed: u ¼ u=h; v ¼ v=h and w ¼ w=h. Then the streamwise, spanwise and transverse displacements in the nondimensional form can be expressed in terms of the mode functions as follows: J
I
u ¼ ∑ ∑ aij ðτÞui ðξÞvj ðηÞ; i¼1j¼1 R
S
v ¼ ∑ ∑ brs ðτÞur ðξÞvs ðηÞ; r ¼1s¼1 K
w ¼ ∑ qk ðτÞφk ðξ; ηÞ:
(6)
k¼1
Based on the von Karman nonlinear strain–displacement relation, the normal strains ϵx ; ϵy in x; y directions and the shear strain γxy in xy plane take the forms as ∂u 1 ∂w 2 ∂v 1 ∂w 2 ∂v ∂u ∂w ∂w : (7) ϵx ¼ þ ; ϵy ¼ þ ; γ xy ¼ þ þ ∂x 2 ∂x ∂y 2 ∂y ∂x ∂y ∂x ∂y Then the elastic energies US and UB induced by the stretching and bending are Z Z Eh 1ν 2 US ¼ ϵ2x þ ϵ2y þ2νϵx ϵy þ γ xy dx dy; 2 2 2ð1 ν Þ UB ¼
D 2
" 2 2 #) 2 Z Z ( 2 ∂ w ∂2 w ∂2 w ∂2 w ∂ w þ 2 ð 1 ν Þ dx dy; ∂x∂y ∂x2 ∂y2 ∂x2 ∂y2
(8)
and T is the kinetic energy: T¼
1 2
Z Z
ρm h
∂w 2 dx dy: ∂t
(9)
6194
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Then L is the Lagrangian taking the form L ¼ T ðU S þ U B Þ, and Lagrange's equations are 8 ∂L > > ¼ 0; > > > ∂aij > > > < ∂L ¼ 0; ∂brs > > > > > d ∂L ∂L > > > : dt ∂ðdq =dtÞ ∂q þ Q k ¼ 0: k k Substituting L ¼ T ðU S þ U B Þ into Eq. (10) yields 8 ∂U S > > ¼ 0; > > > ∂aij > > > < ∂U S ¼ 0; ∂brs > > > h i ∂U > > d ∂U S B > ∂T > > : dt ∂ðdqk =dtÞ þ ∂q þ ∂q þQ k ¼ 0: k k Aerodynamic force is generalized as
Z Z Qk ¼
Δpðx; y; t Þ
∂w dx dy; ∂qk
where the aerodynamic pressure Δp is modeled using first-order piston or quasi-steady aerodynamic theory ! " # 2q ∂w Ma2 2 1 ∂w þ Δp ¼ : β ∂x Ma2 1 V 1 ∂t
(10)
(11)
(12)
(13)
Performing a subsequent substitution of equations as Eqs. (5) (- 6) (- 7) - (8) - (9), along with Eq. (13) - (12), and finally substituting them into Eq. (11) with appropriate nondimensionalization, the following governing equations are yielded: Ca a þCb b ¼ C; Da a þDb b ¼ D;
(14)
2
A
d q dt
2
þ Bqþ F þ Q ¼ 0:
(15)
Eq. (14) is a set of algebraic equations which can be solved for aij and brs in terms of qk, and Eq. (15) is a set of ordinary 2 differential equations in terms of aij ; brs ; qk ; dqk =dt and d qk =dt 2 . Pre-multiplying Eq. (15) by the inverse of matrix A yields 2
d q dt
2
þ A 1 Bq þ A 1 F þ A 1 Q ¼ 0:
(16)
This set of ODEs can be solved numerically by the fourth-order Runge–Kutta method (RK4). From Eq. (6), the panel displacements in time can be calculated. The matrices Ca ; Cb ; Da and Db are all the same as Ref. [5], hence they are omitted here for avoidance of repetition. However, using the numerical derivatives of the POMs, the other coefficients of C; D; A; B; Q and F take different expressions. Based on the mode expansion for the transverse deflection in Eq. (6) compared to Eq. (7) in Ref. [5], we can replace the terms in the expressions of C; D; A; B; Q and F in Ref. [5] as follows: Z Z Z Z M N M N K K ∑ ∑ ∑ ∑ qmn qpl ϕm ϕp dx ψ n ψ l dy- ∑ ∑ qk1 qk2 φk1 φk2 dx dy (17) m n
p
l
k1 k2
And then we can obtain the expressions of C; D; A; B; Q and F given in Appendix A for the present study. 2.3. Numerical derivatives of the POD modes Ye and Dowell [5] used the analytical beam modes for the transverse deflection. The derivatives of the beam modes can be calculated analytically and exactly. By contrast, the numerical POMs are used for the transverse deflection in this study. Hence, the differential manipulation on the POMs is performed numerically. Consequently, the differential and integral operations on the numerical POMs in the transformation of PDEs into ODEs may propagate the numerical errors. To circumvent this problem, a POD method projecting the POMs onto the assumed modes (RRMs or Galerkin modes) semi-analytically is employed [7,24]. In that way, any linear or nonlinear manipulation can be performed analytically to avoid any spatial errors. However, the authors’ experience has shown that the mode projection will lead to a complex math derivation of the governing equations and subsequently increase the computational costs vastly for the cantilever plate.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
6195
Table 1 Chaotic responses for the POMs determination. a=b
λ
Number of RRMs
It
Ts
Tc
Interval
0.5 1 2
90 500 1200
M ¼ 8, N ¼2, I ¼ R ¼ 14 M ¼ 8, N ¼2, I ¼ R ¼ 14 M ¼ 8, N ¼2, I ¼ R ¼ 16
t ¼ 500 t ¼ 200 t ¼ 200
0.005 0.005 0.002
0.4 0.3 0.3
5 2 5
Hence, the derivatives of the POMs are directly calculated by a three-point difference technique to avoid the mode projection in this paper. ~ k of k-th column in Eq. (4) into a matrix form, by sorting φk ðξ; ηÞ with different η First, we rewrite the POD mode vector φ into different columns, and different ξ into different rows: 0 1 φk ðξ1 ; η1 Þ φk ðξ1 ; η2 Þ ⋯ φk ðξ1 ; ηN Þ B C B φk ðξ2 ; η1 Þ φk ðξ2 ; η2 Þ ⋯ φk ðξ2 ; ηN Þ C C: φk ¼ B (18) B C ⋮ ⋮ ⋱ ⋮ @ A φk ðξN ; η1 Þ φk ðξN ; η2 Þ ⋯ φk ðξN ; ηN Þ And then the derivatives with respect to ξ and η can be calculated numerically by a three-point difference technique. The first and second derivatives with respect to ξ are performed in every column of matrix (18) as follows: ∂φk ðξ; ηÞ 3φk ðξi ; ηÞ þ 4φk ðξi þ 1 ; ηÞ φk ðξi þ 2 ; ηÞ ¼ ∂ξi 2dξ ∂φk ðξ; ηÞ φk ðξi þ 2 ; ηÞ φk ðξi ; ηÞ ¼ ∂ξi þ 1 2dξ
∂φk ðξ; ηÞ φk ðξi ; ηÞ 4φk ðξi þ 1 ; ηÞ þ3φk ðξi þ 2 ; ηÞ ¼ ∂ξi þ 2 2dξ ∂2 φk ðξ; ηÞ φk ðξi ; ηÞ 2φk ðξi þ 1 ; ηÞ þ φk ðξi þ 2 ; ηÞ ¼ ∂ξi ðdξÞ2 ∂2 φk ðξ; ηÞ φk ðξi ; ηÞ 2φk ðξi þ 1 ; ηÞ þ φk ðξi þ 2 ; ηÞ ¼ ∂ξi þ 1 ðdξÞ2
∂2 φk ðξ; ηÞ φk ðξi ; ηÞ 2φk ðξi þ 1 ; ηÞ þ φk ðξi þ 2 ; ηÞ ¼ ; ∂ξi þ 2 ðdξÞ2
i ¼ 1; 2; …; N 1:
(19)
In the same way, the derivatives with respect to η are performed in every line of matrix (18). Finally, substitute the numerical derivatives of the POMs into Eqs. (22), (23), (24), (25), (26) and (27) to calculate the coefficients of C; D; A; B; Q and F respectively. 2.4. Numerical integration for the aeroelastic equations Because the PDEs are transformed into low-order ODEs with the numerical POMs, the differential and integral operations on the numerical POMs result in a set of numerical ODEs. Hence, the integrations in the ODEs have to be solved numerically. A one-dimensional trapezoidal integration formula has been used extensively and takes the form " # Z b N1 h f ðxÞ dx ¼ f ðaÞ þ2 ∑ f ðxi Þ þ f ðbÞ ; 2 a i¼1 h¼
ba ; xi ¼ a þ i h; i ¼ 1; 2; …; N 1: N
(20)
In the present study, due to the two-dimensional POMs in xy plane, a corresponding two-dimensional trapezoidal integration formula is proposed and derived as follows: # Z bZ a Z b " N1 h f ð0; yÞ þ 2 ∑ f ðxi ; yÞ þ f ða; yÞ dy f ðx; yÞ dx dy ¼ 0 0 0 2 i¼1 # 2 (" N 1 h ¼ f ð0; 0Þ þ 2 ∑ f ðxi ; 0Þ þ f ða; 0Þ 2 i¼1 " # N1
þ2 ∑
j¼1
"
N1
f ð0; yj Þ þ 2 ∑ f ðxi ; yj Þ þ f ða; yj Þ i¼1
N1
#)
þ f ð0; bÞ þ 2 ∑ f ðxi ; bÞ þ f ða; bÞ i¼1
6196
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
0
4
λ=90
λ=90
3
−1 −2
1
dW/dτ
W(ξ,η,τ)
2
0 −1
−3 −4
−2
−5
−3 −4
0
100
200
300
400
−6 −0.25 −0.2 −0.15 −0.1 −0.05
500
τ
0
0.05
0.1
0.15
W
Fig. 2. Time history and Poincare plot for a panel of a/b ¼0.5 at λ ¼90: (a) time history; (b) Poincare map.
3
10
λ=500
λ=500 5
2
0
dW/dτ
W(ξ,η,τ)
1 0 −1
−10 −15
−2 −3
−5
−20
0
50
100
150
−25
200
−0.5
0
τ
0.5
W
Fig. 3. Time history and Poincare plot for a panel of a/b¼ 1 at λ¼ 500: (a) time history; (b) Poincare map.
60
2
λ=1200
1.5
λ=1200
40
dW/dτ
W(ξ,η,τ)
1 0.5 0 −0.5 −1
20 0 −20
−1.5 −2
0
50
100
150
200
−40
τ
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
W
Fig. 4. Time history and Poincare plot for a panel of a/b ¼2 at λ ¼ 1200: (a) time history; (f) Poincare map.
¼
) 2 ( N1 N1 h f ð0; 0Þ þ f ða; 0Þ þ f ð0; bÞ þ f ða; bÞ þ 4 ∑ ∑ f ðxi ; yj Þ 2 i¼1 j¼1
i N1 N 1h
þ 2 ∑ f ðxi ; 0Þ þ f ðxi ; bÞ þ 2 ∑ f 0; yj þ f a; yj ; i¼1
h¼
j¼1
a b or ; xi ¼ i h; yj ¼ j h; i; j ¼ 1; 2; …; N 1: N N
Here we let N be 100 in the present paper.
(21)
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
6197
Accumulative energy contribution %
100 99 98 97 96 a/b=0.5 a/b=1 a/b=2
95 94
0
2
4 6 POD mode number
8
10
Fig. 5. Accumulative energy of POD eigenvalues versus POD mode number for panels of a/b¼ 0.5, a/b¼ 1 and a/b¼ 2.
0.02
0.05
0.01
0
0 1 0.5 ξ
0 0
0.5 η
1
−0.05 1 0.5 ξ
0.05
0.05
0
0
−0.05 1
−0.05 1
0.5 ξ
0 0
0.5 η
1
0 0
0.5 ξ
0 0
0.5 η
0.5
1
1
η
a/b=0.5, λ =90 Fig. 6. The first 4 POD modes for a panel of a/b ¼ 0.5 determined at λ¼ 90 using 16 RRMs.
3. Results and discussions In this section, the nonlinear dynamics of a cantilever plate in supersonic flow shown in Fig. 1 are investigated, using the POD method free of the mode projection from the POMs to RRMs. The results are compared with those using the Rayleigh– Ritz method. For all the calculations here, h ¼0.01 m, Ma¼2, μ ¼ 0.1 and J¼ S¼2 unless otherwise noted. All the figures are plots with respect to a typical point ξ ¼0.75, η ¼1.0. The initial conditions are q1ðτ ¼ 0Þ ¼ 0:1 with others being zero. Here q1 is q11 in Eq. (6) and deflection Wp is the peak of w. Three panels of a/b¼0.5, 1 and 2 are considered and the POMs are determined from the Rayleigh–Ritz chaotic solutions. The property of the POMs is explored by comparing the POM shapes with the panel deflections. The convergence study is performed based on the modal participation of the POMs, and the efficiency is evaluated in terms of the computational costs. Finally, the chaotic oscillations are analyzed and the robustness of the POD modes to the system parameter variations is investigated. 3.1. Property of the POD modes The results in Ref. [7] suggest that the chaotic oscillations are the best choice to determine the POD modes, which contain the richest information of the dynamic systems. Thus, the general POD modes only need to be calculated once and can be
6198
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
0.02
0.05
0.01
0
0 1
−0.05 1
0.5 ξ
0 0
1
0.5
0.5
0.05
0.05
0
0
−0.05 1
−0.05 1
0.5 ξ
0 0
0 0
ξ
η
1
0.5
0.5 η
0.5 0 0
ξ
η
0.5
1
1
η
a/b=1, λ =500 Fig. 7. The first 4 POD modes for a panel of a/b ¼1 determined at λ ¼500 using 16 RRMs.
0.04
0.05
0.05
0.02
0
0
−0.05 1
−0.05 1
0 1
1
0.5 ξ
0 0
1
0.5
0.5 η
0 0
ξ
ξ
0.05
0.05
0.05
0
0
0
−0.05 1
1
0.5 ξ
0 0
0.5 η
−0.05 1
1
0.5 ξ
0 0
0.5 η
1
0.5
0.5 η
0 0
0.5 η
0 0
0.5 η
−0.05 1
1
0.5 ξ
a/b=2, λ=1200 Fig. 8. The first 6 POD modes for a panel of a/b ¼ 2 determined at λ ¼1200 using 16 RRMs.
2
2
1
1
0 1
(w/h)p
(w/h)p
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
0.5
0 0
0.5 η
0 1
1
1
0.5
0.5
0 1
0.5 ξ
0 0
0.5 ξ
(w/h)p
(w/h)p
ξ
1
1
0.5
0 1
0.5 ξ
η
0 0
0 0
6199
0.5
1
η
0.5 η
1
Fig. 9. Deflection shapes for panels of a/b¼ 0.5, 1, 2; λ¼ 50, 200, 950 are corresponding to LCO deflection shapes and λ ¼ 1200 is about chaotic deflection shape.
Table 2 Modal participation for a/b ¼ 0.5. λ
40 50 55
Mode
RRM POM RRM POM RRM POM
ðw=hÞp
MP (%) q1
q2
q3
q4
q5
q6
0.7825 0.7508 1.0632 1.0253 1.1682 1.1509
80.5228 85.4608 82.1996 89.0744 82.8710 90.8220
8.6812 8.3048 9.6046 5.9960 9.8578 5.1710
8.2011 6.0025 6.2903 4.7139 5.5860 3.7881
1.6080 0.1591 1.2246 0.1494 1.0922 0.1528
0.5959 0.0334 0.4136 0.0356 0.3581 0.0341
0.3910 0.0395 0.2673 0.0307 0.2348 0.0321
ðw=hÞp
MP (%) q1
q2
q3
q4
q5
q6
75.4875 79.8955 75.8395 86.3509 75.6475 86.7527
1.4745 8.9275 8.4764 6.1160 9.1233 6.1182
18.7835 7.5841 11.1256 5.3574 10.7588 5.1964
2.5438 3.1181 2.3365 1.6009 1.9005 1.4971
1.1994 0.3939 1.7402 0.5033 2.0711 0.3619
0.5112 0.0809 0.4819 0.0715 0.4988 0.0737
Table 3 Modal participation for a/b ¼ 1. λ
100 200 250
Mode
RRM POM RRM POM RRM POM
0.3053 0.2963 1.0392 1.0344 1.1689 1.1631
used for the various cases. Specifically, λ ¼90, 500 and 1200 are selected for the panels of a/b¼0.5, 1, and 2, respectively to generate the chaotic responses by the Rayleigh–Ritz method [12]. Herein we choose the parameters as shown in Table 1, where “It” (Integration time of RK4), “Ts” (Time step of RK4), “Tc” (Time chosen for POMs) and “Interval” are all the related time parameters in RK4 using the Rayleigh–Ritz method. a/b¼1, λ ¼500 is taken as an example, for which the time marching (RK4) is carried out from τ ¼ 0 to τ ¼200. Only the last 200 0.3 ¼60 is chosen for the POMs determination, and one solution is selected as a snapshot every 0.005 2 ¼0.01. Figs. 2–4 show the time responses and Poincare plots for the chaotic oscillations aforementioned. Based on Section 2.1, Fig. 5 shows the accumulative energy of the POD eigenvalues. Obviously, 4, 4 and 6 POMs are required for the panels of a/b¼0.5, 1 and 2, respectively. Subsequently, Figs. 6, 7 and 8 show the first four POM shapes for the panels of a/b¼0.5 and a/b¼1, and the first six POMs for the panel of a/b¼ 2, respectively. To demonstrate the POM property, Fig. 9 presents the
6200
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
1.4
12 RRMs 4 POMs
1
1.2
0.8
Wp
Wp
a/b=1
1
0.8 0.6
0.6
0.4
0.4
0.2 0
12 RRMs 4 POMs
a/b=0.5
0.2 20
25
30
35
40
45
50
55
0
60
50
100
150
200
250
300
λ
λ
Fig. 10. Comparison of LCO deflection peak versus dynamic pressure using 12 RRMs and 4 POMs: (a) a panel of a/b ¼0.5; (b) a panel of a/b ¼ 1.
1.5
2.5
12 RRMs 4 POMs
a/b=0.5, λ=50
2
1
1.5
0.5
0.5
12 RRMs 4 POMs
a/b=0.5, λ=50
dW/dτ
W(ξ ,η ,τ)
1
0
0 −0.5 −1
−0.5
−1.5 −2
−1 20
25
30
35
40
45
50
55
−2.5 −1.5
60
−1
−0.5
0
τ
1
1.5
12 RRMs 4 POMs
a/b=0.5,λ=50 1
Amplitude
0.5
W
0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
f/Hz Fig. 11. Comparison of a LCO for a panel of a/b ¼ 0.5 at λ ¼ 50 using 12 RRMs and 4 POMs: (a) time history; (b) phase portrait; (c) frequency spectrum.
panel deflections, including LCO deflections for a/b¼0.5 at λ ¼50, a/b¼1 at λ ¼200, a/b¼2 at λ ¼950 and chaotic deflection for a/b ¼2 at λ ¼1200. Compared with the POM shapes in Figs. 6–8, it can be seen clearly that the first dominant POM holds a great similarity with the panel deflection either for the LCO or the chaotic oscillation. This indicates that the POMs are the inherent modes of the fluttering panel. 3.2. Convergence studies The aeroelastic equations (AEs) using the RRMs or the POMs can be transformed into a smaller number of ODEs. Specifically, how much the AEs can be reduced will be explored with the convergence studies. A modal participation (MP) value that evaluates the contribution of the k-th RRMs or POMs to the total deflection is defined as MP ¼ maxjqk j=∑Kj¼ 1 maxjqj j. Tables 2 and 3 show the LCO deflections and MP values of the first six RRMs and POMs for the panels of a/b¼0.5 and a/b¼1 under three different dynamic pressures. The two methods have a good agreement on the LCO deflections. The MP
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
1.5
12 RRMs 4 POMs
a/b=1, λ=250
10
1
12 RRMs 4 POMs
a/b=1, λ=250
5 dW/d τ
0.5
W(ξ,η ,τ)
6201
0 −0.5
0
−5
−1 10
12
14
τ
16
18
−10 −1.5
20
−1
0
0.5
1
1.5
W
12 RRMs 4 POMs
a/b=1,λ=250 1
Amplitude
−0.5
0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
2.5
3
3.5
4
f/Hz Fig. 12. Comparison of a LCO for a panel of a/b ¼ 1 at λ¼ 250 using 12 RRMs and 4 POMs: (a) time history; (b) phase portrait; (c) frequency spectrum.
values are decreasing rapidly for both the RRMs and the POMs. Specifically, the first two POMs are the most dominant and they contribute more than 90 percent to the deflection, with the last two contributing less than 0.1 percent and 1 percent for a/b¼0.5 and a/b ¼1, respectively. Overall, the POMs converge faster than the RRMs. Based on the convergence studies previously discussed, Fig. 10(a) and (b) plots the LCO deflection versus dynamic pressure using POMs and RRMs for comparison. Obviously, only 4 POMs are necessary to obtain a good agreement with that using 12 ðM ¼ 6; N ¼ 2Þ RRMs. Specially, in Fig. 10(b) there is a departure between the POD and RR results from λ ¼285 to λ ¼300, which is due to the duality bifurcation depending on the initial conditions rather than the inaccurate discrepancy between the two methods. Ref. [12] provides the detail analysis of the duality bifurcation. Figs. 11 and 12 show the time histories, phase portraits and frequency spectra for a/b¼0.5 at λ ¼50 and a/b ¼1 at λ ¼250, respectively. It turns out that the POD method can give accurate LCO solutions with fewer modes than the RR method. 3.3. Efficiency studies In addition to the requirement of fewer modes, a computational cost saving is also our concern using the POD method. For a quantitative evaluation of the cost savings, the efficiency of the POD method is explored in this section. Three panels of a/b¼0.5, 1 and 2 are considered. Both the RR and POD methods are employed to determine the LCO responses. Table 4 provides the comparison between the two methods in terms of the Mode number, Integration time (It), Time step (Ts), and CPU time (Ct) (Intel Corei7-2720QM at 2.20 GHz) in RK4. For the panels of a/b¼0.5 and 1, the CPU time using the RR method with 12 RRMs is approximately 30 times that using the POD method with 4 POMs. For the panel of a/b¼2, the POD method with 6 POMs is approximately 60 times faster than the RR method with 16 RRMs. Hence, the computational cost is vastly reduced with fewer modes and a larger time step using the POD method. By contrast, when the POD method with a mode projection (from the POMs to RRMs) is used, there will be a huge difference. For the sake of brevity, the results are not reported here. However, the numerical simulations indicated that the projection-free POD method is approximately two orders of magnitude faster than the POD method with the mode projection. 3.4. Chaos studies Generally, a linear stability analysis of a nonlinear system is necessary. However, such studies for this nonlinear system are voluminous in the literature and the linear stability results for the presently studied plate configurations are well known.
6202
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Table 4 Efficiency comparison between the Rayleigh–Ritz and the POD methods. Method
a=b
RR POD RR POD RR POD
0.5
λ 50
1
200
2
950
Mode number
It
Ts
Ct (s)
12 4 12 4 16 6
100
0.005 0.01 0.005 0.01 0.002 0.005
715.1 20.8 375.0 11.8 1473.5 24.8
60 30
3 2
5 a/b=0.5,μ/Ma=0.05 I=R=12,J=S=2 M=6,N=2
Local amplitude extrema
Local amplitude extrema
4
1 0 −1 −2 −3 −4 −5 20
40
60
80
100
120
140
a/b=0.5,μ/Ma=0.05 I=R=14,J=S=2 K=4
0
−5 20
160
40
60
80
λ
100
120
140
160
λ
Fig. 13. Comparison of bifurcation diagram for a panel of a/b ¼ 0.5 using 12 RRMs and 4 POMs: (a) bifurcation diagram (12 RRMs); (b) bifurcation diagram (4 POMs).
4 3
3
RR POD
λ=64
RR POD
λ=64
2 1
1
dW/dτ
dW/dτ
2
0 −1
0 −1
−2 −2
−3 −4 −2
−1.5
−1
−0.5
0
0.5
1
1.5
2
−3 −0.1
−0.05
0
0.05
0.1
W
W 1
RR POD
λ=64
0.8 0.6 0.4 0.2 0
0
0.2
0.4
0.6
0.8
1
Fig. 14. Comparison of a quasi-periodic oscillation for the panel of a/b ¼ 0.5 at λ¼ 64 using 12 RRMs and 4 POMs: (a) phase portrait; (b) Poincare plot; (c) frequency spectrum.
The linear stability results agree with the present nonlinear results in the limit as the amplitude of the limit cycle oscillations approaches zero. Therefore, instead of the linear stability study, the bifurcation diagrams are presented to observe the chaotic motions.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
10
6
RR POD
λ=90
λ=90
4
5
2
dW/dτ
dW/dτ
6203
0
0 −2
−5 −4 −10 −4
−3
−2
−1
0
1
2
3
RR POD
−6 −0.25 −0.2 −0.15 −0.1 −0.05
4
0
0.05
0.1
0.15
W
W
λ=90
RR POD
0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
Fig. 15. Comparison of a chaotic oscillation for the panel of a/b ¼ 0.5 at λ¼ 90 using 12 RRMs and 4 POMs: (a) phase portrait; (b) Poincare plot; (c) frequency spectrum.
15
RR POD
λ=110
10
10
5
5
dW/dτ
dW/dτ
15
0 −5
RR S−POD
λ=110
0 −5
−10 −10 −15
−4
−3
−2
−1
0
1
2
3
4
−0.25 −0.2 −0.15 −0.1 −0.05
W
0
0.05 0.1 0.15
W 0.5
RR POD
λ=110
Amplitude
0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
f/Hz Fig. 16. Comparison of a chaotic oscillation for the panel of a/b ¼ 0.5 at λ ¼110 using 12 RRMs and 4 POMs: (a) phase portrait; (b) Poincare plot; (c) frequency spectrum.
6204
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
3
15
λ=132
2 1
5
0
dW/dτ
W(ξ,η,τ)
λ=132
10
−1
0
−2
−5
−3
−10
−4
440
450
460
470
480
490
−15 −5
500
−4
−3
−2
−1
τ
0
1
2
3
W 1
λ=132
Amplitude
0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
f/Hz Fig. 17. A quasi-periodic oscillation for the panel of a/b ¼ 0.5 at λ¼ 132 using 4 POMs: (a) time response; (b) phase portrait; (c) frequency spectrum.
5
λ=132
15
4
10
2
5
dW/dτ
W(ξ,η,τ)
3 1 0 −1
0 −5
−2 −3 −4
λ=132
−10 420
430
440
450
460
470
480
490
−4
−3
−2
−1
0
τ
2
3
4
λ=132
1
Normalized amplitude
1
W
0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
f/Hz Fig. 18. An almost chaotic oscillation for the panel of a/b ¼ 0.5 at λ ¼ 132 using 12 RRMs: (a) time response; (b) phase portrait; (c) frequency spectrum.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
3 a/b=1,μ/Ma=0.05 I=R=14,J=S=2 M=8,N=2
2
Local amplitude extrema
Local amplitude extreme
3
6205
1 0 −1 −2 −3 300
350
400
450
a/b=1,μ/Ma=0.05 I=R=14,J=S=2 K=4
2 1 0 −1 −2 −3 300
500
350
400
450
500
λ
λ
Fig. 19. Comparison of bifurcation diagram for a panel of a/b ¼ 1 using 16 RRMs and 4 POMs: (a) bifurcation diagram (16 RRMs); (b) bifurcation diagram (4 POMs).
2
a/b=2, μ/Ma=0.05 I=R=16,J=S=2 M=8,N=2
1.5 1
Local amplitude extrema
Local amplitude extrema
2
0.5 0 −0.5 −1 −1.5 −2 950
1000
1050
1100
1150
1.5 1 0.5 0 −0.5 −1 −1.5 −2 950
1200
a/b=2,μ/Ma=0.05 I=R=16,J=S=2 K=6
1000
1050
λ
1100
1150
1200
λ
Fig. 20. Comparison of bifurcation diagram for a panel of a/b ¼ 2 using 16 RRMs and 6 POMs: (a) bifurcation diagram (16 RRMs); (b) bifurcation diagram (6 POMs).
1.5
12 RRMs 4 POMs
a/b=0.5
μ/Ma=0.05
1
μ/Ma=0.25
0.8
μ/Ma=0.5
12 RRMs 4 POMs
a/b=1 μ/Ma=0.05 μ/Ma=0.25
1
μ/Ma=0.5
Wp
Wp
1.2
0.6 0.5
0.4 0.2 0
20
25
30
35
40
λ
45
50
55
60
0
100
150
200
250
300
λ
Fig. 21. Comparison of LCO deflection peak versus dynamic pressure at varying μ=Ma ¼ 0:05, 0.25, 0.5 using 12 RRMs and 4 POMs: (a) a panel of a/b ¼ 0.5; (b) a panel of a/b ¼ 1.
Fig. 13 shows the bifurcation diagrams for the panel of a/b¼0.5 using 4 POMs and 12 RRMs. They are in a good agreement with a symmetry property in terms of the amplitude of response. The reason for the symmetry lies in the fact that the aeroelastic system of this cantilever plate is symmetric, and it leads to the symmetric solutions depending on the initial conditions. Specifically, Figs. 14, 15 and 16 show the phase portraits, Poincare plots and frequency spectra for dynamic pressures of λ ¼64, 90 and 110 respectively. Note that the phase portraits using the POD method shown in Figs. 14(a), 15(a) and 16(a) are plotted with the symmetric solutions of the original ones for clear comparison with the RR results. Along with the frequency spectra shown in Figs. 14(c), 15(c) and 16(c), the two methods are in a very good agreement. The Poincare maps using the two methods are approximately symmetric in Figs. 14(b), 15(b) and 16(b). All these indicate that the POD method can be successfully used for the analysis of chaotic oscillations.
6206
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
Exceptionally at λ ¼ 132, Fig. 17 shows a quasi-periodic motion with the POD method. By contrast, Fig. 18 indicates a nearly chaotic oscillation using the Rayleigh–Ritz method. For the panels of a/b¼1 and 2 shown respectively in Figs. 19 and 20, the bifurcation diagrams using the two methods have great similarities but with differences. The possible reasons are twofold. First, the chaotic oscillation is extremely sensitive to the dynamic system, so the approximate POD/ROM may differ from the original system, which would lead to distinct chaotic motions. Another possible reason is that the POMs extracted from the chaotic results are not unique. Maybe the POMs we have determined are not the optimal ones, which could affect the resulting chaotic motions. How to obtain the optimal POMs will be a subject for the future work.
3.5. Robustness of the POD modes for various μ=Ma Both λ and μ=Ma are the typical and important flight parameters. If the POMs can be applied to the variations of the flight parameters in a range of interest, the systems are easy to be analyzed without remodeling by the new POMs determined. The numerical simulations performed previously have shown that the POMs can be applied for different λ. Hence, a question arises here: whether a set of general POMs can be used for various μ=Ma? This is explored next. The general POMs are determined at a fixed μ=Ma ¼ 0:05 with the panel undergoing chaotic oscillations. Fig. 21(a) and (b) shows the LCO deflection versus dynamic pressure using the general POMs at μ=Ma ¼ 0:05, 0.25 and 0.5 for the panels of a/b¼0.5 and 1, respectively. They are all in an excellent agreement with the Rayleigh–Ritz results. Therefore, the question arisen previously can be answered as the general POMs determined at a fixed μ=Ma can be applied for different μ=Ma. Consequently, a conclusion may be drawn that as long as the POMs are determined from the chaotic responses, the POMs are robust for the same panel at varying flight parameters. The same panel is defined as a panel with the same dimensions and support conditions. Said another way, the POMs are the inherent characteristics of the fluttering cantilever panel independent of the flight conditions in this paper.
4. Concluding remarks A computationally fast POD method has been proposed and employed to study the aeroelastic response of a cantilever wing-like plate in supersonic flow. The presently proposed POD method has been verified to be about two orders of magnitude faster than the conventional POD method, due to the fact that the present method avoids the expensive projection procedure from the POD modes to the Rayleigh–Ritz modes by performing direct numerical derivatives to the POD modes. In addition to the computational efficiency improvement, the present projection-free POD method has also been compared with the Rayleigh–Ritz method for investigating the fluttering plate with three length–width ratios, namely a/b¼0.5, 1 and 2. Drawn from numerical simulations are the following conclusions: 1. For a good convergence, the required number of POD modes is much less than that of the Rayleigh–Ritz modes. Concretely, 4, 4 and 6 POD modes in contrast to 12, 12 and 16 Rayleigh–Ritz modes are required in analyzing the panels of a/b¼0.5, 1 and 2, respectively. Hence, the computational effort has been remarkably reduced using the POD method (1/ 60 1/30 that by the Rayleigh Ritz method). Also, numerical simulation indicates that although more POD modes should be included as a=b is increased for both methods, the POD method converges faster. 2. It was shown that the dominant (or first) POD mode extracted from the Rayleigh–Ritz chaotic results resembles the physical deflection of the panel, which accounts for the fact that fewer POD modes are required for capturing the panel deflection. 3. It has been demonstrated that the POD modes determined from the chaotic results for a specified set of parameters are universal modes, and are applicable to a variety of parametric variations (e.g. varying λ, and μ=Ma). This suggests that the POD modes are inherent characteristics of a plate with a given geometry and supporting pattern. 4. Previous studies believed that the POD method only works for capturing the periodic responses of nonlinear dynamical systems. In this study, the POD method has been shown to be able to obtain the periodic solutions as well as the chaotic solutions.
Acknowledgments The first author gratefully acknowledges the helpful discussions with the aeroelasticity research group in the Department of Mechanical Engineering and Materials Science, School of Engineering at Duke University, and colleagues of the College of Astronautics in Northwestern Polytechnical University.
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
6207
Appendix A. Coefficients of C, D, A, B, Q and F
C ij ¼
Drs ¼
3 K K Z Z 2 K K Z Z h ∂φk1 ∂φk2 dui h h ∂φk1 ∂φk2 dui ∑∑ vj dx dy ν v dx dy ∑∑ a k1 k2 b a k1 k2 ∂x ∂x dx ∂y ∂y dx j 2 K K Z Z h h ∂φk1 ∂φk2 dvj ∑∑ u dx dy; ð1 νÞ b a k1 k2 ∂x ∂y i dy 3 K K Z Z 2 K K Z Z h ∂φk1 ∂φk2 dvs h h ∂φk1 ∂φk2 dvs ∑∑ ur dx dy ν ur dx dy ∑∑ b a b ∂y ∂y dy ∂x ∂x dy k1 k2 k1 k2 2 K K Z Z h h ∂φk1 ∂φk2 dur ∑∑ vs dx dy; ð1 νÞ a b k1 k2 ∂x ∂y dx Akk1 ¼
Bkk1 ¼
1 6
Z Z
1 6
Z Z
φk1 φk dx dy;
"
λ K ∑q 6β k1 k1
Z Z
∂φk1 Ma2 2 φk dx dy þ 2 ∂x Ma 1
rffiffiffiffi
μ K dqk1 ∑ λ k1 dt
(23)
(24)
a 4 Z Z ∂ 2 φ ∂ 2 φ
a 2 Z Z ∂2 φ ∂2 φ ∂2 φk1 ∂2 φk ∂2 φk1 ∂2 φk k1 k k1 k dx dy þ dx dyþ ν dx dy þ dx dy b b ∂x2 ∂x2 ∂y2 ∂y2 ∂y2 ∂x2 ∂x2 ∂y2
a 2 Z Z ∂2 φ ∂2 φ k1 k dx dy ; þ 2ð1 νÞ b ∂x∂y ∂x∂y Qk ¼
(22)
Z Z
(25)
#
φk1 φk dx dy ;
(26)
Z Z K K K dui ∂φk1 ∂φk ∂φk1 ∂φk2 ∂φk3 ∂φk vj dx dy þ ∑ ∑ ∑ qk1 qk2 qk3 dx dy h i j k1 dx ∂x ∂x ∂x ∂x ∂x ∂x k1 k2 k3 Z Z
a 3 a R S K dvs ∂φk1 ∂φk þ2 ∑ ∑ ∑ brs qk1 dx dy ur b h r s k1 dy ∂y ∂y Z Z Z Z
a 4 K K K
a 2 a I J K ∂φk1 ∂φk2 ∂φk3 ∂φk dui ∂φk1 ∂φk þ ∑ ∑ ∑ aij qk1 dx dy þ 2ν v dx dy ∑ ∑ ∑ qk1 qk2 qk3 b k1 k2 k3 b h i j k1 ∂y ∂y ∂y ∂y dx j ∂y ∂y Z Z
a 2 K K K ∂φk1 ∂φk2 ∂φk3 ∂φk dx dy ∑ ∑ ∑ qk1 qk2 qk3 þν b k1 k2 k3 ∂y ∂x ∂x ∂y Z Z Z Z
a a R S K
a 2 K K K dvs ∂φk1 ∂φk ∂φk1 ∂φk2 ∂φk3 ∂φk þ 2ν ∑ ∑ ∑ brs qk1 dx dyþ ν dx dy ∑ ∑ ∑ qk1 qk2 qk3 ur b h r s k1 b k1 k2 k3 dy ∂x ∂x ∂x ∂y ∂y ∂x ( Z Z
a 2 a I J K dvj ∂φk1 ∂φk ∂φk1 ∂φk þ ð1 νÞ þ dx dy ui ∑∑∑a q b h i j k1 ij k1 dy ∂y ∂x ∂x ∂y Z Z
a a R S K dur ∂φk1 ∂φk ∂φk1 ∂φk þ ∑ ∑ ∑ brs qk1 vs þ dx dy b h r s k1 dx ∂y ∂x ∂x ∂y ) Z Z
a 2 K K K ∂φk1 ∂φk2 ∂φk3 ∂φk ∂φk1 ∂φk2 ∂φk3 ∂φk þ dx dy : (27) ∑ ∑ ∑ qk1 qk2 qk3 þ b k1 k2 k3 ∂x ∂y ∂y ∂x ∂x ∂y ∂x ∂y
Fk ¼ 2
a
I
J
K
Z Z
∑ ∑ ∑ aij qk1
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
E.H. Dowell, Nonlinear oscillations of a fluttering plate, AIAA Journal 4 (1966) 1267–1275. E.H. Dowell, Nonlinear oscillations of a fluttering plate. II, AIAA Journal 5 (1967) 1856–1862. E.H. Dowell, Aeroelasticity of Plates and Shells, Noordhoff International Publishing, Groningen, Netherlands, 1975. C. Mei, A finite-element approach for nonlinear panel flutter, AIAA Journal 15 (1977) 1107–1110. W.L. Ye, E.H. Dowell, Limit cycle oscillation of a fluttering cantilever plate, AIAA Journal 29 (1991) 1929–1936. J. Zhou, Z.C. Yang, Y.S. Gu, Aeroelastic stability analysis of heated panel with aerodynamic loading on both surfaces, Science China: Technological Sciences 55 (2012) 2720–2726. D. Xie, M. Xu, E.H. Dowell, Proper orthogonal decomposition reduced-order model for nonlinear aeroelastic oscillations, AIAA Journal 52 (2014) 1–13. H.H. Dai, M. Schnoor, S.N. Atluri, A simple collocation scheme for obtaining the periodic solutions of the duffing equation, and its equivalence to the high dimensional harmonic balance method: subharmonic oscillations, Computer Modeling in Engineering and Sciences 84 (2012) 459. H.H. Dai, X.K. Yue, J.P. Yuan, S.N. Atluri, A time domain collocation method for studying the aeroelasticity of a two dimensional airfoil with a structural nonlinearity, Journal of Computational Physics 270 (2014) 214–237. H.H. Dai, J.K. Paik, S.N. Atluri, The global nonlinear Galerkin method for the analysis of elastic large deflections of plates under combined loads: a scalar homotopy method for the direct solution of nonlinear algebraic equations, Computers Materials and Continua 23 (2011) 69–99.
6208
D. Xie et al. / Journal of Sound and Vibration 333 (2014) 6190–6208
[11] H.H. Dai, J.K. Paik, S.N. Atluri, The global nonlinear Galerkin method for the solution of von Karman nonlinear plate equations: an optimal & faster iterative method for the direct solution of nonlinear algebraic equations FðxÞ ¼ 0, using x ¼ λ½αF þ ð1 αÞBT B, Computers Materials and Continua 23 (2011) 155–185. [12] D. Xie, M. Xu, H. Dai, E.H. Dowell, Observation and evolution of chaos for a cantilever plate in supersonic flow, Journal of Fluids and Structures (2014), http://dx.doi.org/10.1016/j.jfluidstructs.2014.05.015, in press. [13] I.R. Dixon, C. Mei, Finite element analysis of large-amplitude panel flutter of thin laminates, AIAA Journal 31 (1993) 701–707. [14] C.E. Gray, C. Mei, Large-amplitude finite element flutter analysis of composite panels in hypersonic flow, AIAA Journal 31 (1993) 1090–1099. [15] P.J. Attar, E.H. Dowell, A reduced order system id approach to the modelling of nonlinear structural behavior in aeroelasticity, Journal of Fluids and Structures 21 (2005) 531–542. [16] X. Guo, C. Mei, Using aeroelastic modes for nonlinear panel flutter at arbitrary supersonic yawed angle, AIAA Journal 41 (2003) 272–279. [17] X. Guo, C. Mei, Application of aeroelastic modes on nonlinear supersonic panel flutter at elevated temperatures, Computers & Structures 84 (2006) 1619–1628. [18] B.I. Epureanu, L.S. Tang, M.P. Paidoussis, Coherent structures and their influence on the dynamics of aeroelastic panels, International Journal of NonLinear Mechanics 39 (2004) 977–991. [19] B.I. Epureanu, L.S. Tang, M.P. Paidoussis, Exploiting chaotic dynamics for detecting parametric variations in aeroelastic systems, AIAA Journal 42 (2004) 728–735. [20] M. Amabili, A. Sarkar, M.P. Paidoussis, Reduced-order models for nonlinear vibrations of cylindrical shells via the proper orthogonal decomposition method, Journal of Fluids and Structures 18 (2003) 227–250. [21] M. Amabili, A. Sarkar, M.P. Paidoussis, Chaotic vibrations of circular cylindrical shells: Galerkin versus reduced-order models via the proper orthogonal decomposition method, Journal of Sound and Vibration 290 (2006) 736–762. [22] D. Xie, M. Xu, A simple proper orthogonal decomposition method for von Karman plate undergoing supersonic flow, Computer Modeling in Engineering and Sciences 93 (2013) 377–409. [23] E.H. Dowell, On asymptotic approximations to beam modal functions, Journal of Applied Mechanics 51 (1984) 439. [24] A. Sarkar, M.P. Païdoussis, A cantilever conveying fluid: coherent modes versus beam modes, International Journal of Non-Linear Mechanics 39 (2004) 467–481.