- Email: [email protected]
Numerical solution of the rotating beam equations I. Flapping equations
IMPACT
OF COMPUTING
IN SCIENCE AND
ENGINEERING
2, 157- 185 ( 1990)
Numerical Solution of the Rotating Beam Equations I. Flapping Equations*
T. F. CHEN, G. J. FIX, AND R. KANNAN Department
of Mathematics,
University Arlington,
of Texas at Arlington, Texas
P.O. Box
19408,
76019
Received August 7, 1989
T. F. Chen, G. .I. Fix, and R. Kannan, Numerical Solution of the Rotating Beam Equations. I. Flapping Equations, IMPACT of Computing in Science and Engineering 2, 157-I 85 ( 1990). In this paper we consider the numerical simulation of the rotating beam equations. The latter are widely used as a model of helicopter rotor blade dynamics. Both spectral and finite element methods are studied, and numerical results are reported with both linear and nonlinear aerodynamic loading. 0 1990 Academic Press, Inc.
1. INTRODUCTION The rotating beam equations provide an effective model for helicopter rotor blade dynamics [l-3]. In this context a number of interesting phenomena which present significant challenges to both the theoretical and the numerical analysis of this problem can arise. First, it is the existence and stability of time periodic solutions that are of greatest interest. The analysis of such cases (theoretical as well as numerical) is inherently more complicated than the analysis of equilibrium solutions or time transient solutions. Second, the spatial operators are fourth order, whereas most of the previous analysis has centered around second-order operators (see, e.g., [4-61). At first glance, it might be hoped that the results from the latter could be carried over to the fourth-order case with some minor technical revisions. Both operators are elliptic, and the only thing that is lost in passing from the second* The work was supported in part by a grant from Bell Helicopter Textron, Inc., under Contract 344499-33. The authors acknowledge the significant contribution to this work made by Dr. J. G. Yen and Dr. Chen Chao of Bell Helicopter Textron, Inc. 157 0899-8248/90
$3.00
Copyright 0 1990 by Academic Press, Inc. All rights of reproduction in any form reserved.
158
CHEN,
FIX.
4ND
KANNAN
to the fourth-order case is the existence of maximum principles. From this research the authors have learned how significant this loss actually is. Moreover, this point is manifested in many unexpected ways. This was also recently indicated in a slightly different context in [ 71. The third aspect of this problem that is worth mentioning is that the complete equations constitute a system of fourth-order time-dependent partial differential equations. Because of rotational effects, i.e., existence of Coriolis forces, there is intriguing coupling between these equations which can lead to physically interesting instabilities [S, 91 in the time periodic solutions. Finally, composite materials are increasingly finding application in the modem helicopter blade. Thus it is of fundamental importance that numerical approximations be able to handle variable material properties and resolve effects due to anisotropy on blade dynamics. In this paper, the first of a planned series, we deal only with the blade flapping equation. A variationally based finite element scheme is introduced, and its effectiveness is studied for a class of model problems. The formulation is amenable to extensive generalizations including variable material properties. nonstandard geometric configurations such as “blade boxes” (see Section 3 ), and the full coupled equations for flapping, lead-lag, and torsional degrees of freedom. These will be considered in a sequel to this paper. In this paper, we shall also consider spectral approximations based on Chebyshev expansions of the solution using the same variational formulation that is used in the finite element scheme. Spectral methods generally find much less application in the analysis of rotor blade dynamics since they in effect require constant material properties in order to be effective. Nevertheless, when they are applicable-the case of a “homogenized blade” [I] is but one example-they can yield extraordinarily accurate approximations. We shall illustrate this in our selected examples. There is also the interesting possibility of using spectral domain decomposition techniques for the inhomogeneous case. See, for example, [ 15 1.
2. THE FLAPPING
EQUATIONS
Consider the rotating beam shown in Fig. 2.1. In this figure Y denotes the nondimensional radial distance (0 G r < 1)) $ the aximuthal angle (0 G $ < 27r), and z the nondimensional out-of-plane displacement. The beam is rotating with angular velocity 0. The associated dynamical equation in nondimensional form [I] is
NUMERICAL
SOLUTION
OF
THE
ROTATING
BEAM
EQUATIONS,
I
159
r=l
r=O FIG. 2.1. Rotating
beam.
where EZ is a sectional stiffness property, M the mass, and R the span of the beam. The source term f represents the aerodynamic loading, and it is a function of $, r, &far, azpq. The boundary conditions relevant to rotor blade dynamics are the following. At r = 0 the blade is pinned and moment free:
z(+, 0) = 0,
2 (I), 0) = 0.
(2.2)
The tip r = 1 is a free end: (2.3)
The transient analysis is determined by specifying initial conditions:
40, r) =
z0(r),
$ (0, r) = z,(r).
However, the case of greatest interest, and the one that will be the focus of this paper, is the time periodic case: z(# + 27r, r) = z($, r).
(2.5)
The normal beam equations admit a variational formulation in terms of strain energy minimization. Here the inertial term a2z/w2 and the aerodynamic term f are regarded as loads, and the strain energy
160
C-HEN.
FlX.
AND
where K = a’:/@’ --J; is minimized condition in this formulation is
KANNAN
at each $. The only essential boundar?
The other conditions in (2.2)-( 2.3) are natural and need not be applied to the trial functions z( * , * ). The latter should also have finite strain energy in the sense that is,’
[&($)‘+(q(q}d.<
a.
(2.8)
See [lo]. In this as well as subsequent papers we shall deal with the equivalent principle of virtual work [lo]. The latter characterizes the solution as a function z(*, .) such that for each fixed rc/both (2.7) and (2.8) are satisfied. In addition, it is required that for fixed $ the relation
hold for all admissible test functions 4 = t(r), i.e., functions 4 which vanish at r = 0 and have finite strain energy. The main advantage of the principle of virtual work is that it is easily generalized to the case of systems, where more than one degree of freedom is involved (e.g., flapping, lead-lag, and torsional modes). The spectral and finite element discretizations derived in Section 3 are based on (2.9). We conclude this section with some comments about the aerodynamic loading termf. Two cases are considered in this paper, one being linear and the other being nonlinear. In industrial design the aerodynamic loading is modeled by [ 111 F = i $
(r + p sin +)%‘Cr,
where p is the density, I* is the advanced ratio, C is the chord length, and C, is the lift force coefficient. The latter is a function of the blade element angle of attack (Y,which in turn depends on r, 1+5, dz/dr, and &/a$. For small angle assumption, the blade element angle of attack can be expressed as a=O(r,$)+
A-p$cos$-* (
all/
)i
(r + p sin $1.
(2.11)
NUMERICAL SOLUTION OF THE ROTATING BEAM EQUATIONS, I
161
The second term in (2.11) captures Coriolis effects, and X is the inflow ratio which represents the flow through the rotor disc plane. The first term represents pilot control input and is usually modeled via a Fourier expansion in 1c/.It is assumed to be independent of z. The relationship between lift C, and angle of attack (Y is empirical and consists of a linear regime (small a) as well as a nonlinear regime. The specific relationship used in this paper is given in Section 6.
3. FINITEELEMENTANDSPECTRALANALYSIS
The principle of virtual work can be used in conjunction with appropriate finite dimensional function spaces to discretize the equations of motion. In particular, this combination will reduce the partial differential equations to a set of ordinary differential equations. Two choices are considered in this paper. First, a spectral method based on expansions in terms of Chebyshev polynomials is considered. Then an alternate finite element formulation using piecewise bicubic Hermite polynomials is introduced. The only essential boundary conditions required in the variational principle (2.9) are that test and trial functions vanish at r = 0. Thus the modified Chebyshev polynomial Tk(r) = rcos{ (k - 1)arc cos(2r - 1))
(3.1)
of degree k ( 1 G k < n) can be used in this formulation. In particular, for a fixed integer n we expand our approximate solution-denoted z,,( * , . )-by
zdr,
$)
=
i j=l
aj(lc/)Tj(r).
(3.2)
Substituting this into (2.9) and requiring that this relation hold for
give the system of n ordinary differential equations M-a + Ka = f($, a, i)
(3.3)
in terms of the vector of weights (3.4)
162
C-HEN.
FIX.
AND
KANN4N
In (3.3), we have put da a=drC3
d’a
(3.5)
a=w’
where M, K are the mass and stiffness matrices (see [lo]). Finite elements offer an alternative version of ( 3.3). In particular, let 7r = { ro, rl, . . . , r,} be a grid on [0, 11: 0 = r. < r, < - - - -c rm = 1
Let ‘pi,0 and pj,l be the two C’ cubic Hermite basis functions associated with thejth node (see [lo]). These are displayed in Fig. 3.1. Then instead of (3.3) we expand the approximate solution-denoted by z,( . , . ) in this case-by
&Jr, $1 =5z,,o(+)P,,o(r) +zz,,l(#)cp,,l(r)> (3.6)
J=I
J=o
where the unknowns now are displacements
z,,o(rl/)= z,(r,, G),
(3.7)
and strains (3.8)
Again we obtain from (2.9) a systsem of IZ = 2m + 1 ordinary differential equations, where in this case the vector a( rc/)of unknowns is
a(+) = (zo,I($),
ZLO(+),
ZLI($),
. . , z,,o(#),
z,,~(lci),
. . . ,JT.
Two different time integration schemes are used in the numerical results presented in Sections 5 and 6. Both are based on the equivalent system
FIG.
3.1. Herrnite
functions.
NUMERICAL
SOLUTION
OF THE ROTATING
il=v 1Adi + Ku =f($,
BEAM EQUATIONS,
I
163 (3.9)
u, v),
which we write as
0 U
Ml++
=
gtik
WI,
W=
v
(3.10) .
In the linear case, the conservative second-order scheme
Ml{W (/+I)_ w(‘)> = Afgth,,z, w c/+1/29
(3.11)
is used, Here a time grid 0 = & < $, < $2 < + * * is introduced, and wu+1/2) = 1
*(W
(I) + w(l+l)
1,
with w’ being the approximation to w( I/‘). A distinguishing feature of this approximation is the absence of any numerical dissipation. Once the 2a periodic solution is approached in time, there is absolutely no attenuation of the waveform (see Section 5 ) . In the nonlinear case the approximation (3.11) is not very useful because the lift vs angle of attack data are given in tabular form. As an alternative, a fourth-order explicit Runga-Kutta scheme is used in Section 6 which requires only table lookups for g( II/, w) to march in time. This scheme goes as follows: Rewriting (3.10) as
6 = K’g($, WI= a+, w) (note that the entries in Mi depend on $), we let wI = w(j) + Arl/g($,, w(I))
(3.12)
w2 = %(rl/[ + $A$, w(‘) + fA$w,)
(3.13)
= ij(#, + iA+, w(‘) + fA$w2)
(3.14)
~3
w4 = P(IC/,+ A+, w(I) + A+w,).
(3.15)
These are the predictors. The value for w ($) at the next step $r+i is
w([+l) = w(I) + y
(w, + 2w2 + 2w3 + w4).
(3.16)
164
CHEN.
FIX.
Unlike (3.1 1), the Runga-Kutta and in particular
ALD
KAhNAh
scheme does have time-step restrictions.
must be of order ( 7 I ‘, where 171 = min (Y,+i - r,I
(see [ 131). Within this restriction, on the other hand, there is no numerical dissipation. This is seen from the waveforms presented in Section 6. Both the spectral and the finite element formulations have applications. For isotropic materials with uniform section properties and smooth displacement fields z(r, $), the former has exceptional accuracy. In particular, the solution (3.2) to (3.3) converges as n --* cc exponentially fast [lo]. It is true that the time-discretized systems (3.11) and (3.12)-( 3.16) have only secondand fourth-order accuracy in rc/.However, the 2~ periodic solution is rapidly formed, and thus $-grid effects are not important (see Sections 5 and 6). The finite element schemes, on the other hand, are less accurate, and in particular ( 3.6) converges like ] r I 4 (see [lo]). However, such schemes are very useful for composite materials, where the elastic properties EZ vary with Y. In addition, quasi two-dimensional structures such as the blade box shown in Fig. 3.2 are amenable to the finite element formulation. Such configurations will be investigated in a sequel to this paper. It is noted that only one of the four boundary conditions need be imposed on the approximate scheme. It is to be emphasized that in approximations based on virtual work, it is generally better not to force the natural boundary conditions. See [ 141. In addition, the magnitudes of
(3.17)
will give a posteriori indicators on the errors present in the second and third derivatives.
1
I
4
FIG. 3.2. Blade box.
NUMERICAL SOLUTION OF THE ROTATING BEAM EQUATIONS, I
165
4. FLOQUET ANALYSIS As mentioned in the Introduction, it is the 27r periodic solutions of (2.9)that are of the greatest practical interest. In this regard, one would want to know in what circumstances a periodic solution does exist, and of equal importance, one would want to know when it is stable. A closely related issue is the determination of those initial conditions (2.4) which will evolve into the periodic solution as 1c/increases. The numerical results presented in Sections 5 and 6 indicated, at least for small angles of attack CY,that there is a stable periodic solution which is globally attractive; i.e., it can be reached with arbitrary initial data. In this section, we confirm these properties under simplifying assumptions. A more elaborate analysis under realistic physical assumptions will be given in a sequel to this paper. In particular, we shall concentrate on the finite dimensional equations (3.3) for the case of linear aerodynamics. We rewrite (3.3) as Mii + Czi + Ka = g($).
(4.1)
lJ in (3.2) for the Let cp1(r), *. . 9 cp,( r) denote the basis (i.e., the functions spectral method or the function (01,~in (3.6) for the finite element method). Then the (j, k) entries of the matrices M, C, K are I Mjk =
s 0 Vj(~)(ok(YW
(4.2)
Cjk = C(ko’+ sin $C$’
(4.3)
Kjk = KjE) +
cos
$[ K$’ + sin $K$‘],
(4.4)
where 1 (0) cjk
_ -
c
s0
K;;)=
((T positive constant)
rvj(r>(pk(r>dr
s,‘[[&]
yv+
(4.5)
(4.6)
(!$)!!$d.+..
The entries of the other matrices C$‘, M$’ (1 = 1, 2 ) are analogous inner products Of pj and pk. Actually these matrices are not important to our analysis. It is only the mean values that are important. Note that &
s
fr C(T)&
= C(O),
+
sz* KIT ?I- 0
= K(O).
(4.7)
166
CHEN.
FIX,
AND
KANNAN
We start be rewriting (4.1) as a first-order system as in ( 3.9). I
0
where u = a and v = da/dt.
0
v +
g ’
I[1[I
-A!-‘C($)
-hr’K($)
U
(4.8)
The solution matrix for (4.8) is
SC+) = em
0
I
-W’Iz(7)
-M-q+)
I +,
(4.9)
where
Indeed, given the initial conditions these initial conditions is
Floquet theory (see [ 121) states that (4.8) has a unique 2~ periodic solution if the matrix I - S(2P)
(4.12)
is invertible. If in addition the eigenvalues of S( 2~) are strictly less than 1 in magnitude, then the periodic solution is asymptotically stable. The latter will hold when the eigenvalues (i.e., the so-called Floquet exponents) of
1
-~-lfp'
0
I -M-'c(O)
1
(4.13)
have a negative real part for it is exactly in that case where the exponential S( 27r) of this matrix has eigenvalues that are less than 1 in magnitude. Note that the eigenvalue problem for (4.3 ) can be written
1
0 -M-l@0
+&(o,][:]
= $:I,
(4.14)
NUMERICAL SOLUTION OF THE ROTATING BEAM EQUATIONS,1
167
and that this is equivalent to
[ x2A4-t xc’o’ + K’O’]u = 0.
(4.15)
Suppose X is an eigenvalue. Then taking the dot product with u gives x2u - (Mu) + xu * ( c’“‘u)
+ u - K’O’U = 0.
(4.16)
Thus
x = -f(c/m) rt d$(~/m)~- k,
(4.17)
c = u. C’O’U2
(4.18)
where m = u*(Mu),
k= u.K’O)u
We recall that the mass matrix M is positive definite. Also, the same is true of the mean value C’(O)of the damping matrix as well as the mean value K(O) of the stiffness matrix. We conclude that each of the numbers in (4.18) is positive and hence ReX < 0 follows from (4.15). This proves the following: THEOREM. Both thefinite element method and the spectral method produce a set of ordinary d$erential equations which has a unique 2~ periodic solution. In addition, the solution associated with any arbitrary set of initial conditions will approach this periodic solution.
5. NUMERICALRESULTS-LINEARCASE
The results presented in this section are for the case of a small angle of attack (Y,where the lift CI is a linear function of the latter. The general functional relationship between these two variables is given in the next section (see Fig. 6.1). Various choices for the number N of trial functions were used for the spectral/Chebyshev method. Ten percent accuracy was achieved with as few as N = 5 functions. The results presented in this section used N = 20 functions, although these results are essentially indistinguishable from the choice for N = 10 functions. This reflects the exponentiallyfast convergence of the spectral approximation.
168
CHEN.
FIX.
AND
KAXNAN
The choice of a G-step represents a more sensitive issue. The second-order implicit scheme described in Section 3 is unconditionally stable. However, il the +-step is too large, than the accurate spatial resolution is compromised. This can typically be seen in the values of
which should be close to zero in spite of the fact that the second of (2.2) and (2.3) have not been imposed on the trial functions. For the case of N = 20 Chebyshev functions, a $-step of r/36 was found to be near optimal. The $ behavior of the approximation is shown in Figs. 5.1 and 5.2 for I = 1 and r = 0.475, respectively. Observe that the blade tip moment (@z/ I%*)( 1, $) is zero to the scale shown for all rc/.The double peaks in the moments at I = f are physical, and the spectral method resolves these perfectly. The results from finite element approximation using piecewise bicubic Hermite functions are shown in Figs. 5.3 and 5.4. The choices of 20 radial intervals, or 4 1 Hermite functions, gave a system which was computationally equivalent to the spectral system for N = 20. The point here is that the finite element method produces a sparse matrix. while the matrix associated with the spectral method is full. Thus for a given computational effort, more functions can be used in the former method. The approximations produced by the finite element approximation are close to those produced by the spectra1 method. The latter is slightly more accurate as can be seen from the noticeable yet slight oscillations in the blade tip moments in Fig. 5.3 for the finite element approximation. Nevertheless. the double peak moment midblade effect is adequately resolved (see Fig. 5.4) by the finite element approximation. It should be noted that in both approximations the 27~periodic solution is quickly established. The results shown in Figs. 5.1 to 5.4 were evolved from a quiescent blade, i.e., z( r, 0) = $
(r, 0) = 0.
(5.2)
Other initial conditions were also used and the periodic state established itself in a similar manner as predicted from the discussion in Section 4. Another point to be noted from Figs. 5.1 to 5.4 is that there is absolutely no numerical dissipation in the waveforms, once they are established. Because of this a Fourier analysis in J/ can be used to represent displacements z, slopes dzldr, and moments a2z/dr2 as w(r,
$) = wO(r) + w,(r)sin
# + w2(r)cos
1+5+ + - -.
(5.3)
NUMERICAL
L? 0 d ‘0.00
SOLUTION
4.00
OF
THE
I 12.00
8’. 00
d(i
00
n
6’. 00
4.00
16.00
BEAM
EQUATIONS,
20.00
24.00
I
169
I 28.00
radians)
16.00
12.00
*(in
3. 0 iz
ROTATING
20.00
1 24.00
1 28.00
radians)
c .-._ 0
2 d‘0.00
4,oo
6.00 b(in
FIG.
5.1. Time
(angle
I 16.00
12.00
2o.oc
24.00
radians)
I/J) history
of blade tip-spectral
method.
1
28.00
CHEN. FIX. AND KANNAN
d
z d ‘0.00
I
I
5.00
4.00
12.00
O(in
z d‘0.00
4.00
6.00
16.00
24.00
,
7 26.00
24.00
L 2b. oc
radians)
12.00
b(in
1
2o.oc
16.00
2O.OC
radtans)
FIG. 5.2. Time (angle $) midblade--spectral
method.
NUMERICAL
SOLUTION
OF THE ROTATING
BEAM EQUATIONS,
I
171
L L N
6 d : s d‘o.oc
4.00
6.00
12.00 Jl(i
‘0’. oc
n
16.00 red
20.00 r
24.00 1
28.00 I
20.00
24.00
28.00
i arts.1
1
4.00
8.00
12.00 6(i
n
16.00 radians)
0, d
3
.
0
z
.
0
s:
N
:: _ d
$(in
radians)
FIG. 5.3. Time (angle 4) history at blade tip-finite
element method.
172
00
4.00
6.00
4.co
6.00 1
12.co
?J.OC
24.00
7 28.00
i6.00
20.0s
21. I oc
26. 1 cc
16.00
2O.CC
21.00
26. cc
16
12.co 1
00
radIansI
$(I"
d
I
1o.oo
4.00
E.09
12.00 JI( f n
FIG.
5.4. Time
(angle
$) history
rad1alTs)
at midblade-finite
element
method.
NUMERICAL
SOLUTION
OF THE
ROTATING
BEAM
EQIJATIONS,
I
(0 0 6
s d L L N s d
Fs d‘o.oc
0.23
0.40
0.60 r
0.80
1.00
1 1.20
0.23
0.40
0.60
0.80
1 1 .oo
I 1.20
;: d L N 8 d
00
FIG.
5.5. Steady
Fourier
component-spectral
method.
173
174
CHEN,
FIG. 5.6. sin Fourier
FIX.
AND
KANNAN
component-spectral
method.
NUMERICAL
SOLUTION
OF THE ROTATING
BEAM EQUATIONS,
;s d
1
0’ d L L N SY d
8 Oo.00
1
0.20
0.40
0.60
I
0.60
I
.oo
1 1.20
I ,.oo
1.20 f
r
2 0 L N s d
8 “0.00
0.20
0.40
0.60
0.80
r
FIG. 5.7.
cos Fourier component-spectral
method.
I
175
CHEN.
FIX.
AND
KANNA>
r\ .oo
0.23
0.40
0.50
0.60
I .oo
I
1.20
zo,, 0.230.50 0.60 0.40
!L-----v
c. 23 FIG.
5.8. Steady
c.40
Fourier
0.50
component-finite
0.60
element
i .!I0
1.23
i .!I0
i .2c
method
NUMERICAL
SOLUTION
d ‘0.00
I
0.23
N l-Iizhd20
OF
0.40
THE
ROTATING
I
0.60
BEAM
0.60
EQUATIONS,
I .oo
I
8 1.20
0.23 r
0 d ‘0.00
0.20
FIG.
0.40
5.9. sin Fourier
0.60 r
component-finite
0.60
element
1 .oo
method.
1
I .20
177
178
CHEN,
FIX.
AND
KANNAN
E d
0" 0 L
zoz
C.23
0.40
0.50
0.60
I .oo
i.2c
c.23
0.40
0.50
0.60
I .oo
1.23
2 d
:: “0.00 FIG.
5.10. cos Fourier
component-finite
element
method.
I
NUMERICAL
SOLUTION
OF
THE
ROTATING
BEAM
EQUATIONS,
I
179
The harmonics Wj(I) are generally of greatest design interest for each of these three variables. The absence of numerical dissipation means that these modes can be captured over any 27r J/-interval, once the periodic solution is established. Results for the first three harmonics are given in Figs. 5.5 to 5.7 for the spectral method, and Figs. 5.8 to 5.10 for the finite element method. The differences in the harmonic representations of these two methods are virtually indistinguishable.
6.
NUMERICAL
RESULTS-NONLINEAR
CASE
The results in this section are for the case of nonlinear aerodynamics. The functional relationship between the lift Cl and angle of attack (Yis shown in Fig. 6.1. This relationship is empirical and is typical of that used in industrial design. The linear regime falls in the range -0.2 < (Y G 0.2. Observe that outside this range there are a number of points where the data become somewhat rough. It is exactly for this reason that the fourth-order, explicit RungaKutta method for 1c,integration is useful, since it requires only numerical values of C, for various CX,and not derivatives with respect to (Y. Since future work will concentrate on cases with variable material properties where finite elements are preferred, only results using this method are given. Ten intervals and a $-step of X/ 120 radians were used. For the explicit scheme a $-step of say 71-/ 36 or even x/72 was found to introduce unacceptable errors. To check the reliability of the approximation the system parameters were chosen so that it should be approximately linear, yet the numerical scheme did not use this. The $ history is shown in Figs. 6.2 to 6.3. Again the 2n
:: 7
-3.23
8
-2.20
-I .20
-0.20 a(in
FIG.
6.1 Empirical
0.60
I
1.60
radians)
data:
lift vs angle of attack.
2.60
3.60
CHEN.
FIX.
4ND
lL4NNAN
0 d‘o.oc
4.00
6.00
12.00
6(i
n
16.00
20.05
24.00
26.CC 3
16.00
20.00
24.00
26.OC I
radians)
0 -
-. 0 ‘0.00
4.00
6.00
12.00
tiCin
FIG. 6.2. Nonlinear
radians)
aerodynamics-blade
tip time history
NUMERICAL
00
4.00
SOLUTION
OF THE ROTATING
1 12.00
6.00
@(in
radians)
$(in
radians)
BEAM EQUATIONS,
4.00
6.00
12.00
@(in FIG.
181
16.00
2o.oc
24.00
1 28.00
16.00
20.OC
24.00
28.00
“I d
‘o.oc
I
radians)
6.3. Nonlinear aerodynamics-midblade
time history.
I
182
CHEN.
FIX.
AND
KANNAN
Nb-:r_-:::I c.23
c.40
0.50 r
0.60
i .oo
N
FIG.
6.4. Nonlinear
aerodynamics-steady
Fourier
component.
I .20
NUMERICAL
SOLUTION
OF THE ROTATING
BEAM EQUATIONS,
I
;: dL N g d-
N d-‘0.00
0.23
0.40
0.60 1
0.60
1 I.oo
1.20 1
0.80
1 .oo
I 1.20
r
d ‘0.00
c.20
0.40
0.60 I-
FIG. 6.5. Nonlinear aerodynamics-sin
Fourier component.
183
CHEN.
cc
:I
C.20
FIX.
0.40
c.23
0.40
AND
FIG. 6.6. Nonlinear
c.4c
I
0.50
0.60
1 .oo
1.23
0.60
0.60
1.3G
i .20
0.60
i .CG
, .20
N;c; c.23
KANNAh
0.50 r
aerodynamics-cos
Fourier
component.
NUMERICAL
SOLUTION
OF
THE
ROTATING
BEAM
EQUATIONS,
I
185
periodic was quickly established, and there was no attenuation of the waveforms due to artificial numerical dissipation. Moreover, the midblade double peak moment effect was satisfactorily resolved. There was, however, slightly larger oscillations in the tip blade moments (compare Fig. 6.2 and Fig. 5.3). This is in part due to the rough data and is very sensitive to the $-step. For example, it was found that the oscillations increased fourfold when the $-step was increased from a/ 120 to a/72. It should also be noted that only ten spatial intervals were used. This is why the results are slightly less accurate than those reported in Section 4. The harmonic representation (5.3) was more robust. The results in Figs. 6.4 to 6.6 are close to the linear case, and it was found that they were not greatly affected by the $-step provided the latter was r/72 or lower.
REFERENCES communication
(Bell
2. R. L. Bisplinghoff
1. F. Harris,
and H. Ashley,
Principles ofderoelasticity.
3. P. P. Friedmann, 1041 (1977).
Recent
4. L. Cesari
Private
developments
and R. Kannan, Solutions ( 1982).
Textron). in rotary
of nonlinear
wing
Wiley,
New York
hyperbolic
equations
( 1962).
.I. Aircraf 14, 1027-
aeroelasticity.
at resonance.
Nonlinear
Anal. 6,751-805 5. A. Haraux,
Nonlinear
evolution
equations-Global New
in Mathematics, Vol. 841. Springer-Verlag,
6. J. K. Hale, Xiao-Biao Lin, and Genevieve approximations of semigroups and partial
behaviour York.
Raugel, differential
Upper semicontinuity equations. Math.
7. A. Friedman and L. Oswald, The blow-up time for higher-order with small leading coefficients. J. Differential Equations 75, 8. K. B. Subrahmanyam, and stability of rotating, 342-352 (1987).
of solutions.
semilinear
In Lecture Notes of attractors parabolic
K. R. V. Kaza, G. V. Brown, and C. Lawrence, pretwisted, preconed blades including Coriolis
Nonlinear vibration effects. J. Aircraft 24,
10. G. Fix and G. Strang, Cliffs, NJ (1973). 11. T. F. Chen, Department
12. L. Cesari,
G. J. Fix, and R. Kannan, of Mathematics, University
Rotor Blade Dynamics Project, Reports of Texas,
Arlington,
equations
239-263 ( 1988).
Helicopter Theory. Princeton Univ. Press Princeton, NJ, ( 1980). An Analysis of the Finite Element Method. Prentice-Hall,
9. W. Johnson,
for
Camp. ( 1987).
and Bell Textron
Englewood 1 and 2. ( 1988).
Asymptotic Behavior and Stability Problems in Ordinary D@erential Equations.
Springer-Verlag,
New York.
13. G. Dahlquist and A. Bjorck, Numerical Analysis. Prentice-Hall, Englewood Cliffs, NJ ( 1978). 14. J. Starch and G. Strang, Paradox lost: Natural boundary conditions in the Ritz-Galerkin method. Internut. J. Numer. Methods. Engrg. 26, 2255-2266 (1988). 15. C. Bernardi
and Y. Maday,
Comput. & Structures 30( l/2),
Spectral methods for the approximation 205-216 ( 1988).
of fourth
problems.
OF COMPUTING
IN SCIENCE AND
ENGINEERING
2, 157- 185 ( 1990)
Numerical Solution of the Rotating Beam Equations I. Flapping Equations*
T. F. CHEN, G. J. FIX, AND R. KANNAN Department
of Mathematics,
University Arlington,
of Texas at Arlington, Texas
P.O. Box
19408,
76019
Received August 7, 1989
T. F. Chen, G. .I. Fix, and R. Kannan, Numerical Solution of the Rotating Beam Equations. I. Flapping Equations, IMPACT of Computing in Science and Engineering 2, 157-I 85 ( 1990). In this paper we consider the numerical simulation of the rotating beam equations. The latter are widely used as a model of helicopter rotor blade dynamics. Both spectral and finite element methods are studied, and numerical results are reported with both linear and nonlinear aerodynamic loading. 0 1990 Academic Press, Inc.
1. INTRODUCTION The rotating beam equations provide an effective model for helicopter rotor blade dynamics [l-3]. In this context a number of interesting phenomena which present significant challenges to both the theoretical and the numerical analysis of this problem can arise. First, it is the existence and stability of time periodic solutions that are of greatest interest. The analysis of such cases (theoretical as well as numerical) is inherently more complicated than the analysis of equilibrium solutions or time transient solutions. Second, the spatial operators are fourth order, whereas most of the previous analysis has centered around second-order operators (see, e.g., [4-61). At first glance, it might be hoped that the results from the latter could be carried over to the fourth-order case with some minor technical revisions. Both operators are elliptic, and the only thing that is lost in passing from the second* The work was supported in part by a grant from Bell Helicopter Textron, Inc., under Contract 344499-33. The authors acknowledge the significant contribution to this work made by Dr. J. G. Yen and Dr. Chen Chao of Bell Helicopter Textron, Inc. 157 0899-8248/90
$3.00
Copyright 0 1990 by Academic Press, Inc. All rights of reproduction in any form reserved.
158
CHEN,
FIX.
4ND
KANNAN
to the fourth-order case is the existence of maximum principles. From this research the authors have learned how significant this loss actually is. Moreover, this point is manifested in many unexpected ways. This was also recently indicated in a slightly different context in [ 71. The third aspect of this problem that is worth mentioning is that the complete equations constitute a system of fourth-order time-dependent partial differential equations. Because of rotational effects, i.e., existence of Coriolis forces, there is intriguing coupling between these equations which can lead to physically interesting instabilities [S, 91 in the time periodic solutions. Finally, composite materials are increasingly finding application in the modem helicopter blade. Thus it is of fundamental importance that numerical approximations be able to handle variable material properties and resolve effects due to anisotropy on blade dynamics. In this paper, the first of a planned series, we deal only with the blade flapping equation. A variationally based finite element scheme is introduced, and its effectiveness is studied for a class of model problems. The formulation is amenable to extensive generalizations including variable material properties. nonstandard geometric configurations such as “blade boxes” (see Section 3 ), and the full coupled equations for flapping, lead-lag, and torsional degrees of freedom. These will be considered in a sequel to this paper. In this paper, we shall also consider spectral approximations based on Chebyshev expansions of the solution using the same variational formulation that is used in the finite element scheme. Spectral methods generally find much less application in the analysis of rotor blade dynamics since they in effect require constant material properties in order to be effective. Nevertheless, when they are applicable-the case of a “homogenized blade” [I] is but one example-they can yield extraordinarily accurate approximations. We shall illustrate this in our selected examples. There is also the interesting possibility of using spectral domain decomposition techniques for the inhomogeneous case. See, for example, [ 15 1.
2. THE FLAPPING
EQUATIONS
Consider the rotating beam shown in Fig. 2.1. In this figure Y denotes the nondimensional radial distance (0 G r < 1)) $ the aximuthal angle (0 G $ < 27r), and z the nondimensional out-of-plane displacement. The beam is rotating with angular velocity 0. The associated dynamical equation in nondimensional form [I] is
NUMERICAL
SOLUTION
OF
THE
ROTATING
BEAM
EQUATIONS,
I
159
r=l
r=O FIG. 2.1. Rotating
beam.
where EZ is a sectional stiffness property, M the mass, and R the span of the beam. The source term f represents the aerodynamic loading, and it is a function of $, r, &far, azpq. The boundary conditions relevant to rotor blade dynamics are the following. At r = 0 the blade is pinned and moment free:
z(+, 0) = 0,
2 (I), 0) = 0.
(2.2)
The tip r = 1 is a free end: (2.3)
The transient analysis is determined by specifying initial conditions:
40, r) =
z0(r),
$ (0, r) = z,(r).
However, the case of greatest interest, and the one that will be the focus of this paper, is the time periodic case: z(# + 27r, r) = z($, r).
(2.5)
The normal beam equations admit a variational formulation in terms of strain energy minimization. Here the inertial term a2z/w2 and the aerodynamic term f are regarded as loads, and the strain energy
160
C-HEN.
FlX.
AND
where K = a’:/@’ --J; is minimized condition in this formulation is
KANNAN
at each $. The only essential boundar?
The other conditions in (2.2)-( 2.3) are natural and need not be applied to the trial functions z( * , * ). The latter should also have finite strain energy in the sense that is,’
[&($)‘+(q(q}d.<
a.
(2.8)
See [lo]. In this as well as subsequent papers we shall deal with the equivalent principle of virtual work [lo]. The latter characterizes the solution as a function z(*, .) such that for each fixed rc/both (2.7) and (2.8) are satisfied. In addition, it is required that for fixed $ the relation
hold for all admissible test functions 4 = t(r), i.e., functions 4 which vanish at r = 0 and have finite strain energy. The main advantage of the principle of virtual work is that it is easily generalized to the case of systems, where more than one degree of freedom is involved (e.g., flapping, lead-lag, and torsional modes). The spectral and finite element discretizations derived in Section 3 are based on (2.9). We conclude this section with some comments about the aerodynamic loading termf. Two cases are considered in this paper, one being linear and the other being nonlinear. In industrial design the aerodynamic loading is modeled by [ 111 F = i $
(r + p sin +)%‘Cr,
where p is the density, I* is the advanced ratio, C is the chord length, and C, is the lift force coefficient. The latter is a function of the blade element angle of attack (Y,which in turn depends on r, 1+5, dz/dr, and &/a$. For small angle assumption, the blade element angle of attack can be expressed as a=O(r,$)+
A-p$cos$-* (
all/
)i
(r + p sin $1.
(2.11)
NUMERICAL SOLUTION OF THE ROTATING BEAM EQUATIONS, I
161
The second term in (2.11) captures Coriolis effects, and X is the inflow ratio which represents the flow through the rotor disc plane. The first term represents pilot control input and is usually modeled via a Fourier expansion in 1c/.It is assumed to be independent of z. The relationship between lift C, and angle of attack (Y is empirical and consists of a linear regime (small a) as well as a nonlinear regime. The specific relationship used in this paper is given in Section 6.
3. FINITEELEMENTANDSPECTRALANALYSIS
The principle of virtual work can be used in conjunction with appropriate finite dimensional function spaces to discretize the equations of motion. In particular, this combination will reduce the partial differential equations to a set of ordinary differential equations. Two choices are considered in this paper. First, a spectral method based on expansions in terms of Chebyshev polynomials is considered. Then an alternate finite element formulation using piecewise bicubic Hermite polynomials is introduced. The only essential boundary conditions required in the variational principle (2.9) are that test and trial functions vanish at r = 0. Thus the modified Chebyshev polynomial Tk(r) = rcos{ (k - 1)arc cos(2r - 1))
(3.1)
of degree k ( 1 G k < n) can be used in this formulation. In particular, for a fixed integer n we expand our approximate solution-denoted z,,( * , . )-by
zdr,
$)
=
i j=l
aj(lc/)Tj(r).
(3.2)
Substituting this into (2.9) and requiring that this relation hold for
give the system of n ordinary differential equations M-a + Ka = f($, a, i)
(3.3)
in terms of the vector of weights (3.4)
162
C-HEN.
FIX.
AND
KANN4N
In (3.3), we have put da a=drC3
d’a
(3.5)
a=w’
where M, K are the mass and stiffness matrices (see [lo]). Finite elements offer an alternative version of ( 3.3). In particular, let 7r = { ro, rl, . . . , r,} be a grid on [0, 11: 0 = r. < r, < - - - -c rm = 1
Let ‘pi,0 and pj,l be the two C’ cubic Hermite basis functions associated with thejth node (see [lo]). These are displayed in Fig. 3.1. Then instead of (3.3) we expand the approximate solution-denoted by z,( . , . ) in this case-by
&Jr, $1 =5z,,o(+)P,,o(r) +zz,,l(#)cp,,l(r)> (3.6)
J=I
J=o
where the unknowns now are displacements
z,,o(rl/)= z,(r,, G),
(3.7)
and strains (3.8)
Again we obtain from (2.9) a systsem of IZ = 2m + 1 ordinary differential equations, where in this case the vector a( rc/)of unknowns is
a(+) = (zo,I($),
ZLO(+),
ZLI($),
. . , z,,o(#),
z,,~(lci),
. . . ,JT.
Two different time integration schemes are used in the numerical results presented in Sections 5 and 6. Both are based on the equivalent system
FIG.
3.1. Herrnite
functions.
NUMERICAL
SOLUTION
OF THE ROTATING
il=v 1Adi + Ku =f($,
BEAM EQUATIONS,
I
163 (3.9)
u, v),
which we write as
0 U
Ml++
=
gtik
WI,
W=
v
(3.10) .
In the linear case, the conservative second-order scheme
Ml{W (/+I)_ w(‘)> = Afgth,,z, w c/+1/29
(3.11)
is used, Here a time grid 0 = & < $, < $2 < + * * is introduced, and wu+1/2) = 1
*(W
(I) + w(l+l)
1,
with w’ being the approximation to w( I/‘). A distinguishing feature of this approximation is the absence of any numerical dissipation. Once the 2a periodic solution is approached in time, there is absolutely no attenuation of the waveform (see Section 5 ) . In the nonlinear case the approximation (3.11) is not very useful because the lift vs angle of attack data are given in tabular form. As an alternative, a fourth-order explicit Runga-Kutta scheme is used in Section 6 which requires only table lookups for g( II/, w) to march in time. This scheme goes as follows: Rewriting (3.10) as
6 = K’g($, WI= a+, w) (note that the entries in Mi depend on $), we let wI = w(j) + Arl/g($,, w(I))
(3.12)
w2 = %(rl/[ + $A$, w(‘) + fA$w,)
(3.13)
= ij(#, + iA+, w(‘) + fA$w2)
(3.14)
~3
w4 = P(IC/,+ A+, w(I) + A+w,).
(3.15)
These are the predictors. The value for w ($) at the next step $r+i is
w([+l) = w(I) + y
(w, + 2w2 + 2w3 + w4).
(3.16)
164
CHEN.
FIX.
Unlike (3.1 1), the Runga-Kutta and in particular
ALD
KAhNAh
scheme does have time-step restrictions.
must be of order ( 7 I ‘, where 171 = min (Y,+i - r,I
(see [ 131). Within this restriction, on the other hand, there is no numerical dissipation. This is seen from the waveforms presented in Section 6. Both the spectral and the finite element formulations have applications. For isotropic materials with uniform section properties and smooth displacement fields z(r, $), the former has exceptional accuracy. In particular, the solution (3.2) to (3.3) converges as n --* cc exponentially fast [lo]. It is true that the time-discretized systems (3.11) and (3.12)-( 3.16) have only secondand fourth-order accuracy in rc/.However, the 2~ periodic solution is rapidly formed, and thus $-grid effects are not important (see Sections 5 and 6). The finite element schemes, on the other hand, are less accurate, and in particular ( 3.6) converges like ] r I 4 (see [lo]). However, such schemes are very useful for composite materials, where the elastic properties EZ vary with Y. In addition, quasi two-dimensional structures such as the blade box shown in Fig. 3.2 are amenable to the finite element formulation. Such configurations will be investigated in a sequel to this paper. It is noted that only one of the four boundary conditions need be imposed on the approximate scheme. It is to be emphasized that in approximations based on virtual work, it is generally better not to force the natural boundary conditions. See [ 141. In addition, the magnitudes of
(3.17)
will give a posteriori indicators on the errors present in the second and third derivatives.
1
I
4
FIG. 3.2. Blade box.
NUMERICAL SOLUTION OF THE ROTATING BEAM EQUATIONS, I
165
4. FLOQUET ANALYSIS As mentioned in the Introduction, it is the 27r periodic solutions of (2.9)that are of the greatest practical interest. In this regard, one would want to know in what circumstances a periodic solution does exist, and of equal importance, one would want to know when it is stable. A closely related issue is the determination of those initial conditions (2.4) which will evolve into the periodic solution as 1c/increases. The numerical results presented in Sections 5 and 6 indicated, at least for small angles of attack CY,that there is a stable periodic solution which is globally attractive; i.e., it can be reached with arbitrary initial data. In this section, we confirm these properties under simplifying assumptions. A more elaborate analysis under realistic physical assumptions will be given in a sequel to this paper. In particular, we shall concentrate on the finite dimensional equations (3.3) for the case of linear aerodynamics. We rewrite (3.3) as Mii + Czi + Ka = g($).
(4.1)
lJ in (3.2) for the Let cp1(r), *. . 9 cp,( r) denote the basis (i.e., the functions spectral method or the function (01,~in (3.6) for the finite element method). Then the (j, k) entries of the matrices M, C, K are I Mjk =
s 0 Vj(~)(ok(YW
(4.2)
Cjk = C(ko’+ sin $C$’
(4.3)
Kjk = KjE) +
cos
$[ K$’ + sin $K$‘],
(4.4)
where 1 (0) cjk
_ -
c
s0
K;;)=
((T positive constant)
rvj(r>(pk(r>dr
s,‘[[&]
yv+
(4.5)
(4.6)
(!$)!!$d.+..
The entries of the other matrices C$‘, M$’ (1 = 1, 2 ) are analogous inner products Of pj and pk. Actually these matrices are not important to our analysis. It is only the mean values that are important. Note that &
s
fr C(T)&
= C(O),
+
sz* KIT ?I- 0
= K(O).
(4.7)
166
CHEN.
FIX,
AND
KANNAN
We start be rewriting (4.1) as a first-order system as in ( 3.9). I
0
where u = a and v = da/dt.
0
v +
g ’
I[1[I
-A!-‘C($)
-hr’K($)
U
(4.8)
The solution matrix for (4.8) is
SC+) = em
0
I
-W’Iz(7)
-M-q+)
I +,
(4.9)
where
Indeed, given the initial conditions these initial conditions is
Floquet theory (see [ 121) states that (4.8) has a unique 2~ periodic solution if the matrix I - S(2P)
(4.12)
is invertible. If in addition the eigenvalues of S( 2~) are strictly less than 1 in magnitude, then the periodic solution is asymptotically stable. The latter will hold when the eigenvalues (i.e., the so-called Floquet exponents) of
1
-~-lfp'
0
I -M-'c(O)
1
(4.13)
have a negative real part for it is exactly in that case where the exponential S( 27r) of this matrix has eigenvalues that are less than 1 in magnitude. Note that the eigenvalue problem for (4.3 ) can be written
1
0 -M-l@0
+&(o,][:]
= $:I,
(4.14)
NUMERICAL SOLUTION OF THE ROTATING BEAM EQUATIONS,1
167
and that this is equivalent to
[ x2A4-t xc’o’ + K’O’]u = 0.
(4.15)
Suppose X is an eigenvalue. Then taking the dot product with u gives x2u - (Mu) + xu * ( c’“‘u)
+ u - K’O’U = 0.
(4.16)
Thus
x = -f(c/m) rt d$(~/m)~- k,
(4.17)
c = u. C’O’U2
(4.18)
where m = u*(Mu),
k= u.K’O)u
We recall that the mass matrix M is positive definite. Also, the same is true of the mean value C’(O)of the damping matrix as well as the mean value K(O) of the stiffness matrix. We conclude that each of the numbers in (4.18) is positive and hence ReX < 0 follows from (4.15). This proves the following: THEOREM. Both thefinite element method and the spectral method produce a set of ordinary d$erential equations which has a unique 2~ periodic solution. In addition, the solution associated with any arbitrary set of initial conditions will approach this periodic solution.
5. NUMERICALRESULTS-LINEARCASE
The results presented in this section are for the case of a small angle of attack (Y,where the lift CI is a linear function of the latter. The general functional relationship between these two variables is given in the next section (see Fig. 6.1). Various choices for the number N of trial functions were used for the spectral/Chebyshev method. Ten percent accuracy was achieved with as few as N = 5 functions. The results presented in this section used N = 20 functions, although these results are essentially indistinguishable from the choice for N = 10 functions. This reflects the exponentiallyfast convergence of the spectral approximation.
168
CHEN.
FIX.
AND
KAXNAN
The choice of a G-step represents a more sensitive issue. The second-order implicit scheme described in Section 3 is unconditionally stable. However, il the +-step is too large, than the accurate spatial resolution is compromised. This can typically be seen in the values of
which should be close to zero in spite of the fact that the second of (2.2) and (2.3) have not been imposed on the trial functions. For the case of N = 20 Chebyshev functions, a $-step of r/36 was found to be near optimal. The $ behavior of the approximation is shown in Figs. 5.1 and 5.2 for I = 1 and r = 0.475, respectively. Observe that the blade tip moment (@z/ I%*)( 1, $) is zero to the scale shown for all rc/.The double peaks in the moments at I = f are physical, and the spectral method resolves these perfectly. The results from finite element approximation using piecewise bicubic Hermite functions are shown in Figs. 5.3 and 5.4. The choices of 20 radial intervals, or 4 1 Hermite functions, gave a system which was computationally equivalent to the spectral system for N = 20. The point here is that the finite element method produces a sparse matrix. while the matrix associated with the spectral method is full. Thus for a given computational effort, more functions can be used in the former method. The approximations produced by the finite element approximation are close to those produced by the spectra1 method. The latter is slightly more accurate as can be seen from the noticeable yet slight oscillations in the blade tip moments in Fig. 5.3 for the finite element approximation. Nevertheless. the double peak moment midblade effect is adequately resolved (see Fig. 5.4) by the finite element approximation. It should be noted that in both approximations the 27~periodic solution is quickly established. The results shown in Figs. 5.1 to 5.4 were evolved from a quiescent blade, i.e., z( r, 0) = $
(r, 0) = 0.
(5.2)
Other initial conditions were also used and the periodic state established itself in a similar manner as predicted from the discussion in Section 4. Another point to be noted from Figs. 5.1 to 5.4 is that there is absolutely no numerical dissipation in the waveforms, once they are established. Because of this a Fourier analysis in J/ can be used to represent displacements z, slopes dzldr, and moments a2z/dr2 as w(r,
$) = wO(r) + w,(r)sin
# + w2(r)cos
1+5+ + - -.
(5.3)
NUMERICAL
L? 0 d ‘0.00
SOLUTION
4.00
OF
THE
I 12.00
8’. 00
d(i
00
n
6’. 00
4.00
16.00
BEAM
EQUATIONS,
20.00
24.00
I
169
I 28.00
radians)
16.00
12.00
*(in
3. 0 iz
ROTATING
20.00
1 24.00
1 28.00
radians)
c .-._ 0
2 d‘0.00
4,oo
6.00 b(in
FIG.
5.1. Time
(angle
I 16.00
12.00
2o.oc
24.00
radians)
I/J) history
of blade tip-spectral
method.
1
28.00
CHEN. FIX. AND KANNAN
d
z d ‘0.00
I
I
5.00
4.00
12.00
O(in
z d‘0.00
4.00
6.00
16.00
24.00
,
7 26.00
24.00
L 2b. oc
radians)
12.00
b(in
1
2o.oc
16.00
2O.OC
radtans)
FIG. 5.2. Time (angle $) midblade--spectral
method.
NUMERICAL
SOLUTION
OF THE ROTATING
BEAM EQUATIONS,
I
171
L L N
6 d : s d‘o.oc
4.00
6.00
12.00 Jl(i
‘0’. oc
n
16.00 red
20.00 r
24.00 1
28.00 I
20.00
24.00
28.00
i arts.1
1
4.00
8.00
12.00 6(i
n
16.00 radians)
0, d
3
.
0
z
.
0
s:
N
:: _ d
$(in
radians)
FIG. 5.3. Time (angle 4) history at blade tip-finite
element method.
172
00
4.00
6.00
4.co
6.00 1
12.co
?J.OC
24.00
7 28.00
i6.00
20.0s
21. I oc
26. 1 cc
16.00
2O.CC
21.00
26. cc
16
12.co 1
00
radIansI
$(I"
d
I
1o.oo
4.00
E.09
12.00 JI( f n
FIG.
5.4. Time
(angle
$) history
rad1alTs)
at midblade-finite
element
method.
NUMERICAL
SOLUTION
OF THE
ROTATING
BEAM
EQIJATIONS,
I
(0 0 6
s d L L N s d
Fs d‘o.oc
0.23
0.40
0.60 r
0.80
1.00
1 1.20
0.23
0.40
0.60
0.80
1 1 .oo
I 1.20
;: d L N 8 d
00
FIG.
5.5. Steady
Fourier
component-spectral
method.
173
174
CHEN,
FIG. 5.6. sin Fourier
FIX.
AND
KANNAN
component-spectral
method.
NUMERICAL
SOLUTION
OF THE ROTATING
BEAM EQUATIONS,
;s d
1
0’ d L L N SY d
8 Oo.00
1
0.20
0.40
0.60
I
0.60
I
.oo
1 1.20
I ,.oo
1.20 f
r
2 0 L N s d
8 “0.00
0.20
0.40
0.60
0.80
r
FIG. 5.7.
cos Fourier component-spectral
method.
I
175
CHEN.
FIX.
AND
KANNA>
r\ .oo
0.23
0.40
0.50
0.60
I .oo
I
1.20
zo,, 0.230.50 0.60 0.40
!L-----v
c. 23 FIG.
5.8. Steady
c.40
Fourier
0.50
component-finite
0.60
element
i .!I0
1.23
i .!I0
i .2c
method
NUMERICAL
SOLUTION
d ‘0.00
I
0.23
N l-Iizhd20
OF
0.40
THE
ROTATING
I
0.60
BEAM
0.60
EQUATIONS,
I .oo
I
8 1.20
0.23 r
0 d ‘0.00
0.20
FIG.
0.40
5.9. sin Fourier
0.60 r
component-finite
0.60
element
1 .oo
method.
1
I .20
177
178
CHEN,
FIX.
AND
KANNAN
E d
0" 0 L
zoz
C.23
0.40
0.50
0.60
I .oo
i.2c
c.23
0.40
0.50
0.60
I .oo
1.23
2 d
:: “0.00 FIG.
5.10. cos Fourier
component-finite
element
method.
I
NUMERICAL
SOLUTION
OF
THE
ROTATING
BEAM
EQUATIONS,
I
179
The harmonics Wj(I) are generally of greatest design interest for each of these three variables. The absence of numerical dissipation means that these modes can be captured over any 27r J/-interval, once the periodic solution is established. Results for the first three harmonics are given in Figs. 5.5 to 5.7 for the spectral method, and Figs. 5.8 to 5.10 for the finite element method. The differences in the harmonic representations of these two methods are virtually indistinguishable.
6.
NUMERICAL
RESULTS-NONLINEAR
CASE
The results in this section are for the case of nonlinear aerodynamics. The functional relationship between the lift Cl and angle of attack (Yis shown in Fig. 6.1. This relationship is empirical and is typical of that used in industrial design. The linear regime falls in the range -0.2 < (Y G 0.2. Observe that outside this range there are a number of points where the data become somewhat rough. It is exactly for this reason that the fourth-order, explicit RungaKutta method for 1c,integration is useful, since it requires only numerical values of C, for various CX,and not derivatives with respect to (Y. Since future work will concentrate on cases with variable material properties where finite elements are preferred, only results using this method are given. Ten intervals and a $-step of X/ 120 radians were used. For the explicit scheme a $-step of say 71-/ 36 or even x/72 was found to introduce unacceptable errors. To check the reliability of the approximation the system parameters were chosen so that it should be approximately linear, yet the numerical scheme did not use this. The $ history is shown in Figs. 6.2 to 6.3. Again the 2n
:: 7
-3.23
8
-2.20
-I .20
-0.20 a(in
FIG.
6.1 Empirical
0.60
I
1.60
radians)
data:
lift vs angle of attack.
2.60
3.60
CHEN.
FIX.
4ND
lL4NNAN
0 d‘o.oc
4.00
6.00
12.00
6(i
n
16.00
20.05
24.00
26.CC 3
16.00
20.00
24.00
26.OC I
radians)
0 -
-. 0 ‘0.00
4.00
6.00
12.00
tiCin
FIG. 6.2. Nonlinear
radians)
aerodynamics-blade
tip time history
NUMERICAL
00
4.00
SOLUTION
OF THE ROTATING
1 12.00
6.00
@(in
radians)
$(in
radians)
BEAM EQUATIONS,
4.00
6.00
12.00
@(in FIG.
181
16.00
2o.oc
24.00
1 28.00
16.00
20.OC
24.00
28.00
“I d
‘o.oc
I
radians)
6.3. Nonlinear aerodynamics-midblade
time history.
I
182
CHEN.
FIX.
AND
KANNAN
Nb-:r_-:::I c.23
c.40
0.50 r
0.60
i .oo
N
FIG.
6.4. Nonlinear
aerodynamics-steady
Fourier
component.
I .20
NUMERICAL
SOLUTION
OF THE ROTATING
BEAM EQUATIONS,
I
;: dL N g d-
N d-‘0.00
0.23
0.40
0.60 1
0.60
1 I.oo
1.20 1
0.80
1 .oo
I 1.20
r
d ‘0.00
c.20
0.40
0.60 I-
FIG. 6.5. Nonlinear aerodynamics-sin
Fourier component.
183
CHEN.
cc
:I
C.20
FIX.
0.40
c.23
0.40
AND
FIG. 6.6. Nonlinear
c.4c
I
0.50
0.60
1 .oo
1.23
0.60
0.60
1.3G
i .20
0.60
i .CG
, .20
N;c; c.23
KANNAh
0.50 r
aerodynamics-cos
Fourier
component.
NUMERICAL
SOLUTION
OF
THE
ROTATING
BEAM
EQUATIONS,
I
185
periodic was quickly established, and there was no attenuation of the waveforms due to artificial numerical dissipation. Moreover, the midblade double peak moment effect was satisfactorily resolved. There was, however, slightly larger oscillations in the tip blade moments (compare Fig. 6.2 and Fig. 5.3). This is in part due to the rough data and is very sensitive to the $-step. For example, it was found that the oscillations increased fourfold when the $-step was increased from a/ 120 to a/72. It should also be noted that only ten spatial intervals were used. This is why the results are slightly less accurate than those reported in Section 4. The harmonic representation (5.3) was more robust. The results in Figs. 6.4 to 6.6 are close to the linear case, and it was found that they were not greatly affected by the $-step provided the latter was r/72 or lower.
REFERENCES communication
(Bell
2. R. L. Bisplinghoff
1. F. Harris,
and H. Ashley,
Principles ofderoelasticity.
3. P. P. Friedmann, 1041 (1977).
Recent
4. L. Cesari
Private
developments
and R. Kannan, Solutions ( 1982).
Textron). in rotary
of nonlinear
wing
Wiley,
New York
hyperbolic
equations
( 1962).
.I. Aircraf 14, 1027-
aeroelasticity.
at resonance.
Nonlinear
Anal. 6,751-805 5. A. Haraux,
Nonlinear
evolution
equations-Global New
in Mathematics, Vol. 841. Springer-Verlag,
6. J. K. Hale, Xiao-Biao Lin, and Genevieve approximations of semigroups and partial
behaviour York.
Raugel, differential
Upper semicontinuity equations. Math.
7. A. Friedman and L. Oswald, The blow-up time for higher-order with small leading coefficients. J. Differential Equations 75, 8. K. B. Subrahmanyam, and stability of rotating, 342-352 (1987).
of solutions.
semilinear
In Lecture Notes of attractors parabolic
K. R. V. Kaza, G. V. Brown, and C. Lawrence, pretwisted, preconed blades including Coriolis
Nonlinear vibration effects. J. Aircraft 24,
10. G. Fix and G. Strang, Cliffs, NJ (1973). 11. T. F. Chen, Department
12. L. Cesari,
G. J. Fix, and R. Kannan, of Mathematics, University
Rotor Blade Dynamics Project, Reports of Texas,
Arlington,
equations
239-263 ( 1988).
Helicopter Theory. Princeton Univ. Press Princeton, NJ, ( 1980). An Analysis of the Finite Element Method. Prentice-Hall,
9. W. Johnson,
for
Camp. ( 1987).
and Bell Textron
Englewood 1 and 2. ( 1988).
Asymptotic Behavior and Stability Problems in Ordinary D@erential Equations.
Springer-Verlag,
New York.
13. G. Dahlquist and A. Bjorck, Numerical Analysis. Prentice-Hall, Englewood Cliffs, NJ ( 1978). 14. J. Starch and G. Strang, Paradox lost: Natural boundary conditions in the Ritz-Galerkin method. Internut. J. Numer. Methods. Engrg. 26, 2255-2266 (1988). 15. C. Bernardi
and Y. Maday,
Comput. & Structures 30( l/2),
Spectral methods for the approximation 205-216 ( 1988).
of fourth
problems.