- Email: [email protected]
Finite element analysis of composite panel flutter
Compufers & Sfrurrurrs Vol. 39, No. l/2, pp. 165-172. 1991 Printed in Great Britain.
FINITE
0
004s7949/91 53.00 + 0.00 1991 Pergmnon Press plc
ELEMENT ANALYSIS OF COMPOSITE PANEL FLUTTER I. Laat and M. H. CHO
Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, P.O. Box 150, Cheonryang, Seoul, Korea (Received
5 March
1990)
Abstract-The finite element method based on the shear deformable theory is developed to analyze the composite panel flutter in supersonic flow. The computational results of the vibration and flutter analysis agree well with those given in the available references. Guyan reduction and the normal mode method are used to reduce the computational time. The plate length ratio, the flow direction and the fiber orientation greatly affect the.flutter boundaries of trapezoidal laminated plates. Flutter boundaries have been calculated for both simply supported and clamped boundaries in order to determine the effect of boundary conditions.
NOTATION
aerodynamic damping and force matrices in finite element thicknessshear stiffnesses and nondimensionalized ones (i,j = 4,5) plate lengths in x- and y-directions bending stiffnesses and nondimensionalized ones (i, j = 1,2,6) nodal forces in finite element nondimensional aerodynamic damping parameter LLtLatrix shear correction coefficients (i = 4,5) nondimensional eigenvalues and critical flutter ones stiffness and condensed stiffness matrices in finite element mass matrix in finite element aerodynamic pressure transformed reduced stiffnesses generalized displacements: displacement component in z-direction, rotations in xz- and yz-planes, respectively airflow velocity Cartesian coordinates length-to-thickness ratio (a/h) and aspect ratio (a/b) nondimensional dynamic pressure parameter and critical one fiber angle nondimensional natural frequency parameter flow angle between flow direction and x-axis densities of plate and air nondimensional Cartesian coordinates direction of a supersonic airstream
1. INTRODUCTION
flutter is the dynamic aeroelastic instability of the external skin of a flight vehicle exposed to an airflow on one side. Fung [l] and Dowel1 [2] have studied the panel flutter for isotropic materials, both Panel
tTo whom correspondence should be addressed.
theoretically and experimentally. The numerical tools of theoretical solutions were Galerkin-type methods. Olson [3,4] and Sander et al. [5] investigated the flutter problem for simply supported and clamped isotropic panels by the finite element method. For the orthotropic plate, Ketter [6] has treated the effects of panel boundary conditions and angle of orthotropicity about the flow direction on the flutter boundaries. In recent years, composite materials have been used for the high performance and weight minimization of advanced aircraft. Srinivasan and Babu [7j investigated composite panel flutter problems for the cross-ply laminates by the use of integral equations. However, the effects of the fiber orientation and the flow direction on the flutter boundaries have not been investigated. In the present paper, the flutter phenomena for a composite plate in supersonic flow have been analyzed by the finite element method based on the shear deformable theory. For trapezoidal composite plates with simply supported and clamped edges, the effects of the geometry shape, the flow direction and the fiber orientation on the flutter boundaries have been obtained. 2.
THEORETICAL
FORMULATION
Figure 1 shows a composite panel and its coordinate systems. The panel is idealized as a thin plate which has one side exposed to a supersonic airflow and the other side to still air. The shear deformable theory is based on the heterogeneous laminated plate theory of Yang, Norris and Stavsky (YNS) [S]. The YNS theory is a generalization of Mindlin’s theory [9] for homogeneous, isotropic plates to arbitrarily laminated anisotropic plates and includes shear deformation and rotary inertia effects. In the present paper, symmetric laminated plates are considered and there will be no coupling between 165
166
I. LEEand M. H.
CHO
Fig. 1. Coordinate systems of composite panel. bending and extension of a laminate. Therefore, the nondimensional governing equations are expressed below in terms of the nondimensional displacements G, S, and $,
where
A; = A,,/E,h
(i,,j = 4, 5)
(i,.; = I. 2. 6)
0: = D,,lW’
s ill
4, = 0, (1)
(Q,,), dz ii?
i, j = 4, 5
I!2 4 =
s Ii 2
(p,,), z’ d-_ i.j=
1,2,6
I!2 (1. z’)p dz = (ph, ph3/12). @, 0 = s A2 Here, q(.u, J-,t) is the transversely distributed load; stiKnesses, D,, are bending stiffnesses; (&,,), are the transformed reduced stiffness for the k th layer of the fiber angle tIk; ki are the shear correction coefficients; p is the density of laminated plate and h is the thickness of laminated plate; and E, is Young’s modulus of any reference layer in the r-direction along the fiber or perpendicular to the fiber. The finite element model of the panel flutter is obtained using the shear deformable plate elements of Reddy [IO]. The midplane fI of the plate is subdivided into a finite number of elements, R, (e = I, 2 .). Over each element R,, the nondimensional generalized displacements (6, s,, s,.) are interpolated by expressions of the form A,, are thickness-shear
(4)
167
Finite element analysis of composite panel flutter where tii are the interpolation functions and n the number of nodes per element. Therefore, the finite element equations can be written as
Thus, the number of equations to be solved is reduced and the final finite element equation is expressed only in terms of {C}
(6) [;z;,
I::;
Z]{
ij}+u4$
[i2.~~Z,IIll{~~}=~i;n9.
(5)
where {F} are nodal forces due to the aerodynamic forces and [K] and [M] are given by
where [a is the condensed stiffness matrix and (M] = [Ml’]. The aerodynamic pressure q at a point (x, y) in supersonic flow is expressed as [l 1]
dx,
[P2] = x2{A:,[S’q + A&y[S@]} [P3] = a2{A,*,[S@j+ A&y[S”O]} [K22]= D:,[see] + D:6y([s@]
t) =
X [g
[K”] = L+{Af,[S”] + A&y([SCg] + [SC”lr) + ‘4$.4y2[P”]}
Y,
2% -(M2
_
1),,2
+;($$)$I
(M > 21’2) (7)
where q. is the dynamic pressure; M is the Mach number; V is the flow velocity; and [ the direction of the supersonic airstream on the surface of the panel (see Fig. 1). Let p = 2q,/(W - I)“2
+ [Ssq]rl + D&y2[S”“] + a2A,*,[S] [K23]= Dfzy[Sa] + Df6[Sss] + D:,y’PJ
+ D&y[SS”Jr+ dA,,[S]
[K3’] = D&[SSr] + D&y([Ssg] + [Seq]l) + D:2y2[S9”] + cc’A&[S]
where C$ is the angle between the [-direction airllow and the x-axis. The nodal forces (F} of eqn (6) are given by
&, {F) = -B*[A,l{@} ,
[M”] = [S]
[fP]
= [P]
pq$ =
Si,=
= &
I ap
z%dtdn
n,
SF=
a4
[S]
(9)
where
(p,q=O,t,tj)
s n,
--LZkl{~tj,
of
_ w4 p=E,hj [A[] = cos t#~ [Sot] + sin by [So7
tit+,dt dq,
where R, is the region of a typical element. Equations (5) can be simplified by neglecting the inertia terms of two rotations (S,, 3,) because [M*‘] and [M3’] become very small compared with [Ml’] for the thin panel (0: > 100). {??,} and {s,} can be written in terms of {@} through the Guyan reduction procedure.
HA = PI = [Ml, where /3* is the nondimensional dynamic pressure parameter. Substituting {F} in eqn (6) and assuming vc = d* e”‘, we can obtain the following equation:
[IQ + B*[A,l -
k**WIl{+*}
= {O>,
(10)
168
I. LEE and M. H. CHO
where
where [E] is an n x m normalized modal matrix obtained from the free vibration analysis and {b} contains m normal coordinates. Here, n is the dimension of the matrix in eqn (lo), and m is the number of normal modes chosen for modal analysis. Therefore, the final finite element equation can be simplified to the following equation:
k*Z= _;*Z . -P’* ;(*2
=w2a4P
h E,h’
/l* =
po.
Here, k * are the nondimensional eigenvalues, and i * are the nondimensional natural frequency parameters. The problem is reduced to that of finding the nondimensional eigenvalues and corresponding modal shapes of the system for a given value of the nondimensional dynamic pressure parameter. When p * = 0, the eigenvalues k* are real and positive, since [KJ and [M] are symmetric and positive definite. HOWever, the aerodynamic matrix [A,] is not symmetric, and hence complex eigenvalues k* are expected for p*>O. As /3* increases monotonically from zero, two of these eigenvalues will approach each other and coalesce to k,t at /I* = Bz and become complex conjugate pairs k*=kj!+ikT
(11)
for /?* > /I:,. Here, p,* corresponds
Wl* + 8*M:l- k*2Pfl*l{b) = (01, where [8* =
[El’[B PI = [A**1
WI* = Pl%W [El = [II &I* = [El%$l[El. 3. NUMERICAL RESULTS AND DISCUSSION
The following elastic properties of T300/5208 graphite/epoxy composite material are used:
to the value of
k* at which the first coalescence occurs. In general, the complex eigenvalues k* may be put in more
E, =20x
convenient
Ez=E,=1.4x
form k*’ = -g,(w/o,,)
- (w/oJ,
(12)
IO’psi
106psi
G,z = G,, = 0.8 x IO6 psi
where
Gz3 = 0.6 x lo6 psi
(M’-2)
p,V
v,* = 0.3
go = (IV* - 1)’ * phw,, ply thickness
h, = 4.92 x lo-’ in.
density p = 0.148 x IO- ’ lb-sec2/in4 Here, g, is the nondimensional aerodynamic damping parameter, and o,, a convenient frequency scale. When flutter occurs, w becomes a pure imaginary number. Hence, the flutter boundary (p,) may be obtained as a function of the aerodynamic damping g, by continuing the calculations for fl* > /I::. In the absence of aerodynamic damping, the flutter boundaries simply correspond to ,!?,*,and i * = ik *. Equation (10) can be solved directly by increasing values of fi* in many steps starting from zero. However, this procedure requires a tremendous amount of computational time. Therefore, the modal analysis is very useful in this case. It is generally seen that only a small number of eigenvectors in the lower modes gives very accurate results, since the first coalescence probably occurs among the first few natural frequencies. For this reason, {@*} can be approximated by the linear combinations of a limited number of normal modes
of the system. {I?*} Lz [E]{b},
(13)
(14)
Fig. 2. Geometry of trapezoidal laminated plate.
169
Finite element analysis of composite panel flutter
Table 1. Convergence of natural frequencies and flutter boundaries for [(O/90),], boron/epoxy rectangular clamped laminate (4 = 0, a/b = 1.0, c/b = 1.0) No. of modes
Mesh 2x2 3x3 4x4 5x5 6x6 7x7 8x8
11
Srinivasan and Babu [7] PIE method Series solution Thornton and Clary [ 121 $RR method FE method Exwrimental
1 24.97 23.61 23.42 23.37 23.35 23.34 23.34
Natural 2 56.71 45.77 43.24 42.65 42.44 42.35 42.30
frequency parameter I * 3 4 5 74.31 90.76 58.67 77.87 95.50 54.98 67.13 84.42 54.12 66.03 78.92 53.81 65.65 77.28 53.68 65.48 76.58 53.62 65.40 76.25
23.33 23.63
42.36 42.28
53.77 53.76
65.56 65.42
76.52 75.89
92.76 92.19
22.71 23.25 17.44
44.14 45.23 35.61
49.59 50.68 43.42
63.58 64.31 54.68
80.84 83.20 64.85
94.10 95.19 79.75
Flutter bounds k2: B:: 55.46 775.62 47.55 494.00 46.43 467.22 46.55 466.29 46.72 468.71 46.83 470.29 46.89 471.16
6 111.60 100.27 95.02 93.44 92.78 92.46
46.09 47.19
446.36 474.60
-
-
-
t IE method: integral equation method. Series solution: solution of double series of beam characteristic functions. $ RR method: Rayleigh-Ritz method. FE method: finite element method. Table 2. Comparison of natural frequencies and flutter boundaries with [7] for [(90/O),], boron/epoxy trapezoidal plate (I#I= 0, a/b = 0.75, c/b = 0.2) Reference
Natural frequency parameter 1* 3 4 5 105.20 114.84 166.63 180.71 105.51 118.85 168.48 187.27
1
Present method Srinivasan and Babu 171
2
58.14 58.14
The stiffness coefficients associated with A, (i,j = 4, 5) are evaluated using a reduced integration technique to circumvent the computational difficulty (i.e. to avoid the shear locking). A nine-node isoparametric, shear deformable plate element is used in this analysis. The geometry of the trapezoidal laminated plate considered in this paper is shown in Fig. 2. Convergence study
Table 1 shows the convergence characteristics of the solution for a typical [(O/90),], boron/epoxy rectangular clamped laminate (c/b = 1.0, a/b = 1.0). The following material properties of boron/epoxy composite are used for the convergence test: E, = clomped
;
191.28 199.73
94.78 96.76
999.07 1007.80
31.0 x lo6 psi, E2 = 2.7 x lo6 psi, Cl1 = 0.75 x lo6 psi, vIZ= 0.28 and p = 0.192 x 10m3lb-sec2/in4. In Table 1, the mesh size was varied from 2 x 2 to 8 x 8, keeping the number of modes constant at 11. It may be observed that the convergence of flutter boundaries is not monotonic, which is in contrast to the self-adjoint problem (free vibration) where the convergence of eigenvalues is consistently monotonic from one side. However, there is a good convergence of solutions for 7 x 7 or 8 x 8 mesh. The present results agree well with those in the available references, as can be seen from Table 1. The effect of number of modes on the flutter boundaries for fixed mesh size is examined in this analysis. There is a good convergence for 10 or
, C$= 0
c/b=0.4
Flutter bounds k,: B,:
6
simply
: c/b=0.4
supported
,@ =70
100.
:I:
F---_=i
80 -
*:I,““;,o”, 0
200
, , q 400
P’
600
800
,“;
mode
4
20 -
1 -
[45/*45/01, 0
I/.
0
I
I 200
I
1
,
400
,
600
]
,
,
10
ok; - -10
600
P’
Fig. 3. Coalescence of eigenvalues for composite trapezoidal plate (graphite/epoxy; a/b = 1.O).
-20
1. LEE and M. H. CHO
170
flutter
mode
(real
mode
2
flutter
shape
part)
mode
(imaginary
shape part)
I and 2 and flutter mode shape for [O/*45/0], graphite/epoxy trapezoidal laminate (clamped: c,+= 0, a/b = 1.0, c/b = 0.4).
Fig. 4. The evolution
of modes
11 modes. Hence, for further analysis, a mesh size of 7 x 7 and 11 modes have been chosen. Table 2 shows
that the present results agree well with those from [7] for a trapezoidal cross-ply panel. Numerical application
For graphite/epoxy trapezoidal laminated plates, flutter boundaries are calculated when the length-tothickness ratio (a/h) is 200, with the panel length (a) in the x-direction held fixed. Figure 3 shows the coalescence of eigenvalues for [O/+45/O], and [45/ +45/O], trapezoidal laminates. For [O/+45/0], and [45! *45/O], panels, the first coalescences occur at the merging points of modes 1 and 2 and of modes 3 and 4, respectively. In Fig. 4. the evolution pattern of modes 1 and 2 and the flutter mode shape with the increase of dynamic pressure up to /I* = 900 for a [O/&-45/O],panel are depicted. Two different modes (modes 1 and 2) become a single flutter mode at B: and complex conjugate pairs for p * > p:.
This shows that the lateral deflection is concentrated near the rear of the panel. Figures 5 and 6 illustrate the effect of the fiber angle for [B/+45/0], trapezoidal laminates with simply supported and clamped edges (4 = 0). The flutter boundary (fix) has the maximum value at 0 = 0 and decreases gradually as the fiber angle (0) increases. Flutter boundaries decrease as the lower to the upper side-length ratio (c/h) increases. Flutter boundaries for laminates with clamped boundaries are higher than those with simply supported boundaries. as we expect. Figures 7 and 8 show the effect of the flow direction for [Q/_t 45/O],) (when 0 = 45 and 90, respectively) laminates with simply supported edges. The maximum value of flutter boundaries appears in the range of C$= 40-70 at the [45/+45/O], laminate, and in the range of 4 = 80-90 at the [90/+45/O], laminate. This illustrates that the maximum flutter boundaries appear in the vicinity of C$= 0 (i.e. when the flow direction is consistent with the fiber angle).
50
k’,,
4o 30
01
0
b:c/b=O4
o.c/b=O.B
q:c,‘b=0.6
‘I : c/b=
15
1
/
I
30
45
60
Fiber
Fig.
5. Effect
of fiber
Angle
I 1 0
I 75
90
0
orientation on flutter boundaries trapezoidal panels (simply
‘“i0,--__! 15
30
Fiber
for [0/*45/O], supported).
45
60
Angle
graphite/epoxy
75
90
8
laminated
171
Finite element analysis of composite panel flutter
60
600
k*cr 40
Fiber
Fig. 6. Effect of fiber orientation
:oc
k*
CT
on flutter boundaries for [0/*45/O], trapezoidal panels (clamped).
[
1
1200
Angle
@
graphite/epoxy
laminated
[
-I
50 40
Flow
Angle
9,
Flow
Angle
@
Fig. 7. Effect of flow direction on flutter boundaries for [45/k 45/O], laminated trapezoidal panels (simply supported; a/b = 1.0).
600 40 400
0
0
10
20
30
Flow
40
50
Angle
60
70
60
90
@
0
0
10
20
30
40
Flow
50
60
Angle
70
60
90
@
Fig. 8. Effect of flow direction on flutter boundaries for [90/+45/O], laminated trapezoidal panels (simply supported; a/b = 1 .O). Finally, Figs 9 and 10 show the flutter boundary (j?:) as a function of the aerodynamic damping g, for [60/ _+45/O], and [90/&45/O], laminates, respectively. Generally, as the aerodynamic damping increases, the flutter boundaries increase. If the actual flight conditions are known, these graphs allow one to define the flutter boundaries including the effect of the aerodynamic damping. From these figures, we see that the boundary condition affects the flutter boundary
significantly. The flutter dynamic pressure for a clamped boundary is much higher than that for a simply supported boundary. 4. CONCLUSIONS
The application of the finite element method based on the shear deformable theory to supersonic flutter problems has been presented. The results of the
172
I. LEE and
M. H.
CHO [ 60/*45/O]s
[60/*45/O]s
1200
02 /
0
20
Fig. 9. Effect of aerodynamic
40
60
damping
0
on flutter dynamic
pressure
20
00
60
c/h (Q, = 0, o/h = 1.O)
with various
1600
1200
I”, p;
800
0
fl;
800
400 simply
supported
[90/145/O], i-
boundary
(mped
bou,ndc;y
/
,
1 : Fig.
10. Effect
of aerodynamic
damping on flutter dynamic (c/h = 0.4, a/b = 1.0).
present method agree well with those in the available references. It is found that the influence of the panel length ratio, the flow direction, the fiber orientation and boundary conditions greatly affect the flutter boundaries of trapezoidal laminated plates. For [e/f 45/O], trapezoidal laminates, the maximum values of the flutter dynamic pressure appear near the point at which the fiber angle (0) is consistent with the flow angle (4), and the flutter boundaries increase as the ratio of the lower side length to the upper side length (c/b) decreases. The flutter dynamic pressure for laminates with a clamped boundary is much higher than that for those with a simply supported boundary.
REFERENCES
1. Y. C. Fung, Some recent contributions to panel flutter research. AIAA Jnl 1, 898-909 (1963). 2. E. H. Dowell, Panel flutter: a review of the aeroelastic stability of plates and shells. AIAA Jnl 8, 386-399 (1970).
pressure
20 with
40
60
various
How
angles
3. M. D. Olson, Finite element applied to panel flutter. AIAA Jnl 5, 2261-2270 (1967). 4. M. D. Olson, Some flutter solution using finite elements. AIAA Jnl8, 747-152 (1970). Finite element 5. G. Sander, C. Bon and M. Geradin, analysis of supersonic panel flutter. In!. J. Numer. Mefh. Engng I, 379-394 (1973). orthotropic 6. D. J. Ketter, Flutter of flat, rectangular, panels. AIAA Jnl 8, 116-124 (1967). 7. R. S. Srinivasan and B. J. Babu, Free vibration and flutter of laminated quadrilateral plates. Compur. Struct. 27, 2977304 (1987). 8. C. Yang, C. H. Norris and Y. Stavsky, Elastic wave propagation in heterogeneous plates. Int. J. Solids Strucf. 2, 665-684 (1966). 9. R. D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates. J. uppl. Mech. 18,31. 38 (1951). element for the 10. J. N. Reddy, A penalty plate-bending analysis of laminated anisotropic composite plates. Ini. J. Numer. Meth. Engng 15, 118771206 (1980). Principles of Aero11. H. Ashley and R. L. Bisplinghoff, elasticity. Dover, New York (1975). and R. R. Clary, Correlation study of 12. E. A. Thornton finite element modeling for vibrations of composite material panels. In Composite Materials: Testing and Design, ASTM STP 546, pp. 111-129. American Society for Testing Materials (1974).
FINITE
0
004s7949/91 53.00 + 0.00 1991 Pergmnon Press plc
ELEMENT ANALYSIS OF COMPOSITE PANEL FLUTTER I. Laat and M. H. CHO
Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, P.O. Box 150, Cheonryang, Seoul, Korea (Received
5 March
1990)
Abstract-The finite element method based on the shear deformable theory is developed to analyze the composite panel flutter in supersonic flow. The computational results of the vibration and flutter analysis agree well with those given in the available references. Guyan reduction and the normal mode method are used to reduce the computational time. The plate length ratio, the flow direction and the fiber orientation greatly affect the.flutter boundaries of trapezoidal laminated plates. Flutter boundaries have been calculated for both simply supported and clamped boundaries in order to determine the effect of boundary conditions.
NOTATION
aerodynamic damping and force matrices in finite element thicknessshear stiffnesses and nondimensionalized ones (i,j = 4,5) plate lengths in x- and y-directions bending stiffnesses and nondimensionalized ones (i, j = 1,2,6) nodal forces in finite element nondimensional aerodynamic damping parameter LLtLatrix shear correction coefficients (i = 4,5) nondimensional eigenvalues and critical flutter ones stiffness and condensed stiffness matrices in finite element mass matrix in finite element aerodynamic pressure transformed reduced stiffnesses generalized displacements: displacement component in z-direction, rotations in xz- and yz-planes, respectively airflow velocity Cartesian coordinates length-to-thickness ratio (a/h) and aspect ratio (a/b) nondimensional dynamic pressure parameter and critical one fiber angle nondimensional natural frequency parameter flow angle between flow direction and x-axis densities of plate and air nondimensional Cartesian coordinates direction of a supersonic airstream
1. INTRODUCTION
flutter is the dynamic aeroelastic instability of the external skin of a flight vehicle exposed to an airflow on one side. Fung [l] and Dowel1 [2] have studied the panel flutter for isotropic materials, both Panel
tTo whom correspondence should be addressed.
theoretically and experimentally. The numerical tools of theoretical solutions were Galerkin-type methods. Olson [3,4] and Sander et al. [5] investigated the flutter problem for simply supported and clamped isotropic panels by the finite element method. For the orthotropic plate, Ketter [6] has treated the effects of panel boundary conditions and angle of orthotropicity about the flow direction on the flutter boundaries. In recent years, composite materials have been used for the high performance and weight minimization of advanced aircraft. Srinivasan and Babu [7j investigated composite panel flutter problems for the cross-ply laminates by the use of integral equations. However, the effects of the fiber orientation and the flow direction on the flutter boundaries have not been investigated. In the present paper, the flutter phenomena for a composite plate in supersonic flow have been analyzed by the finite element method based on the shear deformable theory. For trapezoidal composite plates with simply supported and clamped edges, the effects of the geometry shape, the flow direction and the fiber orientation on the flutter boundaries have been obtained. 2.
THEORETICAL
FORMULATION
Figure 1 shows a composite panel and its coordinate systems. The panel is idealized as a thin plate which has one side exposed to a supersonic airflow and the other side to still air. The shear deformable theory is based on the heterogeneous laminated plate theory of Yang, Norris and Stavsky (YNS) [S]. The YNS theory is a generalization of Mindlin’s theory [9] for homogeneous, isotropic plates to arbitrarily laminated anisotropic plates and includes shear deformation and rotary inertia effects. In the present paper, symmetric laminated plates are considered and there will be no coupling between 165
166
I. LEEand M. H.
CHO
Fig. 1. Coordinate systems of composite panel. bending and extension of a laminate. Therefore, the nondimensional governing equations are expressed below in terms of the nondimensional displacements G, S, and $,
where
A; = A,,/E,h
(i,,j = 4, 5)
(i,.; = I. 2. 6)
0: = D,,lW’
s ill
4, = 0, (1)
(Q,,), dz ii?
i, j = 4, 5
I!2 4 =
s Ii 2
(p,,), z’ d-_ i.j=
1,2,6
I!2 (1. z’)p dz = (ph, ph3/12). @, 0 = s A2 Here, q(.u, J-,t) is the transversely distributed load; stiKnesses, D,, are bending stiffnesses; (&,,), are the transformed reduced stiffness for the k th layer of the fiber angle tIk; ki are the shear correction coefficients; p is the density of laminated plate and h is the thickness of laminated plate; and E, is Young’s modulus of any reference layer in the r-direction along the fiber or perpendicular to the fiber. The finite element model of the panel flutter is obtained using the shear deformable plate elements of Reddy [IO]. The midplane fI of the plate is subdivided into a finite number of elements, R, (e = I, 2 .). Over each element R,, the nondimensional generalized displacements (6, s,, s,.) are interpolated by expressions of the form A,, are thickness-shear
(4)
167
Finite element analysis of composite panel flutter where tii are the interpolation functions and n the number of nodes per element. Therefore, the finite element equations can be written as
Thus, the number of equations to be solved is reduced and the final finite element equation is expressed only in terms of {C}
(6) [;z;,
I::;
Z]{
ij}+u4$
[i2.~~Z,IIll{~~}=~i;n9.
(5)
where {F} are nodal forces due to the aerodynamic forces and [K] and [M] are given by
where [a is the condensed stiffness matrix and (M] = [Ml’]. The aerodynamic pressure q at a point (x, y) in supersonic flow is expressed as [l 1]
dx,
[P2] = x2{A:,[S’q + A&y[S@]} [P3] = a2{A,*,[S@j+ A&y[S”O]} [K22]= D:,[see] + D:6y([s@]
t) =
X [g
[K”] = L+{Af,[S”] + A&y([SCg] + [SC”lr) + ‘4$.4y2[P”]}
Y,
2% -(M2
_
1),,2
+;($$)$I
(M > 21’2) (7)
where q. is the dynamic pressure; M is the Mach number; V is the flow velocity; and [ the direction of the supersonic airstream on the surface of the panel (see Fig. 1). Let p = 2q,/(W - I)“2
+ [Ssq]rl + D&y2[S”“] + a2A,*,[S] [K23]= Dfzy[Sa] + Df6[Sss] + D:,y’PJ
+ D&y[SS”Jr+ dA,,[S]
[K3’] = D&[SSr] + D&y([Ssg] + [Seq]l) + D:2y2[S9”] + cc’A&[S]
where C$ is the angle between the [-direction airllow and the x-axis. The nodal forces (F} of eqn (6) are given by
&, {F) = -B*[A,l{@} ,
[M”] = [S]
[fP]
= [P]
pq$ =
Si,=
= &
I ap
z%dtdn
n,
SF=
a4
[S]
(9)
where
(p,q=O,t,tj)
s n,
--LZkl{~tj,
of
_ w4 p=E,hj [A[] = cos t#~ [Sot] + sin by [So7
tit+,dt dq,
where R, is the region of a typical element. Equations (5) can be simplified by neglecting the inertia terms of two rotations (S,, 3,) because [M*‘] and [M3’] become very small compared with [Ml’] for the thin panel (0: > 100). {??,} and {s,} can be written in terms of {@} through the Guyan reduction procedure.
HA = PI = [Ml, where /3* is the nondimensional dynamic pressure parameter. Substituting {F} in eqn (6) and assuming vc = d* e”‘, we can obtain the following equation:
[IQ + B*[A,l -
k**WIl{+*}
= {O>,
(10)
168
I. LEE and M. H. CHO
where
where [E] is an n x m normalized modal matrix obtained from the free vibration analysis and {b} contains m normal coordinates. Here, n is the dimension of the matrix in eqn (lo), and m is the number of normal modes chosen for modal analysis. Therefore, the final finite element equation can be simplified to the following equation:
k*Z= _;*Z . -P’* ;(*2
=w2a4P
h E,h’
/l* =
po.
Here, k * are the nondimensional eigenvalues, and i * are the nondimensional natural frequency parameters. The problem is reduced to that of finding the nondimensional eigenvalues and corresponding modal shapes of the system for a given value of the nondimensional dynamic pressure parameter. When p * = 0, the eigenvalues k* are real and positive, since [KJ and [M] are symmetric and positive definite. HOWever, the aerodynamic matrix [A,] is not symmetric, and hence complex eigenvalues k* are expected for p*>O. As /3* increases monotonically from zero, two of these eigenvalues will approach each other and coalesce to k,t at /I* = Bz and become complex conjugate pairs k*=kj!+ikT
(11)
for /?* > /I:,. Here, p,* corresponds
Wl* + 8*M:l- k*2Pfl*l{b) = (01, where [8* =
[El’[B PI = [A**1
WI* = Pl%W [El = [II &I* = [El%$l[El. 3. NUMERICAL RESULTS AND DISCUSSION
The following elastic properties of T300/5208 graphite/epoxy composite material are used:
to the value of
k* at which the first coalescence occurs. In general, the complex eigenvalues k* may be put in more
E, =20x
convenient
Ez=E,=1.4x
form k*’ = -g,(w/o,,)
- (w/oJ,
(12)
IO’psi
106psi
G,z = G,, = 0.8 x IO6 psi
where
Gz3 = 0.6 x lo6 psi
(M’-2)
p,V
v,* = 0.3
go = (IV* - 1)’ * phw,, ply thickness
h, = 4.92 x lo-’ in.
density p = 0.148 x IO- ’ lb-sec2/in4 Here, g, is the nondimensional aerodynamic damping parameter, and o,, a convenient frequency scale. When flutter occurs, w becomes a pure imaginary number. Hence, the flutter boundary (p,) may be obtained as a function of the aerodynamic damping g, by continuing the calculations for fl* > /I::. In the absence of aerodynamic damping, the flutter boundaries simply correspond to ,!?,*,and i * = ik *. Equation (10) can be solved directly by increasing values of fi* in many steps starting from zero. However, this procedure requires a tremendous amount of computational time. Therefore, the modal analysis is very useful in this case. It is generally seen that only a small number of eigenvectors in the lower modes gives very accurate results, since the first coalescence probably occurs among the first few natural frequencies. For this reason, {@*} can be approximated by the linear combinations of a limited number of normal modes
of the system. {I?*} Lz [E]{b},
(13)
(14)
Fig. 2. Geometry of trapezoidal laminated plate.
169
Finite element analysis of composite panel flutter
Table 1. Convergence of natural frequencies and flutter boundaries for [(O/90),], boron/epoxy rectangular clamped laminate (4 = 0, a/b = 1.0, c/b = 1.0) No. of modes
Mesh 2x2 3x3 4x4 5x5 6x6 7x7 8x8
11
Srinivasan and Babu [7] PIE method Series solution Thornton and Clary [ 121 $RR method FE method Exwrimental
1 24.97 23.61 23.42 23.37 23.35 23.34 23.34
Natural 2 56.71 45.77 43.24 42.65 42.44 42.35 42.30
frequency parameter I * 3 4 5 74.31 90.76 58.67 77.87 95.50 54.98 67.13 84.42 54.12 66.03 78.92 53.81 65.65 77.28 53.68 65.48 76.58 53.62 65.40 76.25
23.33 23.63
42.36 42.28
53.77 53.76
65.56 65.42
76.52 75.89
92.76 92.19
22.71 23.25 17.44
44.14 45.23 35.61
49.59 50.68 43.42
63.58 64.31 54.68
80.84 83.20 64.85
94.10 95.19 79.75
Flutter bounds k2: B:: 55.46 775.62 47.55 494.00 46.43 467.22 46.55 466.29 46.72 468.71 46.83 470.29 46.89 471.16
6 111.60 100.27 95.02 93.44 92.78 92.46
46.09 47.19
446.36 474.60
-
-
-
t IE method: integral equation method. Series solution: solution of double series of beam characteristic functions. $ RR method: Rayleigh-Ritz method. FE method: finite element method. Table 2. Comparison of natural frequencies and flutter boundaries with [7] for [(90/O),], boron/epoxy trapezoidal plate (I#I= 0, a/b = 0.75, c/b = 0.2) Reference
Natural frequency parameter 1* 3 4 5 105.20 114.84 166.63 180.71 105.51 118.85 168.48 187.27
1
Present method Srinivasan and Babu 171
2
58.14 58.14
The stiffness coefficients associated with A, (i,j = 4, 5) are evaluated using a reduced integration technique to circumvent the computational difficulty (i.e. to avoid the shear locking). A nine-node isoparametric, shear deformable plate element is used in this analysis. The geometry of the trapezoidal laminated plate considered in this paper is shown in Fig. 2. Convergence study
Table 1 shows the convergence characteristics of the solution for a typical [(O/90),], boron/epoxy rectangular clamped laminate (c/b = 1.0, a/b = 1.0). The following material properties of boron/epoxy composite are used for the convergence test: E, = clomped
;
191.28 199.73
94.78 96.76
999.07 1007.80
31.0 x lo6 psi, E2 = 2.7 x lo6 psi, Cl1 = 0.75 x lo6 psi, vIZ= 0.28 and p = 0.192 x 10m3lb-sec2/in4. In Table 1, the mesh size was varied from 2 x 2 to 8 x 8, keeping the number of modes constant at 11. It may be observed that the convergence of flutter boundaries is not monotonic, which is in contrast to the self-adjoint problem (free vibration) where the convergence of eigenvalues is consistently monotonic from one side. However, there is a good convergence of solutions for 7 x 7 or 8 x 8 mesh. The present results agree well with those in the available references, as can be seen from Table 1. The effect of number of modes on the flutter boundaries for fixed mesh size is examined in this analysis. There is a good convergence for 10 or
, C$= 0
c/b=0.4
Flutter bounds k,: B,:
6
simply
: c/b=0.4
supported
,@ =70
100.
:I:
F---_=i
80 -
*:I,““;,o”, 0
200
, , q 400
P’
600
800
,“;
mode
4
20 -
1 -
[45/*45/01, 0
I/.
0
I
I 200
I
1
,
400
,
600
]
,
,
10
ok; - -10
600
P’
Fig. 3. Coalescence of eigenvalues for composite trapezoidal plate (graphite/epoxy; a/b = 1.O).
-20
1. LEE and M. H. CHO
170
flutter
mode
(real
mode
2
flutter
shape
part)
mode
(imaginary
shape part)
I and 2 and flutter mode shape for [O/*45/0], graphite/epoxy trapezoidal laminate (clamped: c,+= 0, a/b = 1.0, c/b = 0.4).
Fig. 4. The evolution
of modes
11 modes. Hence, for further analysis, a mesh size of 7 x 7 and 11 modes have been chosen. Table 2 shows
that the present results agree well with those from [7] for a trapezoidal cross-ply panel. Numerical application
For graphite/epoxy trapezoidal laminated plates, flutter boundaries are calculated when the length-tothickness ratio (a/h) is 200, with the panel length (a) in the x-direction held fixed. Figure 3 shows the coalescence of eigenvalues for [O/+45/O], and [45/ +45/O], trapezoidal laminates. For [O/+45/0], and [45! *45/O], panels, the first coalescences occur at the merging points of modes 1 and 2 and of modes 3 and 4, respectively. In Fig. 4. the evolution pattern of modes 1 and 2 and the flutter mode shape with the increase of dynamic pressure up to /I* = 900 for a [O/&-45/O],panel are depicted. Two different modes (modes 1 and 2) become a single flutter mode at B: and complex conjugate pairs for p * > p:.
This shows that the lateral deflection is concentrated near the rear of the panel. Figures 5 and 6 illustrate the effect of the fiber angle for [B/+45/0], trapezoidal laminates with simply supported and clamped edges (4 = 0). The flutter boundary (fix) has the maximum value at 0 = 0 and decreases gradually as the fiber angle (0) increases. Flutter boundaries decrease as the lower to the upper side-length ratio (c/h) increases. Flutter boundaries for laminates with clamped boundaries are higher than those with simply supported boundaries. as we expect. Figures 7 and 8 show the effect of the flow direction for [Q/_t 45/O],) (when 0 = 45 and 90, respectively) laminates with simply supported edges. The maximum value of flutter boundaries appears in the range of C$= 40-70 at the [45/+45/O], laminate, and in the range of 4 = 80-90 at the [90/+45/O], laminate. This illustrates that the maximum flutter boundaries appear in the vicinity of C$= 0 (i.e. when the flow direction is consistent with the fiber angle).
50
k’,,
4o 30
01
0
b:c/b=O4
o.c/b=O.B
q:c,‘b=0.6
‘I : c/b=
15
1
/
I
30
45
60
Fiber
Fig.
5. Effect
of fiber
Angle
I 1 0
I 75
90
0
orientation on flutter boundaries trapezoidal panels (simply
‘“i0,--__! 15
30
Fiber
for [0/*45/O], supported).
45
60
Angle
graphite/epoxy
75
90
8
laminated
171
Finite element analysis of composite panel flutter
60
600
k*cr 40
Fiber
Fig. 6. Effect of fiber orientation
:oc
k*
CT
on flutter boundaries for [0/*45/O], trapezoidal panels (clamped).
[
1
1200
Angle
@
graphite/epoxy
laminated
[
-I
50 40
Flow
Angle
9,
Flow
Angle
@
Fig. 7. Effect of flow direction on flutter boundaries for [45/k 45/O], laminated trapezoidal panels (simply supported; a/b = 1.0).
600 40 400
0
0
10
20
30
Flow
40
50
Angle
60
70
60
90
@
0
0
10
20
30
40
Flow
50
60
Angle
70
60
90
@
Fig. 8. Effect of flow direction on flutter boundaries for [90/+45/O], laminated trapezoidal panels (simply supported; a/b = 1 .O). Finally, Figs 9 and 10 show the flutter boundary (j?:) as a function of the aerodynamic damping g, for [60/ _+45/O], and [90/&45/O], laminates, respectively. Generally, as the aerodynamic damping increases, the flutter boundaries increase. If the actual flight conditions are known, these graphs allow one to define the flutter boundaries including the effect of the aerodynamic damping. From these figures, we see that the boundary condition affects the flutter boundary
significantly. The flutter dynamic pressure for a clamped boundary is much higher than that for a simply supported boundary. 4. CONCLUSIONS
The application of the finite element method based on the shear deformable theory to supersonic flutter problems has been presented. The results of the
172
I. LEE and
M. H.
CHO [ 60/*45/O]s
[60/*45/O]s
1200
02 /
0
20
Fig. 9. Effect of aerodynamic
40
60
damping
0
on flutter dynamic
pressure
20
00
60
c/h (Q, = 0, o/h = 1.O)
with various
1600
1200
I”, p;
800
0
fl;
800
400 simply
supported
[90/145/O], i-
boundary
(mped
bou,ndc;y
/
,
1 : Fig.
10. Effect
of aerodynamic
damping on flutter dynamic (c/h = 0.4, a/b = 1.0).
present method agree well with those in the available references. It is found that the influence of the panel length ratio, the flow direction, the fiber orientation and boundary conditions greatly affect the flutter boundaries of trapezoidal laminated plates. For [e/f 45/O], trapezoidal laminates, the maximum values of the flutter dynamic pressure appear near the point at which the fiber angle (0) is consistent with the flow angle (4), and the flutter boundaries increase as the ratio of the lower side length to the upper side length (c/b) decreases. The flutter dynamic pressure for laminates with a clamped boundary is much higher than that for those with a simply supported boundary.
REFERENCES
1. Y. C. Fung, Some recent contributions to panel flutter research. AIAA Jnl 1, 898-909 (1963). 2. E. H. Dowell, Panel flutter: a review of the aeroelastic stability of plates and shells. AIAA Jnl 8, 386-399 (1970).
pressure
20 with
40
60
various
How
angles
3. M. D. Olson, Finite element applied to panel flutter. AIAA Jnl 5, 2261-2270 (1967). 4. M. D. Olson, Some flutter solution using finite elements. AIAA Jnl8, 747-152 (1970). Finite element 5. G. Sander, C. Bon and M. Geradin, analysis of supersonic panel flutter. In!. J. Numer. Mefh. Engng I, 379-394 (1973). orthotropic 6. D. J. Ketter, Flutter of flat, rectangular, panels. AIAA Jnl 8, 116-124 (1967). 7. R. S. Srinivasan and B. J. Babu, Free vibration and flutter of laminated quadrilateral plates. Compur. Struct. 27, 2977304 (1987). 8. C. Yang, C. H. Norris and Y. Stavsky, Elastic wave propagation in heterogeneous plates. Int. J. Solids Strucf. 2, 665-684 (1966). 9. R. D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates. J. uppl. Mech. 18,31. 38 (1951). element for the 10. J. N. Reddy, A penalty plate-bending analysis of laminated anisotropic composite plates. Ini. J. Numer. Meth. Engng 15, 118771206 (1980). Principles of Aero11. H. Ashley and R. L. Bisplinghoff, elasticity. Dover, New York (1975). and R. R. Clary, Correlation study of 12. E. A. Thornton finite element modeling for vibrations of composite material panels. In Composite Materials: Testing and Design, ASTM STP 546, pp. 111-129. American Society for Testing Materials (1974).