An Alternating Frequency/Time Domain Method for

Calculating the Steady-State Response of Nonlinear

Dynamic Systems
T. M. Cameron, J.H. Griffin

To cite this version:

T. M. Cameron, J.H. Griffin. An Alternating Frequency/Time Domain Method for Calculating
the Steady-State Response of Nonlinear Dynamic Systems. Journal of Applied Mechanics,
American Society of Mechanical Engineers, 1989, .

HAL Id: hal-01333697

https://hal.archives-ouvertes.fr/hal-01333697
Submitted on 18 Jun 2016

HAL is a multi-disciplinary open access Larchive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinee au depot et a la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publies ou non,
lished or not. The documents may come from emanant des etablissements denseignement et de
teaching and research institutions in France or recherche francais ou etrangers, des laboratoires
abroad, or from public or private research centers. publics ou prives.
 T. M. Cameron
An Alternating Frequency/Time
Domain Method for Calculating the J. H. Griffin
Mem. ASME

Steady-State Response of Department of Mechanical Engineering,

Carnegie-Mellon University,

Nonlinear Dynamic Systems Pittsburgh, PA 15213

A method is proposed for analyzing the steady-state response of nonlinear dynamic

systems. The method iterates to obtain the discrete Fourier transform of the system
response, returning to the time domain at each iteration to take advantage of the
ease in evaluating nonlinearities there- rather than analytically describing the
nonlinear terms in the frequency domain. The updated estimates of the nonlinear
terms are transformed back into the frequency domain in order to continue iterating
on the frequency spectrum of the steady-state response. The method is demonstrated
by solving a problem with friction damping in which the excitation has multiple
discrete frequencies.

solution is assumed to have one dominant frequency compo-

1 Introduction nent- or at most two. 1 In situations where it can be applied,
A variety of methods are presently available for calculating and trusted to give accurate solutions (where a single harmonic
the response of nonlinear dynamic systems, (see, for example, description of the nonlinearity is adequate), the HB method is
Iwan (1974) for analytical approaches, and Press et al. (1986) by far the most computationally efficient method for obtain-
for numerical approaches). The method introduced in this ing steady-state solutions to nonlinear dynamics problems. It
paper offers advantages over most other methods in its has been used extensively to study friction damping (e.g.,
generality, efficiency, accuracy, and ease-of-use for Griffin, 1980; Sinha and Griffin, 1983 and 1984; Dowell,
calculating the steady-state response of forced nonlinear 1983; and Jezequel, 1983). Generally, however, this method is
systems. The method involves an iterative scheme that alter- limited to harmonically excited systems and it becomes in-
nates from the frequency to the time domain in order to take creasingly inaccurate as higher harmonics in the expansion of
advantage of the convenience in performing certain tasks in the nonlinear force become significant. A fair amount of work
each domain. Consequently, this Alternating Frequency /Time is also required to expand the nonlinear force, although DF's
(AFT) method easily provides either time trajectory or fre- for a large number of nonlinearities have been tabulated (Gelb
quency spectrum solutions. and VanderVelde, 1968, Appendix B).
The purpose of this paper is to introduce the AFT method An Incremental Harmonic Balance (IHB) method which ac-
and illustrate its utility with a friction damping problem: the commodates multiple harmonic components in the description
response of a single degree-of-freedom (sdf) system with an of the nonlinear force was introduced by Lau et al. in 1982.
elastic/perfectly-plastic friction damper subject to three The method was originally applied to continuous (cubic)
discrete frequencies of excitation is calculated. nonlinearities (Lau et al., 1982, 1983), but Pierre et al. (1985)
extended the method to analyze the rigid/perfectly-plastic fric-
2 Alternative Spectral Methods tion damper problem in which the nonlinear force is discon-
tinuous. The IHB method appears to handle strong
One of the most popular methods for approximating the nonlinearities well, and provides accurate solutions over a
frequency response of nonlinear systems is known as the Har- much broader range than the HB method since it accom-
monic Balance (HB) method in the vibration literature (e.g., modates multiple harmonics. The method is ideally suited to
Nayfeh and Mook, 1979), and as the Describing Function parametric studies since it can step from a known state of
(DF) approach in the control literature (e.g., Gelb and Vander vibration to a neighboring state which corresponds to an in-
Velde, 1968). This is essentially a spectral method in which the cremental change in one of the governing parameters of the

