- Email: [email protected]
Flutter analysis for control surface of launch vehicle with dynamic stiffness
Compurers
Pergamon
00457949(95)00425-4
& Srruclures Vol. 60. No. 4. pp. 593-599, 1996 Copyright 0 1996 Elsev~er Science Ltd Printed in Great Britain. All nghts reserved 0045-7949/96 %I 5.00 + 0.00
FLUTTER ANALYSIS FOR CONTROL SURFACE OF LAUNCH VEHICLE WITH DYNAMIC STIFFNESS Seung-Kil Paek and In Lee? Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, 373-l Kusong-dong, Yusong-gu, Taejon 305-701, Korea (Received 20 April 1995)
Abstract-The root locus and iterative V-g method have been applied to analyze the flutter for a control surface of a launch vehicle with control actuators. The actuator is considered as a spring with dynamic stiffness. The fictitious mass method is adopted for an efficient modal flutter analysis. The methods arc applied to the flutter analysis of a control wing with a pneumatic actuator at a rotating axis in the supersonic region. The effect of the sweep angle on the flutter characteristics of the wing with dynamic stiffness is investigated and is compared with that of the wings with several values of static stiffness. Copyright Q 1996 Elsevier Science-Ltd.
INTRODUCTION The flutter analysis for the launch vehicle control surface to which control actuators are connected needs an approach different from those of conventional aircraft wings because the dynamic characteristics of the wing and the actuator must be considered simultaneously. Barai and Durvasula [I] had made a flutter analysis for the small aspect ratio control surfaces of launch vehicle, assuming the center of the wing root to be fixed. However, most studies are those for a fixed elastic wing and few works are devoted to the analyses considering the actuators. Actuators must have proper stiffness to negate external disturbance in order to satisfy the requirements for accuracy and stability imposed on its control and stability loop, and also to meet the requirements for the resonance property of the control wing including its actuator [2]. That stiffness is called dynamic stiffness, since it is given as a function of the Laplace variable s or frequency. The concept of dynamic stiffness gives an idea that the actuator can be regarded as a spring with variable stiffness in the flutter analysis of the system where the control surface and actuator are coupled. The idea can be easily applied to the system with dynamic stiffness, since the classical flutter analysis is made in the frequency or Laplace domain. To obtain flutter solutions, high-order eigenvalue problems need to be solved repeatedly. This requires us to consume much computing time and large computer memory storage, so we must undertake a modal approach to reduce those factors. However, care must be taken in the selection of vibration modes for a modal approach. Since the local system par-
t Author addressed.
to
whom
all correspondence
should
be
ameters vary continuously in the system which has dynamic stiffness, the selected vibration modes have to be able to accommodate these local structural changes. Karpel and Wieseman [3] have suggested a fictitious mass method, obtaining vibration modes that give very accurate solutions in cases where large local structural changes occur. Karpel [4] had used this method as a technique of component mode synthesis. Livne [5] applied it to the computation of eigen derivatives of control-augmented structures. In this study, we will introduce a flutter analysis technique for the aeroelastic system with dynamic stiffness and explain how to select proper modes for efficient calculation. These are applied to the control surface with an actuator. The results are compared with those of the case where the actuator is assumed to be a spring with constant stiffness.
ANALYSIS
Equations
cf motion
The general aeroelastic written as
equation
of motion
can be
where [GM], [CC], [GK], p, [Q(k)] and {q j are the generalized mass, damping, stiffness matrices, dynamic pressure, generalized oscillatory aerodynamic influence coefficient matrix and generalized coordinates, respectively. The unsteady aerodynamics by lifting surface theories, such as the doublet point method (DPM) or the doublet lattice method (DLM) are generally functions of the Mach number and reduced frequency, k = bw/V, where b, w and V are the reference length, oscillatory frequency, and free stream velocity, respectively.
Seung-Kil Paek and In Lee
594
[Q(k)1g [&(jk)l= [&I + tA,l(jk) + [.41(jk)2
Root locus method
The dynamic stiffness of an actuator can be expressed as a transfer function between the external forces and the displacement of its actuating part. This transfer function can be written as a rational function of s. Neglecting the damping term, the Laplace transform of eqn (1) gives the following equation:
6 L%l(jk) +,T3(jk)+8_
r&s)1= t44 + M(;)s
[GKI= F&I + kbITWA~W~I,
+ M(;)i
Substituting eqns (4) and (6) into eqn (2) and multiplying the least common multiple of denominators on both sides of eqn (2) yields ([F,]P + [F,_
,]srn-’ +
‘.. + ]F,I){4} = 0,
(7)
where m is the highest order of s. If we define {X} as
(4) eqn (7) results in a complex eigenvalue problem:
where [G&I and [@,I are the stiffness matrix of the nominal structure and the modal displacement matrix for degrees of freedom included in set B. To obtain [Q(s)], the generalized oscillatory aerodynamic influence coefficient matrices, [Q(k)], are
0
(j)
(6)
where NB is the number of degrees of freedom with dynamic stiffness, [K,,(s)] is the rational function matrix for dynamic stiffness of actuators, Xi is the Laplace transform of the elements of set B and 3 (i = 1 _ N,) is the concentrated forces applied to set B, respectively. Assuming that the structural behavior can be represented as a linear combination of vibration modes, the stiffness matrix in eqn (2) can be written as
VI
’
(3)
1
I
2
where j = fl and p, _ 2 is an arbitrary constant. As described in Ref. [6] the form of eqn (5) permits an approximation of the time delays inherent in unsteady aerodynamics subject to the following requirements: complex conjugate symmetry and denominator roots in the left-hand plane. Thus, good approximating coefficients, [Ai] (i = O-6) in eqn (5) are evaluated by a least-squares curve fit with the tabulated [Q(k)]s. Then [Q(s)] is obtained by inserting jk = sb/V into eqn (5);
where (41 is the Laplace transform of {q} and [Q(s)] is the unsteady aerodynamic influence coefficient matrix in the Laplace (s) domain, respectively. Let a set of degrees of freedom with dynamic stiffness be called B. Then the following relation can be obtained:
[-F,‘F,_,]
I
six}=L%,I{X), where
[-F,‘F,_,]
.‘.
0
0 -
computed at various tabulated reduced frequency (k) values and approximated by the rational interpolation function in the Laplace domain. There are many approximation techniques-rational function approximation (RFA) is one such technique. In this study, Roger’s method [6] is used. Roger suggested the following approximating function of polynomials and simple poles for unsteady aerodynamic transfer function:
(9)
[-F,‘F,l
[-K’FJ
...
0
0
...
0
0
...
il
The eigenvalues of [A,] are the roots of the characteristic flutter equation. Since [F,](i = 0 - m) in [A,] are functions of dynamic pressure for a constant Mach number, the loci of roots as a function of dynamic pressure can be constructed. These loci correspond to the variation in the eigenvalues of each flexible mode as dynamic pressure is varied. The point where a locus meets the imaginary axis of the s-plane is the flutter point.
Flutter analysis of launch vehicle Iterative
V-g method
Neglecting the damping term in eqn (1) and assuming that jq) = {4}eJo',eqn (1) becomes
([GKI- 4GMl){~l =~K?(k)1{4).
(10)
595
structure where large local structural changes occur. Let the set of degrees of freedom where local structural changes occur be B and the complimentary set of B be I. The equation of motion of undamped free vibration system can be written as follows:
With the use of an appropriate transformation, eqn (10) becomes a typical eigenvalue problem:
K[GKl- @@)I)- 1.*[GKlIid = 0,
(11) where M and K are the finite element mass and stiffness matrices, respectively, which can be divided into the sub matrices for I and B. When the fictitious masses are added to the B set, eqn (14) then becomes
where
@(~)I =
$nk)l,
+l+jg
02 ’
and g is the structural damping. The stiffness of the actuator is assumed to be the stiffness at the flutter frequency. Assuming that the flutter frequency is w,,,, the stiffness matrix can be written as where MF is a diagonal matrix whose diagonal el-
(12) WI = [GKI,+ [@‘B1~~~A.i~est)l Pd. ements are concentrated fictitious masses. Assuming Substituting eqn (12) into eqn (1 I), the characteristic flutter equation becomes
IWKI - [!&)I) - J.*W&I
[@I=
From the complex eigenvalue, the oscillatory frequency, structural damping and air speed V can be computed: Im(l*) U=&’
that NM modes are computed from the mass-added system, the vibration modal matrix is
[
I’ 6’
I
(NJ+ Nd x NM,
(16)
where N, is the number of the element of the I set and NB is that of the B set. The generalized matrices can be obtained as follows:
I&@.
g=Re(i*)’
k
When the sign of g is negative, the structure is dynamically stable. Otherwise, the structure is unstable and hence g = 0 is the flutter boundary. The flutter speed can be computed from V-g diagram and the corresponding frequency is the flutter frequency. The procedure is repeated until the difference between the computed flutter frequency and the assumed value satisfies the pre-defined tolerance. FICTITIOUS MASS METHOD
A temporary mass matrix is constructed by adding very large fictitious masses to the diagonal elements of a nominal mass matrix, whose corresponding degrees of freedom are expected to have local structural changes. From the undamped free vibration analysis with this mass matrix, the finite element vibration modes are obtained. The obtained modes have large distortions in the positions where the fictitious masses are loaded. These distortions cause the modes to clearly represent the behavior of the
where [GK] is a diagonal matrix but generally [GM] is not. To obtain the vibration modes that give a good representation for the behavior of a locally varying system, the reduced-order bases selection techniques, which are the component mode synthesis methods, can be also used as well as the fictitious mass method. In Ref. [5], Livne indicated that the results when the fictitious mass method was applied to his controlaugmented structures were the same as those when the Craig-Bampton method [7] was applied. The Craig-Bampton method, where each substructure is represented by its constraint mode set and fixed interface normal mode set. is representative of one of the component mode synthesis methods and is wellknown for its good accuracy. The modes by the fictitious mass method consist of a very low frequency mode group and the other mode group. As the fictitious mass increases, the low frequency mode
Seung-Kil Paek and In Lee
596
6
0.1
-
--_--
i;
;&
%.-.
t
\-
0.01
.Y=_,___ __~___-------o .*------=--___________ I
I
I
10
12
14
0.001 4
6
6
16
Number of Modes Used
Fig. 1. The shape. and dimension of the sample model (mm). group spans the same subspace as the constraint mode set of the Craig-Bampton method, and the other group approaches the fixed interface normal mode set. Weissenburger [8] observed the effect of local structure changes on linear vibration properties. He reported that concentrated masses show the above-mentioned effect. It must be noted that only one natural mode analysis that is very intuitive, has to be undertaken in the fictitious mass method, while a static analysis and a natural mode analysis for different boundary conditions have to be done in the Craig-Bampton method. This is the most attractive
0.0
-0.5
x 3
-1.0
Dato points
0
co’
r
2%:,” -1.5
Fig. 3. Error vs number of modes used (fictitious mass = 10,000 x nominal value, reference Mach number = 1.4). side when engaged in the aeroelastic analysis, like in the present case where large local structural changes due to dynamic stiffness occur in the dynamic behavior.
NUMERICALRESULTSANDDISCUSSION Problem description
The previously mentioned methods are applied to a simple plate wing with a control actuator as shown in Fig. 1. The wing will be denoted from now on as ‘sample model’. The dimensions of wing root chord length, root thickness, wing tip chord length. tip thickness and half span length are 150, 10, 78, 4 and 156 mm, respectively. Mass density, Young’s modulus and shear modulus of the material are 2713 kg mm3, 72.4 GPa and 26.2 GPa, respectively. Air density and air speed are assumed to be 1.225 kg m-3 and 340.3 m s-l, respectively. The structure is modeled as 4 x 4 eight-node shear deformable elements [9]. The actuating shaft of a pneumatic control actuator (PCA) is connected to the center of the wing root. The dynamic stiffness of PCA is represented as
b-1. ,.-:..
-2.0
T-
-1.0
-1.2
-0.6
Re(dd
$p(s)
-0.6
(m2)
gwlz,
0.0
W
-
Approximation function
1)(S/G(S/P,
I)(s/z,-
l(S/Pz--
l)G/z4-
1)
1)
(19) '
-0.5 0 loo
7
10
-l.O-
^
& k -1.5-
-
s
1
0.1
-
Flutter speed Flutterreducedfrequency
t
h
k-1.00 -2.0
, 0.7
0.6
0.9
1.0
1.1
Multiplesof fictitiousmass of DOF Fig. 2. Aerodynamic curve fit for (a) $,, at M = 1.4,(b) oz2 at M= 1.4.
B
Fig. 4. Error vs multiples of fictitious mass (with six modes used, reference Mach number = 1.4).
Flutter analysis of launch vehicle where T, is the applied torque to the actuator and 0 is the rotation angle of the actuator. z, (i = l-4) are - 1353.92 (215.48 Hz), -713.53 (11356Hz) and - 60.14 + j143.65 (24.79 Hz, c = 0.386), respectively. p, (i=l,2) are -666.67 (106.104 Hz), -34.11 (5.43 Hz), respectively. Ver$cation qf the Jictitious mass method To verify the effectiveness of the fictitious mass method, a flutter analysis is performed for the sample model. The fictitious mass method is used to obtain vibration modes for modal approach, where the rotation degrees of freedom, 0 in Fig. 1, has dynamic stiffness by the connected actuator, and fictitious masses are added to the corresponding diagonai element of finite element mass matrix which is called the nominal value from now on. The amount of fictitious mass can be represented by the multitude of the nominal value. Supersonic DPM (lo] is used to calculate unsteady aerodynamics in supersonic flow region. The interpolation between aerodynamic grid and structural grid is performed by surface spline [I 11. The aerodynamic mesh for DPM consists of 7 x 30 elements. The unsteady aerodynamic force is approximated by eqn (6) for the root locus method. /I,s (i = 14) are obtained by the simplex direct search method [12] so that the errors due to approximation can be minimized. The used aerodynamic data and
597
the approximation functions for o,, and $,, at A4 = 1.4 are compared in Fig. 2, which shows the approximation is accurate enough to use. d,, and & are selected because they are especially important in flutter analyses. The flutter solutions using the vibration modes computed by the fictitious mass method are compared with those by direct full-size calculation. The iterative V-g method is used for the full-size calculation. The percentage error E is defined as E
Smodrl - sciir.ect
=
--
x 100 (O/o),
&WC,
(20)
where Smodalis the result by modal approach and Sdirec, is the result by direct full-size calculation. Figure 3 shows how much the errors decrease as the number of vibration modes increases when the magnitude of fictitious mass used to obtain the modes is 10,000 times larger than the nominal value. The word FM in the figure indicates that it is the results by fictitious mass method. In Fig. 3, there are also the results by vibration modes of the nominal structure. Figure 4 shows how the results change as the magnitude of fictitious mass increases. In this study, six modes have been used. It is observed that it is difficult (a) 4.5
-
2nd (316.3 Hz)
b
3rd (1212.2
Direct calculation Q
Modal approach (FM)
Hz)
1.5
2.0
2.5
3.0
I!& Reference 4th (1353.0
Hz)
5th (2150.7
Hz)
6th (2651.7
Mach
Number
Hz)
0.‘) 285
(b)
I3 1 st (0 Hz)
4th (1015.7 Hz)
,-------
-
Direct calwlation Modal approach (FM)
b 2nd (307.3 Hz)
5th (1216.6 Hz)
3rd (998.4 Hz)
6th (2152.7 Hz)
Fig. 5. The nodal lines of the (a) nominal vibration modes. (b) vibration modes by the fictitious mass method.
Reference
Mach
Number
Fig. 6. (a) Flutter speed and (b) flutter frequency of the sample model with actuator.
Seung-Kil Paek and In Lee
598
to obtain an accurate result with the modes from the nominal structure even if the number of modes increases, and that sufficiently accurate results can be obtained if the magnitude of fictitious mass is 10,000 times larger than the nominal value and the number of modes is six. Figure 5 shows the nodal lines of the nominal modes and the modes by the fictitious mass method. The results by the modal approach with these modes are compared with those by direct fullsize calculation as shown in Fig. 6. It is observed that both results are nearly the same. From the above results, the fictitious mass method can be applied effectively to obtain accurate results while reducing the size of the aeroelastic system with dynamic stiffness. The general wings have the property of coupled mode flutter, where a bending mode and a torsion mode coalesce as the air speed increases. When the actuator is considered however, it is observed that the mode coupling does not appear and a unique mode becomes unstable as the air speed increases. Parametric studies
The effects of wing sweep back angle on the flutter speed are investigated for the present model. The
sweep back angle at mid chord changes from -30” to 30.0”, while the wing chord length and span remain constant. The result is shown in Fig. 7, which indicates that there exists a specific angle where the flutter speed has the maximum value and that the angle is about 0” in this case. When the sweep back angle changes to negative, that is, when the sweep forward angle becomes positive, the flutter mode changes. The flutter speed decreases as the sweep forward angle increases. The flutter frequency is high for the sweep back wing and low for the sweep forward wing. The flutter modes are different in these wings. It must be also noted that divergence never occurs though the wings are swept forward, while conventional forward swept wings generally undergo divergence. Flutter analyses are performed in the case where the actuator is replaced by some springs of appropriate static stiffnesses. The results are compared for four different values of stiffness in Fig. 8: 0.2k,; k,; 2k, and lOk,, where k, is referred to eqn (19). It demonstrates that sweep forward causes divergence and that when the stiffness is k, the flutter speed has the minimum value about the sweep angle of 0,. unlike in the previous case. It also shows that the e-
(a) .5-
4-
3 2’.
2 A’
compatibility
line 1 -t compatibility
I
0 -30
-20
-10
0
10
20
30
-30
Sweep angle at mid chord,0 (degree)
500
(b) 300 -
77
d G 250
Lz
-10
0
10
20
3c
T(b)
400
G
G :s
-20
Sweep angle at mid chord,@ (degree)
1-1
z5
line
is E
200 -
300
G 150 -
&
200
‘i 5 i;
2 i;
100 -
100 50 -
0
0 , -30
-20
-10
0
10
Sweep angle at mid chord,
20
30
(degree)
Fig. 7. Efkcts of sweep angle on the (a) flutter speed and (b) flutter frequency of the sample model with actuator (reference Mach number = 1.4).
1
-30
-20
-10
0
10
20
30
Sweep angle at mid chard,0 (degree)
Fig. 8. Effects of sweep angle on the (a) Butter speed and (b) flutter frequency of the sample model with torsional spring (reference Mach number = 1.4).
Flutter
analysis
(a) 5
compatibility
Stiffness
line
Stiffness
Ratio,r(k/ko)
Fig. 9. Effects of the torsional spring magnitude on the (a) flutter speed and (b) flutter frequency for various sweep angles (reference Mach number = 1.4).
sweep angle vs flutter
speed
graph
changes
substan-
The flutter frequency increases as the stiffness increases. Figure 9 shows how the flutter speed changes as the stiffness changes at a particular sweep angle: 0”; 10” and 20”. There exists a particular stiffness magnitude where flutter speed becomes a minimum at each configuration. the
599
the dynamic properties of the control actuator for the elastic launch vehicle wing. It is shown that the fictitious mass method can be applied effectively in obtaining accurate results efficiently in the flutter analysis for the aeroelastic system with an actuator. It is observed that there exists a sweep back angle at which the flutter speed becomes the maximum value for the actuator used in this study. For the case with static stiffness, there exists a sweep angle where the flutter speed becomes a minimum. It must be noted that the cases with static stiffness show a significantly different trend from those with an actuator and that an actuator cannot be replaced by a static spring with an appropriate stiffness. The dynamic characteristics of actuators must be properly considered for an accurate flutter prediction. REFERENCES
5oT-----l
as
vehicle
Ratio.r(k/kg)
(b)
tially
of launch
magnitude
of
stiffness
changes.
CONCLUSIONS
In this paper, the root locus method and iterative V-g method are used for flutter analysis considering
1. A. Barai and S. Durvasula. Flutter of shaft-supported low aspect-ratio control surfaces. J. Aeronaut. Sot. India 42, 19-26 (1990). 2. A. B. Kondrat’ev and V. A. Chashchin, On the stiffness of a flight vehicle pneumatic control actuator. IZO. Vi/Z. Auiarsion. Tekh. 33, 7--10 (1990). of 3. M. Karpel and C. D. Wieseman, Time simulation flutter with large stiffness changes. In: Proc. 33th AIAA/ASME/ASCEIAHS/ACS Structural Dynamics and Materials Corzference, paper 92-2394, Dallas, TX, pp. 1878-1886 (1992). 4. M. Karpel, Efficient vibration mode analysis of aircraft with multiple external store configuration. J. Aircr. 25, 747-751 (1988). 5. E. Livne, Accurate calculation of control-augmented structural eigenvalue sensitivities using reduced-order models. AIAA .I. 27, 947-954 (1989). math modeling methods for 6. K. L. Roger, Airplane active control design, AGARD-CP-228. pp. 4-14-l I (1977) 7. R. R. Craig Jr and M. C. C. Bampton, Coupling of substructures for dynamic analysis. AIAA. J. 6, 13131319 (1968). 8. J. T. Weissenburger, Effect of local modification on the vibration characteristics of linear systems. Trans. ASME Ser. E 35, 327.-332 (1968). 9. 1. Lee and J. J. Lee, Vibration analysis of composite plate wings. Compur. S/ruci. 37, 1077.-1085 (1990). 10. T. Ueda and E. H. Dowell, Doublet-point method for supersonic unsteady lifting surfaces. AIAA J. 22, 179.-186 (1984). II. R. L. Harder and R. N. Desmarais, Interpolation using surface splines. J. Aircr. 9, 189-191 (1972). 12. J. A. Nelder and R. Mead, A simplex method for function minimization. Compur. J. 7, 279-285 (1965).
Pergamon
00457949(95)00425-4
& Srruclures Vol. 60. No. 4. pp. 593-599, 1996 Copyright 0 1996 Elsev~er Science Ltd Printed in Great Britain. All nghts reserved 0045-7949/96 %I 5.00 + 0.00
FLUTTER ANALYSIS FOR CONTROL SURFACE OF LAUNCH VEHICLE WITH DYNAMIC STIFFNESS Seung-Kil Paek and In Lee? Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, 373-l Kusong-dong, Yusong-gu, Taejon 305-701, Korea (Received 20 April 1995)
Abstract-The root locus and iterative V-g method have been applied to analyze the flutter for a control surface of a launch vehicle with control actuators. The actuator is considered as a spring with dynamic stiffness. The fictitious mass method is adopted for an efficient modal flutter analysis. The methods arc applied to the flutter analysis of a control wing with a pneumatic actuator at a rotating axis in the supersonic region. The effect of the sweep angle on the flutter characteristics of the wing with dynamic stiffness is investigated and is compared with that of the wings with several values of static stiffness. Copyright Q 1996 Elsevier Science-Ltd.
INTRODUCTION The flutter analysis for the launch vehicle control surface to which control actuators are connected needs an approach different from those of conventional aircraft wings because the dynamic characteristics of the wing and the actuator must be considered simultaneously. Barai and Durvasula [I] had made a flutter analysis for the small aspect ratio control surfaces of launch vehicle, assuming the center of the wing root to be fixed. However, most studies are those for a fixed elastic wing and few works are devoted to the analyses considering the actuators. Actuators must have proper stiffness to negate external disturbance in order to satisfy the requirements for accuracy and stability imposed on its control and stability loop, and also to meet the requirements for the resonance property of the control wing including its actuator [2]. That stiffness is called dynamic stiffness, since it is given as a function of the Laplace variable s or frequency. The concept of dynamic stiffness gives an idea that the actuator can be regarded as a spring with variable stiffness in the flutter analysis of the system where the control surface and actuator are coupled. The idea can be easily applied to the system with dynamic stiffness, since the classical flutter analysis is made in the frequency or Laplace domain. To obtain flutter solutions, high-order eigenvalue problems need to be solved repeatedly. This requires us to consume much computing time and large computer memory storage, so we must undertake a modal approach to reduce those factors. However, care must be taken in the selection of vibration modes for a modal approach. Since the local system par-
t Author addressed.
to
whom
all correspondence
should
be
ameters vary continuously in the system which has dynamic stiffness, the selected vibration modes have to be able to accommodate these local structural changes. Karpel and Wieseman [3] have suggested a fictitious mass method, obtaining vibration modes that give very accurate solutions in cases where large local structural changes occur. Karpel [4] had used this method as a technique of component mode synthesis. Livne [5] applied it to the computation of eigen derivatives of control-augmented structures. In this study, we will introduce a flutter analysis technique for the aeroelastic system with dynamic stiffness and explain how to select proper modes for efficient calculation. These are applied to the control surface with an actuator. The results are compared with those of the case where the actuator is assumed to be a spring with constant stiffness.
ANALYSIS
Equations
cf motion
The general aeroelastic written as
equation
of motion
can be
where [GM], [CC], [GK], p, [Q(k)] and {q j are the generalized mass, damping, stiffness matrices, dynamic pressure, generalized oscillatory aerodynamic influence coefficient matrix and generalized coordinates, respectively. The unsteady aerodynamics by lifting surface theories, such as the doublet point method (DPM) or the doublet lattice method (DLM) are generally functions of the Mach number and reduced frequency, k = bw/V, where b, w and V are the reference length, oscillatory frequency, and free stream velocity, respectively.
Seung-Kil Paek and In Lee
594
[Q(k)1g [&(jk)l= [&I + tA,l(jk) + [.41(jk)2
Root locus method
The dynamic stiffness of an actuator can be expressed as a transfer function between the external forces and the displacement of its actuating part. This transfer function can be written as a rational function of s. Neglecting the damping term, the Laplace transform of eqn (1) gives the following equation:
6 L%l(jk) +,T3(jk)+8_
r&s)1= t44 + M(;)s
[GKI= F&I + kbITWA~W~I,
+ M(;)i
Substituting eqns (4) and (6) into eqn (2) and multiplying the least common multiple of denominators on both sides of eqn (2) yields ([F,]P + [F,_
,]srn-’ +
‘.. + ]F,I){4} = 0,
(7)
where m is the highest order of s. If we define {X} as
(4) eqn (7) results in a complex eigenvalue problem:
where [G&I and [@,I are the stiffness matrix of the nominal structure and the modal displacement matrix for degrees of freedom included in set B. To obtain [Q(s)], the generalized oscillatory aerodynamic influence coefficient matrices, [Q(k)], are
0
(j)
(6)
where NB is the number of degrees of freedom with dynamic stiffness, [K,,(s)] is the rational function matrix for dynamic stiffness of actuators, Xi is the Laplace transform of the elements of set B and 3 (i = 1 _ N,) is the concentrated forces applied to set B, respectively. Assuming that the structural behavior can be represented as a linear combination of vibration modes, the stiffness matrix in eqn (2) can be written as
VI
’
(3)
1
I
2
where j = fl and p, _ 2 is an arbitrary constant. As described in Ref. [6] the form of eqn (5) permits an approximation of the time delays inherent in unsteady aerodynamics subject to the following requirements: complex conjugate symmetry and denominator roots in the left-hand plane. Thus, good approximating coefficients, [Ai] (i = O-6) in eqn (5) are evaluated by a least-squares curve fit with the tabulated [Q(k)]s. Then [Q(s)] is obtained by inserting jk = sb/V into eqn (5);
where (41 is the Laplace transform of {q} and [Q(s)] is the unsteady aerodynamic influence coefficient matrix in the Laplace (s) domain, respectively. Let a set of degrees of freedom with dynamic stiffness be called B. Then the following relation can be obtained:
[-F,‘F,_,]
I
six}=L%,I{X), where
[-F,‘F,_,]
.‘.
0
0 -
computed at various tabulated reduced frequency (k) values and approximated by the rational interpolation function in the Laplace domain. There are many approximation techniques-rational function approximation (RFA) is one such technique. In this study, Roger’s method [6] is used. Roger suggested the following approximating function of polynomials and simple poles for unsteady aerodynamic transfer function:
(9)
[-F,‘F,l
[-K’FJ
...
0
0
...
0
0
...
il
The eigenvalues of [A,] are the roots of the characteristic flutter equation. Since [F,](i = 0 - m) in [A,] are functions of dynamic pressure for a constant Mach number, the loci of roots as a function of dynamic pressure can be constructed. These loci correspond to the variation in the eigenvalues of each flexible mode as dynamic pressure is varied. The point where a locus meets the imaginary axis of the s-plane is the flutter point.
Flutter analysis of launch vehicle Iterative
V-g method
Neglecting the damping term in eqn (1) and assuming that jq) = {4}eJo',eqn (1) becomes
([GKI- 4GMl){~l =~K?(k)1{4).
(10)
595
structure where large local structural changes occur. Let the set of degrees of freedom where local structural changes occur be B and the complimentary set of B be I. The equation of motion of undamped free vibration system can be written as follows:
With the use of an appropriate transformation, eqn (10) becomes a typical eigenvalue problem:
K[GKl- @@)I)- 1.*[GKlIid = 0,
(11) where M and K are the finite element mass and stiffness matrices, respectively, which can be divided into the sub matrices for I and B. When the fictitious masses are added to the B set, eqn (14) then becomes
where
@(~)I =
$nk)l,
+l+jg
02 ’
and g is the structural damping. The stiffness of the actuator is assumed to be the stiffness at the flutter frequency. Assuming that the flutter frequency is w,,,, the stiffness matrix can be written as where MF is a diagonal matrix whose diagonal el-
(12) WI = [GKI,+ [@‘B1~~~A.i~est)l Pd. ements are concentrated fictitious masses. Assuming Substituting eqn (12) into eqn (1 I), the characteristic flutter equation becomes
IWKI - [!&)I) - J.*W&I
[@I=
From the complex eigenvalue, the oscillatory frequency, structural damping and air speed V can be computed: Im(l*) U=&’
that NM modes are computed from the mass-added system, the vibration modal matrix is
[
I’ 6’
I
(NJ+ Nd x NM,
(16)
where N, is the number of the element of the I set and NB is that of the B set. The generalized matrices can be obtained as follows:
I&@.
g=Re(i*)’
k
When the sign of g is negative, the structure is dynamically stable. Otherwise, the structure is unstable and hence g = 0 is the flutter boundary. The flutter speed can be computed from V-g diagram and the corresponding frequency is the flutter frequency. The procedure is repeated until the difference between the computed flutter frequency and the assumed value satisfies the pre-defined tolerance. FICTITIOUS MASS METHOD
A temporary mass matrix is constructed by adding very large fictitious masses to the diagonal elements of a nominal mass matrix, whose corresponding degrees of freedom are expected to have local structural changes. From the undamped free vibration analysis with this mass matrix, the finite element vibration modes are obtained. The obtained modes have large distortions in the positions where the fictitious masses are loaded. These distortions cause the modes to clearly represent the behavior of the
where [GK] is a diagonal matrix but generally [GM] is not. To obtain the vibration modes that give a good representation for the behavior of a locally varying system, the reduced-order bases selection techniques, which are the component mode synthesis methods, can be also used as well as the fictitious mass method. In Ref. [5], Livne indicated that the results when the fictitious mass method was applied to his controlaugmented structures were the same as those when the Craig-Bampton method [7] was applied. The Craig-Bampton method, where each substructure is represented by its constraint mode set and fixed interface normal mode set. is representative of one of the component mode synthesis methods and is wellknown for its good accuracy. The modes by the fictitious mass method consist of a very low frequency mode group and the other mode group. As the fictitious mass increases, the low frequency mode
Seung-Kil Paek and In Lee
596
6
0.1
-
--_--
i;
;&
%.-.
t
\-
0.01
.Y=_,___ __~___-------o .*------=--___________ I
I
I
10
12
14
0.001 4
6
6
16
Number of Modes Used
Fig. 1. The shape. and dimension of the sample model (mm). group spans the same subspace as the constraint mode set of the Craig-Bampton method, and the other group approaches the fixed interface normal mode set. Weissenburger [8] observed the effect of local structure changes on linear vibration properties. He reported that concentrated masses show the above-mentioned effect. It must be noted that only one natural mode analysis that is very intuitive, has to be undertaken in the fictitious mass method, while a static analysis and a natural mode analysis for different boundary conditions have to be done in the Craig-Bampton method. This is the most attractive
0.0
-0.5
x 3
-1.0
Dato points
0
co’
r
2%:,” -1.5
Fig. 3. Error vs number of modes used (fictitious mass = 10,000 x nominal value, reference Mach number = 1.4). side when engaged in the aeroelastic analysis, like in the present case where large local structural changes due to dynamic stiffness occur in the dynamic behavior.
NUMERICALRESULTSANDDISCUSSION Problem description
The previously mentioned methods are applied to a simple plate wing with a control actuator as shown in Fig. 1. The wing will be denoted from now on as ‘sample model’. The dimensions of wing root chord length, root thickness, wing tip chord length. tip thickness and half span length are 150, 10, 78, 4 and 156 mm, respectively. Mass density, Young’s modulus and shear modulus of the material are 2713 kg mm3, 72.4 GPa and 26.2 GPa, respectively. Air density and air speed are assumed to be 1.225 kg m-3 and 340.3 m s-l, respectively. The structure is modeled as 4 x 4 eight-node shear deformable elements [9]. The actuating shaft of a pneumatic control actuator (PCA) is connected to the center of the wing root. The dynamic stiffness of PCA is represented as
b-1. ,.-:..
-2.0
T-
-1.0
-1.2
-0.6
Re(dd
$p(s)
-0.6
(m2)
gwlz,
0.0
W
-
Approximation function
1)(S/G(S/P,
I)(s/z,-
l(S/Pz--
l)G/z4-
1)
1)
(19) '
-0.5 0 loo
7
10
-l.O-
^
& k -1.5-
-
s
1
0.1
-
Flutter speed Flutterreducedfrequency
t
h
k-1.00 -2.0
, 0.7
0.6
0.9
1.0
1.1
Multiplesof fictitiousmass of DOF Fig. 2. Aerodynamic curve fit for (a) $,, at M = 1.4,(b) oz2 at M= 1.4.
B
Fig. 4. Error vs multiples of fictitious mass (with six modes used, reference Mach number = 1.4).
Flutter analysis of launch vehicle where T, is the applied torque to the actuator and 0 is the rotation angle of the actuator. z, (i = l-4) are - 1353.92 (215.48 Hz), -713.53 (11356Hz) and - 60.14 + j143.65 (24.79 Hz, c = 0.386), respectively. p, (i=l,2) are -666.67 (106.104 Hz), -34.11 (5.43 Hz), respectively. Ver$cation qf the Jictitious mass method To verify the effectiveness of the fictitious mass method, a flutter analysis is performed for the sample model. The fictitious mass method is used to obtain vibration modes for modal approach, where the rotation degrees of freedom, 0 in Fig. 1, has dynamic stiffness by the connected actuator, and fictitious masses are added to the corresponding diagonai element of finite element mass matrix which is called the nominal value from now on. The amount of fictitious mass can be represented by the multitude of the nominal value. Supersonic DPM (lo] is used to calculate unsteady aerodynamics in supersonic flow region. The interpolation between aerodynamic grid and structural grid is performed by surface spline [I 11. The aerodynamic mesh for DPM consists of 7 x 30 elements. The unsteady aerodynamic force is approximated by eqn (6) for the root locus method. /I,s (i = 14) are obtained by the simplex direct search method [12] so that the errors due to approximation can be minimized. The used aerodynamic data and
597
the approximation functions for o,, and $,, at A4 = 1.4 are compared in Fig. 2, which shows the approximation is accurate enough to use. d,, and & are selected because they are especially important in flutter analyses. The flutter solutions using the vibration modes computed by the fictitious mass method are compared with those by direct full-size calculation. The iterative V-g method is used for the full-size calculation. The percentage error E is defined as E
Smodrl - sciir.ect
=
--
x 100 (O/o),
&WC,
(20)
where Smodalis the result by modal approach and Sdirec, is the result by direct full-size calculation. Figure 3 shows how much the errors decrease as the number of vibration modes increases when the magnitude of fictitious mass used to obtain the modes is 10,000 times larger than the nominal value. The word FM in the figure indicates that it is the results by fictitious mass method. In Fig. 3, there are also the results by vibration modes of the nominal structure. Figure 4 shows how the results change as the magnitude of fictitious mass increases. In this study, six modes have been used. It is observed that it is difficult (a) 4.5
-
2nd (316.3 Hz)
b
3rd (1212.2
Direct calculation Q
Modal approach (FM)
Hz)
1.5
2.0
2.5
3.0
I!& Reference 4th (1353.0
Hz)
5th (2150.7
Hz)
6th (2651.7
Mach
Number
Hz)
0.‘) 285
(b)
I3 1 st (0 Hz)
4th (1015.7 Hz)
,-------
-
Direct calwlation Modal approach (FM)
b 2nd (307.3 Hz)
5th (1216.6 Hz)
3rd (998.4 Hz)
6th (2152.7 Hz)
Fig. 5. The nodal lines of the (a) nominal vibration modes. (b) vibration modes by the fictitious mass method.
Reference
Mach
Number
Fig. 6. (a) Flutter speed and (b) flutter frequency of the sample model with actuator.
Seung-Kil Paek and In Lee
598
to obtain an accurate result with the modes from the nominal structure even if the number of modes increases, and that sufficiently accurate results can be obtained if the magnitude of fictitious mass is 10,000 times larger than the nominal value and the number of modes is six. Figure 5 shows the nodal lines of the nominal modes and the modes by the fictitious mass method. The results by the modal approach with these modes are compared with those by direct fullsize calculation as shown in Fig. 6. It is observed that both results are nearly the same. From the above results, the fictitious mass method can be applied effectively to obtain accurate results while reducing the size of the aeroelastic system with dynamic stiffness. The general wings have the property of coupled mode flutter, where a bending mode and a torsion mode coalesce as the air speed increases. When the actuator is considered however, it is observed that the mode coupling does not appear and a unique mode becomes unstable as the air speed increases. Parametric studies
The effects of wing sweep back angle on the flutter speed are investigated for the present model. The
sweep back angle at mid chord changes from -30” to 30.0”, while the wing chord length and span remain constant. The result is shown in Fig. 7, which indicates that there exists a specific angle where the flutter speed has the maximum value and that the angle is about 0” in this case. When the sweep back angle changes to negative, that is, when the sweep forward angle becomes positive, the flutter mode changes. The flutter speed decreases as the sweep forward angle increases. The flutter frequency is high for the sweep back wing and low for the sweep forward wing. The flutter modes are different in these wings. It must be also noted that divergence never occurs though the wings are swept forward, while conventional forward swept wings generally undergo divergence. Flutter analyses are performed in the case where the actuator is replaced by some springs of appropriate static stiffnesses. The results are compared for four different values of stiffness in Fig. 8: 0.2k,; k,; 2k, and lOk,, where k, is referred to eqn (19). It demonstrates that sweep forward causes divergence and that when the stiffness is k, the flutter speed has the minimum value about the sweep angle of 0,. unlike in the previous case. It also shows that the e-
(a) .5-
4-
3 2’.
2 A’
compatibility
line 1 -t compatibility
I
0 -30
-20
-10
0
10
20
30
-30
Sweep angle at mid chord,0 (degree)
500
(b) 300 -
77
d G 250
Lz
-10
0
10
20
3c
T(b)
400
G
G :s
-20
Sweep angle at mid chord,@ (degree)
1-1
z5
line
is E
200 -
300
G 150 -
&
200
‘i 5 i;
2 i;
100 -
100 50 -
0
0 , -30
-20
-10
0
10
Sweep angle at mid chord,
20
30
(degree)
Fig. 7. Efkcts of sweep angle on the (a) flutter speed and (b) flutter frequency of the sample model with actuator (reference Mach number = 1.4).
1
-30
-20
-10
0
10
20
30
Sweep angle at mid chard,0 (degree)
Fig. 8. Effects of sweep angle on the (a) Butter speed and (b) flutter frequency of the sample model with torsional spring (reference Mach number = 1.4).
Flutter
analysis
(a) 5
compatibility
Stiffness
line
Stiffness
Ratio,r(k/ko)
Fig. 9. Effects of the torsional spring magnitude on the (a) flutter speed and (b) flutter frequency for various sweep angles (reference Mach number = 1.4).
sweep angle vs flutter
speed
graph
changes
substan-
The flutter frequency increases as the stiffness increases. Figure 9 shows how the flutter speed changes as the stiffness changes at a particular sweep angle: 0”; 10” and 20”. There exists a particular stiffness magnitude where flutter speed becomes a minimum at each configuration. the
599
the dynamic properties of the control actuator for the elastic launch vehicle wing. It is shown that the fictitious mass method can be applied effectively in obtaining accurate results efficiently in the flutter analysis for the aeroelastic system with an actuator. It is observed that there exists a sweep back angle at which the flutter speed becomes the maximum value for the actuator used in this study. For the case with static stiffness, there exists a sweep angle where the flutter speed becomes a minimum. It must be noted that the cases with static stiffness show a significantly different trend from those with an actuator and that an actuator cannot be replaced by a static spring with an appropriate stiffness. The dynamic characteristics of actuators must be properly considered for an accurate flutter prediction. REFERENCES
5oT-----l
as
vehicle
Ratio.r(k/kg)
(b)
tially
of launch
magnitude
of
stiffness
changes.
CONCLUSIONS
In this paper, the root locus method and iterative V-g method are used for flutter analysis considering
1. A. Barai and S. Durvasula. Flutter of shaft-supported low aspect-ratio control surfaces. J. Aeronaut. Sot. India 42, 19-26 (1990). 2. A. B. Kondrat’ev and V. A. Chashchin, On the stiffness of a flight vehicle pneumatic control actuator. IZO. Vi/Z. Auiarsion. Tekh. 33, 7--10 (1990). of 3. M. Karpel and C. D. Wieseman, Time simulation flutter with large stiffness changes. In: Proc. 33th AIAA/ASME/ASCEIAHS/ACS Structural Dynamics and Materials Corzference, paper 92-2394, Dallas, TX, pp. 1878-1886 (1992). 4. M. Karpel, Efficient vibration mode analysis of aircraft with multiple external store configuration. J. Aircr. 25, 747-751 (1988). 5. E. Livne, Accurate calculation of control-augmented structural eigenvalue sensitivities using reduced-order models. AIAA .I. 27, 947-954 (1989). math modeling methods for 6. K. L. Roger, Airplane active control design, AGARD-CP-228. pp. 4-14-l I (1977) 7. R. R. Craig Jr and M. C. C. Bampton, Coupling of substructures for dynamic analysis. AIAA. J. 6, 13131319 (1968). 8. J. T. Weissenburger, Effect of local modification on the vibration characteristics of linear systems. Trans. ASME Ser. E 35, 327.-332 (1968). 9. 1. Lee and J. J. Lee, Vibration analysis of composite plate wings. Compur. S/ruci. 37, 1077.-1085 (1990). 10. T. Ueda and E. H. Dowell, Doublet-point method for supersonic unsteady lifting surfaces. AIAA J. 22, 179.-186 (1984). II. R. L. Harder and R. N. Desmarais, Interpolation using surface splines. J. Aircr. 9, 189-191 (1972). 12. J. A. Nelder and R. Mead, A simplex method for function minimization. Compur. J. 7, 279-285 (1965).