----r:rhe DF approach has been extended to allow two input frequencies in the
case of odd, memoryless (nonhysteretic) nonlinearities (Gelb and VanderVelde,
1968, Chapter 5 and Appendix D), but the nonlinear force and the system
response are described only in terms of the two input frequencies.
system. However, there are also significant disadvantages to time domains. A fast Fourier transform 2 algorithm makes this
the IHB method. cost reasonable, especially in obtaining steady-state solutions
Computationally the IHB method is no more efficient than to periodically excited systems.
numerical integration, at least for the discontinuous friction Consider a nonlinear dynamic system in which the terms
damping problem (cf., Pierre et a!., 1985, "Conclusions"), containing nonlinearities are grouped on the right hand side of
and the computational effort grows geometrically with the the equation with the external forcing term. The equation of
number of harmonics considered. As currently implemented, motion for this system may be represented by:
the frequencies are integer multiples of the excitation frequen- (1)
L(x(t)J =f(t) +Nj(X,X)
cy, which will hide the subharmonic resonances that may oc-
cur in nonlinear systems. Another disadvantage of the IHB where L is the differential operator representing the linear
method is the amount of analytical work required. Functions terms, f(t) is the external forcing term, Nj( ) represents the
in the equation are expanded in a first order Taylor series nonlinear force terms, and x(t) is the desired system response.
(which may be difficult to obtain for discontinuous functions) These are vector functions for multi-degree-of-freedom (mdf)
and substituted into the original equation to obtain a linear systems, or scalar functions for sdf systems. Taking the
"incrementalized" equation. The incrementalized equation is Fourier transform (FT) of both sides yields:
solved by the Galerkin method, which requires further expan- A(w)X(w)=F(w)+1j(w,X(w)) (2)
sions and integrations. Finally, although the approach is
general, the expansions are specialized, so each new dynamic where A(w), the transform of L, is a (matrix) function of the
system requires its own formulation and implementation. frequency, w, F(w) is the transform of f(t), 1j(w,X(w)) is the
Other numerical methods, using techniques similar to those transform of Nj(), and X(w) is the transform of x(t). In the
in the AFT method, describe functions in both the frequency transform domain the nonlinear forces depend only on w and
and time domains. In computational fluid dynamics, "Spec- X(w) because derivatives of x(t) transform into linearly scaled
tral'' and "Pseudo-Spectral" methods (cf., Zang and Hus- forms of X(w): 3
saini, 1985) use Fourier transform techniques to expand func- x(t) #X(w)
tions into (trigonometric or Chebyshev) polynomials. This
x(t)#-jwX(w)
permits derivatives and other operations on those functions to
be performed with higher order accuracy than finite difference x(t) # -w 2 X(w) (3)
methods achieve for similar effort. Cesari (1963) appears to be (" # " denotes a Fourier Transform pair).
the first to obtain steady-state solutions to nonlinear dynamic
equations by using a trigonometric polynomial expansion, Since an FFT will be used to compute and represent the fre-
which is equivalent to a truncated Fourier series or a discrete quency dependent functions, (2) is required to hold at the
Fourier transform. Urabe (1965) expands Cesari's results and discrete frequencies, wk, introduced by the FFT. Thus, (2)
(Urabe and Reiter, 1966) shows how the Fourier coefficients becomes a system of nonlinear algebraic equations 4 :
may be computed at each iteration with a discrete Fourier A(w)X(w) =F(w) + 1j (w,X(w)) (4)
transform, once the nonlinear functions have been evaluated
in the time domain. Urabe identifies this approach as an ap- where
plication of the Galerkin method. Ling and Wu (1987) use a w=(wd (5)
fast Fourier transform (FFT) to implement Urabe's approach, and
calling it the "Fast Galerkin" (FG) method. Although there 27rk N
1 .
are fundamental similarities in the concepts underlying the FG wk k=0, ... , (6)
and AFT methods, there are differences in how the methods ATN
are formulated and implemented. These differences are Here, AT is the sampling period and N is the number of
discussed later. samples in the FFT of each time domain function. In general,
A is a matrix whose components are functions of the indepen-
dent variable (wk), while X, F, and 7l are vector functions of
wk. For example, at a given frequency, say w4 , there will be
3 An Alternating Frequency /Time (AFT) Method X 1 (w 4 ), , Xm(w 4 ) if there are m degrees-of-freedom.
Therefore, the length of X(w) is the product of the number of
3.1 Fundamental Concepts. Frequency domain methods frequencies in the FFT (N/2+ 1) times the number of degrees-
have long been used to analyze linear systems because the fre- of-freedom (m). To minimize the bandwidth of L, the equa-
quency content of a signal is often of greater interest than the tions should be ordered so that X(w) appears as {X1 (w 0 ),
time history and because convolution- a convenient solution X2(wo), . .. 'Xm (wo), Xl(wl), . . . 'Xm (wl), . . . 'Xl(wNn),
method- is more efficiently accomplished in the frequency Xm(wNn )J.
domain . In general, nonlinear forces are nonlinear functions Various iterative methods are available for solving
of the past and present states of the system. Since they are not nonlinear algebraic equations, like (4). One problem facing
linear functions of the states or explicit functions of time, they any method is evaluation of the nonlinear term 1j(w, X (w)). To
cannot be transformed directly into the frequency domain. add it directly to F (w) this term must be known explicitly in
This makes nonlinear systems difficult to analyze in the fre- terms of w, as 7l (w), but this form is difficult to determine in
quency domain, as evidenced by the complicated expansions
required in the frequency domain methods discussed previous- ~Fast Fourier transform" refers to a class of algorithms used to calculate the
ly. Nonlinear forces are easy to evaluate in the time domain, discrete Fourier transform (OFT- cf., Jackson, 1986, Chapter 7; and Op-
penheim and Schafer, 1975, Chapter 3). The notation "FFT" will be used here
however, since the state history of the system is readily to refer to the OFT as well as its implementation. The authors use an FFT
available. This is evidenced by the ease of implementation and subroutine from Press, et al. (1986) .
greater complexity allowed by numerical integration. Th.e
AFT method employs the advantages of both domains by per-
3 The Fourier Transform definition used here

dt and the inverse transform is x(t) =

~ X(w) = j: 00
x(t)fiwt
1/ 2,- j_ 00 X(w)e-Jwtdw . Anal-
forming convolution (as well as filtering and convergence ogous definitions are employed in the FFT.
testing) in the frequency domain and by switching to the time 4 The "w" notation in (4) is for conceptual purposes. For digital implementa-
domain to evaluate the nonlinear forces as functions of the
tion (4) is cast in stricter vector notation- A(k)X(k) =F(k) +11(k)- where the
state trajectory. The cost involved in this technique is the time kth component of each vector represents a sampled version of the corresponding
required to switch back and forth between the frequency arid continuous function at the frequency "'k (ignoring scale factors and errors).
 F(m)-t50--1X(l!l) Equation (4) represents m(1 + N/2) equations (for m

YJ(l!) Converged''
Yes
degrees-of-freedom) which for most methods is unwieldy, ex-
cept for small values of N. The computational effort of the
direct iteration method, O[mN(m + log(N)] multiplications
FFT per iteration, is much less affected by "N' than Newton-
Raphson, which has O[(mN) 2 log(N)] multiplications per
iteration. In the approach used here, rather than limiting N,
which would greatly reduce the number of frequencies con-
sidered in the solution, the direct iteration approach is used
first to estimate which frequency components are significant.
Then, if Newton-Raphson is employed, equation (4) is im-
Fig. 1 Direct iteration implementation of the AFT method plemented only for those values of k for which X(') (wk) is a
significant component in the last estimate of the solution
the frequency domain. However, if N1 ( o) can be expressed as vector obtained from the direct iteration method,
an explicit function of time, N 1 ( t), then it can be transformed (_x{O (w 0 ), , _x(l) (wN12 )). This is discussed further in Sec-
to provide 77(w). To find N 1 (t), the most recent estimate of tion 5. This approach, by solving only the significant frequen-
X(w) is inverse transformed to provide an estimate of x(t). cy components, keeps the number of equations and unknowns
Equation (3)- or the analogous FFT relationships- may be at a minimum while still allowing a full spectrum of frequen-
used to obtain other state variables so that N/ o) may be cies to be considered in the solution.
evaluated as N 1 (t), and transformed to 77(w). For There are two types of error that may corrupt FFT values.
nonlinearities that involve products of time domain functions Aliasing error arises when the FT is nonzero for frequencies
this avoids frequency domain convolution. For history- above the maximum frequency considered by the FFT
dependent nonlinearities, this procedure allows the nonlineari- (wmax = 1r I !J.T). Power at the higher frequencies is aliased or
ty to be evaluated directly from the state history as long as the folded back and added to the FFT values for w < 1r/!J.T.
time period used in the analysis is longer than the ''term of Leakage error occurs when the frequency values of the FFT
memory" of the nonlinearity. With short term memory (k!J.w, or wk) do not coincide precisely with each frequency in
nonlinearities, values of the nonlinear terms at the beginning the FT that has nonzero power. Power at a frequency not
of the period may be obtained from periodicity requirements. represented in the FFT is leaked to frequencies that are
Another problem for nonlinear equation solvers is obtain- represented. Leakage can only be eliminated if the sampled
ing an initial estimate of the solution. An initial estimate time signal is periodic.
(designated .XO(w)) may be obtained from the linear system, ig- For the following sample problem aliasing is insignificant
noring 77( o ), or from a known solution corresponding to a (as verified by taking smaller time steps), and both the excita-
closely related state. For example, consider a sdf system with tion on the system and the system response are periodic, so no
f(t) =A cos(w.t) where a parametric study of the effect of A leakage will occur. Most nonlinear forces are periodic when
on the system response is desired. Once the solution, X(w), is subject to a periodic input, so the assumption of a periodic
obtained for one value of A, it may be used as the initial response to a periodic excitation is typically valid. Discussions
estimate for a neighboring case with a slightly different value of the DFT, FFT algorithms, aliasing, and leakage may be
of A. found in Jackson (1986), Oppenheim and Schafer (1975), and
Press et a!. ( 1986).
3.2 Implementation Considerations. The authors
presently use two methods to solve (4). The first method,
direct iteration, is an iterative implementation of (4):
4 Sample Problem
A(w)o_x{O(w) =F(w) + 17(w,XU- 1 (w)). (7)
In many mechanical applications the most significant source
A(w) 5 and F(w) are fixed by the system and the external forcing of vibrational damping is provided by friction. For example,
function, respectively, once the w/s are selected. New genera- friction damping occurs naturally at joints in fabricated struc-
tions of XUl (w) may be found by evaluating 77( ), as described tures, friction damping is used to control the earthquake
previously, adding it to F(w), and multiplying by A - 1 Since A response of buildings, and friction dampers are used in high
is a constant matrix representing the linear part of the system, speed turbomachinery to reduce resonant stresses in blades.
it may be inverted once and saved. A flow chart of the direct An elastic/perfectly-plastic friction model may also be used to
iteration approach is shown in Fig. 1. analyze the elasto-plastic behavior of materials. Interest in
The second approach used by the authors for solving (4) is these applications has motivated considerable research into
the Newton-Raphson method (cf., Stoer and Bulirsch, 1980), how friction affects the dynamic response of structures (cf.,
with a difference quotient approximation of the Jacobian. For Plunkett's review of friction damping research, 1980).
the example problem in this paper, Newton-Raphson is more In the problem considered the nonlinear element consists of
robust than direct ithation (it converges at times that direct a spring in series with a friction constraint (Fig. 2). The fric-
iteration fails) and it converges in fewer iterations. However, tion contact is assumed to exhibit simple Coulomb-type
when both methods converge, computation time is greater for behavior, i.e., when slip occurs the force is constant in
Newton-Raphson than for direct iteration because of the ef- magnitude and directed to oppose relative motion. Since the
fort required to recompute the Jacobian at each iteration. damper force in this model depends on the history of the
Newton-Raphson is also more difficult to implement. In some response, it is easiest to illustrate the damper force by exam-
cases Newton-Raphson fails to converge, so more robust ple. In Fig. 2 (b), for instance, the damper force exhibits
methods (cf., Wolfe, 1959; Broyden, 1969; Stoer and bilinear hysteresis, which occurs when the input to the damper
Bulirsch, 1980; Shacham, 1986) may be required if a solution is harmonic. For more complicated motion the damper force
is desired for all regions of the parameter space. could take on any value within a similar hysteresis loop, where
the upper and lower limits of the input, which represent the ex-
treme amplitudes of motion, need not be symmetric about
---sp;ofessor Alok Sinha of Pennsylvania State University has pointed out that
in the absence of damping A may be singular. This may be remedied by adding a
zero. In the sample problem the equations have been nor-
linear damping term to both sides of any equation at fault. The damping term malized so that the friction force, 11-C1 , is unity and does not
may be lumped with the nonlinear term(s) on the right-hand side of an equation. appear explicitly.
 Input x(t), x(t) 2.0
R-K
1.5 AFf
1.0
Contact Force, Cr

0.5

x(t) 0.0

-0.5

-1.0
(a)
Output Force -1.5

-2.0
725 730 735 740 745 750 755

time
Fig. 4 Amplitude and phase comparison between AFT method and
fourth-order RungeKutta solutions lor three frequency excitation

103

(b) 102 "

" "
Fig. 2 Elastic/perfectly-plastic friction damper; (a) graphic representa 101
lion and (b) sample output for harmonic input " "
100 " "l!J'l,
f(t)
10- 1
" "
i x(t)

_j
I X(ro) I
10- 2
" ~"
10-3

10" 4

10-5

10"
.I
(0
10 " " 100

Fig. 5 Discrete Fourier transform of response lor three frequency

excitation

Fig. 3 Elastic/perfectlyplastic damper and system

where f(t) =0.75 sin (0.5 1)+0.5 sin(0.75t)+0.375 sin(t), .1.
the fraction of critical viscous damping, is 0.01 (1 percent),
The system used to illustrate the AFT method has only one and E, a measure of the relative spring stiffnesses, is 0.25. The
degree-of-freedom. From the formulation in Section 3.1 it is distance the mass must displace for slip of the friction damper
clear that the method easily extends to multiple degrees-of- to occur is unity. The nonlinear force on the mass is equal to
freedom (mdt). The computational savings of the AFT the stiffness, E, of the elastic friction damper times the
method over numerical integration should be even more displacement of the friction contact point, y( ). The contact
dramatic for mdf systems since the number of frequencies to slips when the magnitude of the friction force equals E; the
consider depends only on the excitation and the nonlinearities, resulting force resists the motion. The specific excitation was
not on the number of degrees-of-freedom. chosen to provide an interesting response.
4.1 An Elastic/Perfectly-Plastic Friction Damper. The One period of the steady-state response of this system is
sample problem was previously unsolved, although it can be shown in Fig. 4. The AFT and numerical integration solutions
solved easily by numerical integration. It would be either dif- compare almost perfectly in both amplitude and phase. Since
ficult or impossible to solve by other methods. The accuracy slip occurs at an amplitude of unity, there is considerable par-
of the solution is judged by comparison to numerical integra- ticipation of the friction damper in this response. The frequen-
tion. The comparison is made in the time domain to show the cy spectrum of the response (FFT) is shown in Fig. 5. This
accuracy of the AFT method in describing both the amplitude problem is used in the following section to make accuracy and
and phase of the steady-state response. efficiency comparisons between the AFT method and fourth-
Multiple discrete frequencies of excitation arise in high order Runge-Kutta integration.
speed turbomachinery applications because as the turbine
blades rotate they pass through a nonuniform, circumferen- 4.2 Convergence, Error, and Computational Efficiency.
tially distributed pressure field. This introduces excitation fre- In the absence of aliasing and leakage the FFT of a function is
quencies that are integer multiples of the rotation speed of the equivalent to a truncated Fourier series representation.
turbomachine. It is also frequently necessary to model friction Therefore, the AFT method should exhibit the convergence
dampers as having compliance, rather than being perfectly and error characteristics of Fourier series. Given a periodic ex-
rigid. These physical considerations motivate the sdf model citation, with period II, and assuming the time sampling is
with an elastic/perfectly-plastic friction damper subjected to over one period, the maximum frequency in the FFT represen-
three discrete frequencies of excitation (Fig. 3). The non- tation is:
dimensionalized equation of motion for this system is: 7r
Wmax=ITN. (9)
x+2tx+(1-E)x=f(t) -Ey(x,x) (8)
 Depending on damping, for example, the system may require
any number of periods of oscillation before steady-state is
reached. Another concern for R-K is that with complicated
wave forms, like the solution to this problem (Fig. 4), it is very
difficult to determine when steady-state has been achieved.
The AFT method tests for convergence in the frequency do-
main where it is easy to discern. However, for a sample com-
parison, the following test was performed. This sample prob-
lem- a sdf system with three input frequencies and an
elastic/perfectly-plastic friction damper- was analyzed
through the 30th period of oscillation. The 30th period was
selected as representative of an average number of cycles to
reach steady-state. In order to achieve the same accuracy (i.e.,
similar error in Fig. 6) the AFT method was implemented in
10 100 1000 10000
single precision with 256 points per period, and the R-K
N method was implemented in double precision with 512 points
Fig. 6 Error convergence comparison between AFT and fourth-order per period. In this case the AFT method took 21.7 seconds 6
R-K methods and Runge-Kutta took 389 seconds on a Macintosh II, a factor
of 18 difference.
Since the highest frequency is proportional to the number of
samples per period, N, a continuous function which exhibits 5 Features of the AFT Method
w- 2 convergence, like A(w)- 1 for second order systems, will
exhibit N- 2 convergence in its FFT components. For the sam- Certain fundamental concepts underlying the AFT method
ple problem, the nonlinear force on the system, N 1 (t), is con- are also employed in the FG (Fast Galerkin) method intro-
tinuous, but has a discontinuous derivative, so 11(w) will also duced by Ling and Wu (1987). Ling and Wu implement
converge O(N- 2 ). Therefore, X(w) is expected to converge Urabe's Galerkin method with an FFT, making the method
O(N- 4 ), which is confirmed by Fig. 5. (F(w) has infinite computationally efficient. In Urabe's approach the system of
order convergence so it does not affect the convergence of nonlinear equations is put into state variable form, a solution
X(w).) is asssumed in the form of a truncated Fourier series, and the
In order to compare convergence of the time domain errors, coefficients of the Fourier series are solved for
as N increases, for the AFT and Runge-Kutta methods the er- iteratively- with nonlinear terms evaluated in the time do-
ror is assumed to be of the form: main. Since the coefficients of the Fourier series are equivalent
to the components of the DFT, the AFT and FG methods are
EN=XN-Xex =aNc (10) similar in this respect. There are, however, differences in how
where xN is the value of x( t) at a specified time using N sample the methods formulate and implement solutions of equations.
points, Xex is an extrapolated "exact" answer, EN is the These differences highlight some of the features of AFT
estimated error, N is the number of sample points per period, method.
"a" is a constant, and "c" is the constant that gives the order The first difference between the AFT and Galer kin methods
of convergence of the error. Using three values of N gives lies in the form of the equations to be solved. The AFT
three equations that allow the parameters Xex a, and c to be method permits the equations to be of any order, while
evaluated. The consistency of the error convergence is checked Urabe's formulation, which Ling and Wu repeat, assumes the
with two of the original three values of N, plus and an addi- equations are in first order state variable form. The AFT
tional, fourth value. These permit another calculation of the method only finds the linearly independent functions, thereby
parameters in (10). reducing the number of equations to be solved.
For the AFT method, the O(N- 4 ) convergence of the FFT Another feature of the AFT formulation is the general man-
should produce similar convergence in the time domain error. ner in which it handles nonlinearities. In the problems con-
This is the case if the FFT components are known precisely. sidered by Urabe, as well as Ling and Wu, the nonlinearities
However, changing N affects truncation and aliasing errors, are expressed analytically. In many applications, such as the
so for different values of N the FFT components for the same elastic-plastic damper problem, the nonlinearities are
frequency vary slightly. Consequently, additional errors ap- hysteretic and analytical expressions may be difficult or im-
pear in the time domain. Using N = 64, 128, and 256, in single possible to obtain. In the AFT approach the user is provided
precision, the AFT time domain error for equation (10) was with the full state history in order to compute the nonlinear
found to be O(N- 3 1). Since the sampling period, !::..T, is II/N, forces in a user supplied subprogram, in a manner similar to
this is equivalent to O(!::..Tl 1) convergence in the error. The er- that used in standard time integration codes.
ror was consistent when estimated with N= 128, 256 and 512. Perhaps the most significant computational advantage of
When the functioli.s in a problem are infinitely differen- the AFT method is its technique for reducing the number of
tiable a fourth-order Runge-Kutta (R-K) integration has error equations. As discussed previously, increasing the number of
O(!::..T'), or O(N- 4 ). When a function has discontinuous frequencies included in the solution can significantly increase
derivatives, like Ey(t) (the nonlinear force in problem 4.1), that
error estimate is invalid. Since the R-K method steps pro- ~would also be possible to solve this problem using Runge-Kutta time in-
gressively forward in time, accumulated roundoff error can tegration in conjunction with a two point boundary value formulation. In this
also dominate truncation error. In fact, for this problem, no approach the displacement and velocity of the mass at the beginning of the
consistent error convergence (using (10)) was observed for R-K period are unknowns whose values must be determined so that the solution
satisfies the periodicity requirements. Such an approach requires an iterative
in single precision, with N = 64, 128, 256, 512, and 1024. Using method for solving the resulting nonlinear constraint equations (for a discussion
the same values of N in double precision R-K simulations, the of this approach see, for example, Press et al., 1986, Chap. 16). This procedure
error was O(!::..TL 6 ), or O(N-1. 6 ). The error results for both is less standard than time integration and, consequently, was not investigated
R-K and AFT methods are plotted in Fig. 6. here. However, it could be significantly more efficient than time integration and
It is impossible to make a general comparison of computa- more competitive with the computational efficiency of the AFT method. Again,
the outcome of such a comparison would depend on the specific physical system
tional efficiency between the AFT and R-K methods. The effi- considered, as well as on the complexity of the frequency content of the
ciencies of both depend strongly on system parameters. excitation.
the computational effort, but in many problems it may be Cesari, L, 1963, "Functional Analysis and Periodic Solutions of Nonlinear
necessary to consider both low and high frequencies. The AFT Differential Equations," Contributions to Differential Equations, Vol. I, pp.
149-187.
implementation permits the .user to include a relatively large Dowell, E. H., 1983, "The Behavior of a Linear, Damped Modal System with
number of frequencies and to set a "cutoff" fraction. While a Nonlinear Spring-Mass-Dry Friction Damper System Attached," Journal of
using direct iteration (which is computationally efficient), the Sound and Vibration, Vol. 89, No. I, pp, 65-84.
algorithm scans all the frequency components to estimate Gelb, A., and VanderVelde, W. E., 1968; Multiple-Input Describing Func-
tions and Nonlinear System Design, McGraw-Hill Book Co., New York.
which components lie above the cutoff fraction relative to the Griffin, J, H., 1980, "Friction Damping of Resonant Stresses in Gas Turbine
largest component. These components are considered ~'signifi Engine Airfoils," ASME Journal of Engineeringjor Power, Vol. 102, pp.
cant.'' If it becomes necessary to employ an intensive equation 329-333.
solver, like Newton-Raphson, only the equations corre- !wan, W. D., 1974, "Application of Nonlinear Analysis Techniques," Ap-
sponding to the significant components are included. This ap- plied Mechanics in Earthquake Engineering, W, D. !wan, ed,, ASME AMD
Vol. 8,
proach also makes it easy for the user to examine tradeoffs Jackson, L. B., 1986, Digital Filters and Signal Processing, Kluwer Academic
between higher accuracy (low cutoff) and lower computation Publishers, Boston, Chapter 7.
time (high cutoff). Jezequel, L, 1983, "Structural Damping by Slip in Joints," ASME Journal
of Vibration, Acoustics, Stress, and Reliability in Design, Vol. 105, No. 4, pp.
6 Concluding Remarks 497-504.
Lau, S. L., Cheung, Y. K., and Wu, S. Y., 1982, "A Variable Parameter In-
In the absence of aliasing, leakage, and roundoff errors the cremental Method for Dynamic Instability of Linear and Nonlinear Elastic
AFT method converges uniformly to the exact answer as N - Systems," ASME JouRNAL oF APPLIED MECHANics, Vol. 49, No. 4, pp.
849-853.
oo. Even when it cannot be eliminated, means are available for Lau, S. L., Cheung, Y. K., and Wu, S. Y., 1983, "Incremental Harmonic
limiting the error so that good approximate solutions may be Balance Method with Multiple Time Scales for Aperiodic Vibration of
obtained. For example, in many problems the power in a Nonlinear Systems," ASME JouRNAL OF APPLIED MECHANICS, Vol. 50, No.4,
signal decays asymptotically to zero as the frequency in- pp, 871-876.
Ling, F. H., and Wu, X. X., 1987, "Fast Galerkin Method and Its Applica-
creases, so aliasing cannot be eliminated (without analog tion to Determine Periodic Solutions of Non-Linear Oscillators," International
filtering). However, it is still possible to obtain exact solutions Journal of Non-Linear Mechanics, Vol. 22, No. 2, pp. 89-98.
to within the precision of the computer by selecting il T so that Nayfeh, A. H., and Mook, D. T., 1979, Nonlinear Oscillations, John Wiley
power at frequencies exceeding (1r/ ilT) cannot numerically af- and Sons, New York.
fect the solution. In the sample problem the error in the AFT Oppenheim, A. V., and Schafer, R. W., 1975, Digital Signal Processing,
Prentice-Hall, Inc., Englewood Cliffs, N.J., Chapters 1, 3, and 6.
method dropped off much more quickly than the error in Pierre, C., Ferri, A. A., and Dowell, E. H., 1985, "Multi-Harmonic Analysis
numerical integration. of Dry Friction Damped Systems Using an Incremental Harmonic Balance
Solutions generated by the AFT method are conveniently Method, ASME JOURNAL OF APPLIED MECHANICS, Vol. 52, No.4, pp. 958-964.
presented in either the frequency or time domain, or both. The Plunkett, R., 1980, "Friction Damping," DampingApplicationsfor Vibra-
AFT method also lends itself to parametric and sensitivity tion Control, AMD Vol. 38, P. J. Torvik, ed., ASME, New York.
Press, W. H., Flannery, B. P., Teukolsky, S, A., and Vetterling, W. T., 1986,
studies because the solution for one set of parameters supplies Numerical Recipes, Cambridge University Press, New York, Chapters 9 and 12.
a good estimate of the solution for a "neighboring" set of Shacham, M., 1986, "Numerical Solution of Constrained Nonlinear
parameters, which aids in convergence of the iterations. Algebraic Equations," International Journal for Numerical Methods in
Relative to other frequency domain methods, the AFT Engineering, Vol. 23, pp. 1455-1481.
Sinha, A., and Griffin, J. H., 1983, "Friction Damping of Flutter in Gas Tur-
method is amenable to more general, complicated bine Engine Airfoils," AIAA Journal of Aircraft, Vol. 20, No.4, pp. 372-376.
nonlinearities and excitations. By switching to the time do- Sinha, A., and Griffin, J. H., 1984, "Effects of Static Friction on the Forced
main to evaluate the nonlinear forces, the AFT method can Response of Frictionally Damped Turbine Blades," ASME Journal of
use the same description of the nonlinear forces used for Engineering for Gas Turbines and Power, Vol. 106, No. 1, pp. 65-69.
numerical integration. The only constraint on the excitation is Stoer, J., and Bulirsch, R., 1980, Introduction to Numerical Analysis,
Springer-Verlag, New York, Chapter 5.
that its Fourier transform can be approximated using the FFT Urabe, M., 1965, "Galerkin's Procedure for Nonlinear Periodic Systems,"
approach. Archive for Rational Mechanics and Analysis, Vol. 20, pp, 120-152.
Urabe, M., and Reiter, A., 1966, "Numerical Computation of Nonlinear
Acknowledgment Forced Oscillations by Galerkin's Procedure," Journal of Mathematical
Analysis and Applications, Vol. 14, pp, 107-140.
This work was supported by the National Science Founda- Wolfe, P., 1959, "The Secant Method for Simultaneous Nonlinear Equa-
tion, grant number MSM 8703887. tions," Communications of the ACM, Vol. 2, pp. 12-13,
Zang, T. A., and Hnssaini, M. Y., 1985, "Recent Applications of Spectral
References Methods in Fluid Dynamics," Lectures in Applied Mathematics, Vol. 22: Large-
Scale Computations in Fluid Mechanics, edited by B. E. Engquist, S. Osher, and
Broyden, C, G., 1969, "A New Method of Solving Nonlinear Simultaneous R. C. J. Somerville; eds., American Mathematical Society. Providence, R.I.,
Equations," Computer Journal, Vol. 12, pp, 94-99. Part 2, pp. 379-409.