NBNotation
complex variable
complex variable
dimensionless composite parameter
complex variable
error function
imaginary error function
linear sweep rate in Hz per minute
nonzero start frequency for an octave sweep rate
natural frequency in Hz
complex variable
complex variable
octave sweep rate in octaves per minute
complex variable
time (also used as a dummy integration variable)
upper limit of search for peak values
time at which instantaneous frequency of excitation for linear sweep equals 
time at which instantaneous frequency of excitation for octave sweep equals 
new independent variable for octave sweep
value at which instantaneous frequency of excitation for octave sweep equals 
single-degree-of-freedom system displacement response
single-degree-of-freedom system velocity response
single-degree-of-freedom system acceleration response
initial displacement
initial velocity
complex variable
octave sweep rate in octaves per second
composite parameter for closed-form solution for linear sweep
composite parameter for closed-form solution for linear sweep
composite parameter for closed-form solution for linear sweep
composite parameter for closed-form solution for octave sweep
general phase function
composite parameter for closed-form solution for octave sweep
initial phase value
incomplete gamma function
dummy integration variable
composite variable proportional to 
composite parameter for closed-form solution for linear sweep
composite parameter for closed-form solution for linear sweep
natural frequency in radians per second
linear sweep rate in radians per second per second
generalized sweep forcing function
multiplication
composite parameter for closed-form solution for octave sweep
composite parameter for closed-form solution for octave sweep
critical damping ratio
dummy integration variable
1. Introduction
Harmonic excitation is a fact of life in systems with rotating machinery, such as liquid rocket engine turbopumps, spacecraft momentum wheels, aircraft turbojet engines, electric plant steam turbines and liquid-transport turbine compressor trains. Associated with high performance are high shaft speeds and the resulting excitation caused by imbalances in the rotating components and imperfections in the shafts and ball bearings. Furthermore, phenomena such as shaft whirl and rotor dynamic instability are critical design aspects. Although performance requirements dictate design parameters such as shaft speed, avoiding certain speeds due to dynamic interactions within the system is also a critical design consideration. Completely avoiding critical speeds may not be possible. For example, if the critical speeds are below the operational shaft speed, then at startup and shutdown, the rotation rate sweeps through them. The magnitude of the response is a function of the sweep rate, system damping and modal gains at the excitation and response locations. In addition, bearing imperfection can produce excitation above and below the operational frequency, and responses to these imperfections are also a function of the sweep rate associated with the startup and shutdown of the system. In addition to rotating machinery considerations, frequency sweep effects are a critical aspect of harmonic base shake vibration testing, as employed in the aerospace industry, for example. Therefore, it has been recognized that being able to predict the vibration response of systems to swept-frequency excitation is critical (e.g. [1–7]).
In 1932 Lewis presented the first response of a single-degree-of-freedom system to linear frequency sweep excitation [1]. He derived an expression for the envelope functions that contained the peak values. The limited quantitative results presented by Lewis were obtained by graphical integration for various levels of damping and sweep rate. Lewis concluded that the greater the sweep rate, the larger the attenuation relative to steady-state response, and the higher the instantaneous frequency of excitation would be at which the peak envelope response occurs. Fearn developed in 1967 [2] an algebraic expression for the time at which the peak displacement response of a single-degree-of-freedom system subjected to a linear frequency sweep would occur, and an approximate magnitude of the displacement response. Until Cronin’s dissertation [3], published in 1965, analytical studies were generally restricted to linear frequency sweep, and exponential sweep-excitation studies were mostly experimental in nature. Cronin did provide results for relatively slow sweep rates; his work included analog studies involving linear and exponential excitation frequency sweeps. In addition to spring-mass single-degree-of-freedom systems, work has also been done on unbalanced flexible rotors whose spin rate swept through its critical speeds, e.g., [4]. In these types of systems the modes of vibration would be a function of the spin rate and the resulting gyroscopic moments. In 1964 Hawkes [5] described an approach for obtaining the envelope function of the response of single-degree-of-freedom systems subjected to octave sweep rates. He credits the solution approach to an unpublished document written in April 1961 by T. J. Harvey. From the publication, it is unclear how all required initial conditions were obtained for the resulting differential equations that were solved by numerical integration. The results, however, are consistent with subsequent work published by Lollock [6], who extended the work for both linear and octave sweep rates to useful damping and natural frequency ranges.
In approaches where the envelope function is used to identify the peak response, several factors need to be considered. First, the peak of the envelope function may not coincide with the peak of the time history response; this could lead to an incorrect estimate of the instantaneous excitation frequency that coincides with the peak response. The discrepancy would be greatest for low-frequency systems and decrease as the natural frequency increases relative to the starting frequency of the sweep. Another peculiar feature of this approach is that, whereas the original equation of motion is a second-order differential equation with two initial conditions (say, on the function and its derivative), the envelope equations turn out to be two coupled second-order differential equations, each of which requires two initial conditions, and there does not appear to be any way to derive these four necessary initial conditions from the original two for the equation of motion. There are physical arguments that one could make regarding what the initial conditions ought to be, but there does not appear to be any way to mathematically derive them from the original initial conditions.
It is the purpose of this article to extend and complement previously published work by proposing explicit closed-form solutions to both linear and octave frequency-sweep excitation. This allows the computation of the peak response, not just the peak of the envelope function. The closed-form solutions involve error functions and incomplete gamma functions of complex arguments, computations of which require numerical precision exceeding that which today’s computers can provide. The approach used to overcome this will be described. The closed-form solutions are compared to solutions obtained by numerical integration of the equations of motion. Having the ability to compute closed-form solutions, studies were performed to explore the impact of the frequency separation between the start frequency of the sweeps and the natural frequency of the system. In addition, results are presented showing the fine structure of the peak response in relation to the steady-state resonance response as a function of natural frequencies and critical damping ratios. This includes some unexpected results, in that the peak response curves exhibit highly nonlinear behavior with discontinuities in the derivative.
2. Equations of Motion
The differential equation for the motion of a single-degree-of-freedom system driven by harmonic excitation with a linear frequency sweep is given by
 |
(1) |
where 
is the critical damping ratio, 
is the natural frequency, 
is the sweep rate in radians per second per second, and the dots indicate differentiation with respect to time. Assume, without any loss of generality, a sweep starting frequency of zero, a force magnitude equal to the mass of the system and initial conditions of 
and 
. The differential equation of motion of a single-degree-of-freedom system driven by harmonic excitation with an octave frequency sweep is
 |
(2) |
where 
, 
is the octave sweep rate in octaves/sec, and 
is the nonzero start frequency of the sweep. As for the linear sweep case, assume a force magnitude equal to the mass of the system and initial conditions of 
and 
.
It is helpful to also write both the linear sweep and the octave sweep equations in the following form:
 |
(3) |
where 
is a general phase function and 
is the initial phase. Both the linear and octave sweep equations of motion can be put into the following more general form, which will be useful for constructing closed-form solutions:
 |
(4) |
3. Closed-Form Solution: Linear Sweep
The solution to equation (4) can be expressed as
 |
(5) |
For linear sweep, this becomes
 |
(6) |
If the sine terms are expanded in terms of complex exponentials, then the resulting integrals can be computed in terms of the error function, 
, and the imaginary error function, 
, each with complex argument 
, where 
. Conceptually, the process proceeds as follows:
- After converting the sine terms to complex exponentials, expand out the products of sums of exponentials, splitting the integral accordingly into a sum of several integrals of exponentials and pulling the parts of each integrand that do not depend on the integration variable outside the integral; the resulting integrals will all have the form
. - With some algebraic manipulation, these integrals can be put into the form
or 
, where 
, 
and 
are, in general, complex valued. - Choosing
as the new integration variable, the first of these integrals becomes:
An identical procedure can be applied to the second of these integrals, leading to an expression involving imaginary error functions. Performing the indicated calculations (including the associated algebra) gives the following closed-form solution for the linear frequency-sweep excitation case. In the interests of compactness, it is helpful to first introduce the following auxiliary parameters:
Then the closed-form solution for the linear sweep case can be written as
 |
(7) |
In order to verify that this equation for 
does in fact satisfy the equation of motion, we make use of the fact that the derivatives of the error function and the imaginary error function are given by the exponentials 
and 
. Then substituting equation (7) and its derivatives into the equation of motion yields an expression involving all of the original 
and 
functions, plus a number of terms that do not contain any error functions. Collecting terms with respect to the various error functions, which is relatively straightforward although algebraically tedious, verifies that the coefficient of each of the error functions is zero, and that the terms that do not contain any error functions sum to 
, which is the forcing function on the right-hand side of the equation. Since we are interested in the peak acceleration response, the second derivative of the solution, equation (7), is the sought-after response time history.
4. Closed-Form Solution: Octave Sweep
For the case of octave sweep, it is helpful to make a change of independent variable in equation (2) and let 
, where 
is the octave sweep rate in octaves per minute. With this change of variable, the equation of motion for octave sweep in the 
domain becomes
 |
(8) |
The initial conditions become 
and 
. Similarly, the expression for the second derivative with respect to time becomes (in terms of derivatives with respect to 
),
 |
(9) |
where we define 
. The advantage of making this change of variable, from the perspective of numerical integration, is the absence of exponential functions of time in the forcing function in equation (8); rather, the forcing function is a constant-frequency sine wave, and the coefficients in the equation are at most quadratic in 
. This greatly improves the stability and reliability of the numerical integration.
It is helpful to write equation (8) for the octave sweep in the following more general form:
 |
(10) |
where 
. Then using the variation of parameters method, we obtain the following expression for 
:
 |
(11) |
Substituting 
for 
and then expanding the sines in terms of complex exponentials yields integrals of the form 
, which are readily expressed in terms of incomplete gamma functions after algebraic transformation. The incomplete gamma function is given by 
. For compactness, it is helpful to first introduce the following auxiliary parameters: 
, 
, 
, 
and 
. Then the resulting expression for 
reduces to
 |
(12) |
Substituting 
yields the corresponding solution in the time domain:
 |
(13) |
Computing the first and second derivatives of equation (13) and substituting them into the original equation of motion, equation (2), one discovers, after some algebra and collecting terms with respect to the various incomplete gamma functions, that the resulting equation can be put into the form 
. Since we are interested in oscillatory motion, which implies 
, it follows that 
reduces to zero, thereby showing that equation (13) does indeed satisfy equation (2).
5. Challenges in Separating Real and Imaginary Parts of Closed-Form Solutions
The sought-after solutions are the real parts of equation (7) and equation (13). For the linear sweep, series expressions exist for the real and imaginary parts of both 
and 
: functions.wolfram.com/GammaBetaErf/Erf/19 for 
and functions.wolfram.com/GammaBetaErf/Erfi/19 for 
contain series expressions in terms of Hermite polynomials as well as hypergeometric functions. In practice, these series have very slow and highly nonmonotonic convergence properties, with the partial sums fluctuating over many orders of magnitude as successive terms are added. Furthermore, numerical evaluation of these partial sums using exact numbers as inputs is extremely slow and computation time increases nonlinearly with the number of terms, while evaluation using finite-precision numbers yields erroneous results. Since one does not know ahead of time how many terms will be needed for an accurate computation, this approach is impractical. As with the error function, there are similar numerical challenges in computing the incomplete gamma function of complex arguments. Accordingly, the closed-form solutions will be computed using equations (7) and (13) directly.
6. Challenges in Numerical Evaluation of the Exact Closed-Form Solutions
There are also numerical challenges associated with the exact solutions because of the complex arguments of the error and gamma functions. Recall that the error function is given by 
and observe that the magnitude (i.e. the absolute value) of 
is the same as the magnitude of 
, since 
. However, once the argument 
becomes complex, we would need to integrate expressions of the form 
, and the presence of the 
term in the exponent means that the real part of the exponent grows very quickly with 
, that is, as 
. Since 
is analytic in the complex plane, we can use the Cauchy integral theorem for line integrals [8] to break the integral from 0 to (complex) 
into two parts: the integral from 
to 
plus the integral from 
to 
. In the integral from 
to 
, we are in effect integrating 
from 
to 
. Thus, both 
and 
increase very quickly, as shown in the plots in Figures 1 and 2. In order for the end result of the combinations of 
and 
that appear in the exact solution to sum to an oscillatory function, very precise cancellations are needed, meaning that extremely high precision is needed in order to do the numerical evaluations correctly.
Figure 1. Plot of 
. Observe that over the range 
and 
, the magnitude of 
increases to about 
.
Figure 2. Plot of 
. The behavior of 
is similar to that of 
.
Because of the extremely high numerical precision requirements, Mathematica, which implements arbitrary-precision arithmetic, was chosen to compute the closed-form solutions. This made it possible to experiment with different levels of computational precision. Some results were computed with hundreds or thousands of digits of precision. Depending on the values of the input parameters (sweep rate, natural frequency, damping coefficient, etc.), it was found that different levels of precision were needed in order to get reliable results—not a very attractive idea, since it is impossible to know ahead of time how much precision would be needed for any particular set of inputs. Fortunately, Mathematica also allows exact arithmetic (using rational and/or exact symbolic numbers as inputs), and this made it possible to use the exact analytic solutions in a computationally tractable form. More specifically, one can evaluate functions numerically using exact arithmetic by means of the following steps:
- Convert all of the inputs to integers, rational numbers or exact symbolic numbers such as
or 
, or rational multiples thereof, all of which are treated as having infinite precision. - Set the global variable
, which specifies the maximum number of extra digits of precision, to 
. This enables as much extra precision as possible. - Evaluate the function of interest with the desired (exact) inputs. This will, in general, yield a very complicated exact expression.
- Evaluate this exact value to the desired number of digits of precision for the output in order to get a recognizable numerical value, with the understanding that any imaginary “dust” arising from this numerical truncation will be ignored. (For the results presented later, we used 30 digits of output precision.)
7. Comparison of Exact and Numerical Solutions
To build confidence in the closed-form solution, the equation of motion was also solved by direct numerical integration. For the linear sweep case, the results presented herein were obtained from the closed-form solution, equation (7), as well as direct integration of the differential equation of motion, equation (1). The closed-form solution was evaluated by first rationalizing all of the inputs to (7) (other than integers, rational numbers and multiples of 
and 
) using 
(which converts any number to rational form), and then evaluating the real part of the result (to eliminate any very small imaginary numbers) to the desired number of digits of precision (typically 30, 50 or 100) with the 
function. The numerical solution was obtained by integrating the equation of motion (1) with 
out to a some desired maximum time (typically some time after the sweep frequency hits the natural frequency of the system), with 
, 
, and 
.
Figure 3 shows the response time histories for a system with a natural frequency of 5 Hz and a critical damping ratio of 1%. The sweep frequency was started at zero Hz and the sweep rate was 150 Hz/min, or 
. In the figure, the dashed orange line is the closed-form solution and the dotted blue line is the direct numerical integration solution. Clearly, the differences are imperceptible. Table 1 shows the numerical values for both solutions for a randomly selected subset of the time points used in plotting Figure 3. Again, it is evident that for all practical purposes, the solutions are identical.
Figure 3. Acceleration response time histories of a single-degree-of-freedom system, 
, excited by a harmonic force with a linear sweep rate frequency of 150 Hz/min.
For the octave sweep case, the results in Table 2 were obtained from the closed-form solution, equation (13), as well as by direct integration of the differential equation of motion in the 
-domain, equation (8). The procedure for evaluating the closed-form solution for the octave sweep case was identical to that described for the closed-form solution in the linear sweep case. The numerical solution was obtained by integrating equation (8) in the 
-domain with 
from 
out to some desired maximum value of 
(typically corresponding to some time beyond the time at which the sweep frequency hits the natural frequency of the system), with 
and 
, and then using equation (9) to transform the acceleration back to the 
-domain. Figure 4 shows the response time histories for a system with a circular natural frequency of 1/4 Hz and critical damping ratio of 0.01. The sweep frequency was started at 1/8 Hz and the sweep rate was 1/2 octaves/min. The orange dashed line is the closed-form solution and the blue dotted line is the direct numerical integration solution. Again, the differences are imperceptible. Table 2 provides the numerical values for both solutions for a randomly selected subset of the time points used in plotting Figure 4; for all practical purposes, the results are identical.
Figure 4. Acceleration response time histories of a single-degree-of-freedom system, 
, excited by a harmonic force with an octave sweep rate frequency of 0.5 octaves/min.
8. Construction of Peak Response Curves
The construction of the peak response curves involved two steps. First, the times at which the peak of the absolute value of the acceleration occurred were obtained via numerical integration for the desired combinations of 
, 
and 
for linear sweep or 
for octave sweep. These times were then used as the starting points for a very fine-grained search of the exact analytical solutions in order to determine the peak acceleration in each case. Development of a generic algorithm to accomplish this was not trivial, as will be discussed. However, the effort was made easier by previously published results that indicate that the peak envelope values, which would contain the peak response values, would occur after the instantaneous frequency of excitation was equal to the natural frequency of the system. Hence, the search for the peaks was started at the point in the response time history where the instantaneous frequency of excitation was equal to the circular natural frequency of the system. For the linear sweep excitation, the time was computed as
 |
(14) |
and for the octave sweep excitation, the value was computed as
 |
(15) |
In the case of the numerical approach, we sorted the list of computed acceleration values generated via integration, starting at 
in order to find an initial approximation to the peak acceleration, and then did a more refined local search around this peak using standard local optimization techniques. In the case of the analytical approach, much smaller increments in 
were used in order to get a sharper picture of some unusual phenomena that emerge at low frequencies. Accordingly, interpolations were generated of the times at which the numerically generated peak responses occurred, as a function of 
for combinations of 
and 
for linear sweep, or 
for octave sweep. Thus, for any value of 
, we could use this interpolated time value as the starting point for a refined numerical search that involved evaluating the exact analytical solution at very closely spaced time points in a neighborhood of this time. For this, we chose time points that were equally spaced in the phase of the forcing function, that is, at 0.25° phase increments, which provided precise, although computationally intense results. In addition, care was taken to search for the peak sufficiently past the start of the search, given by equation (14) or equation (15), to guarantee that the global maximum peak had been found.
Associated with the question of at which point in a time history to start the search for the peak value, that is, 
and 
, is the question of how far past this point the search should be conducted to guarantee that the global peak has been identified. Unfortunately, the only way we found to reliably accomplish this was through trial and error. For linear sweep, we found experimentally that it was very helpful to divide the 
range into two parts: 
and 
.
For relatively low natural frequencies, that is 
, it was found experimentally that evaluating the function out to 
gives reliable results in most cases, with the peak response typically occurring about 20% of the way out to 
. At low values of 
, however, sometimes the peak response occurred about 45 to 50% of the way out to 
. Although with hindsight, we could have obtained the peak response without going out so far in time, we wanted to be sure that the peak response found was in fact the true global peak response. We observed that in some cases, what looked like a global peak value eventually got “dethroned” by a peak that occurred quite a few cycles later, due to the beating of the frequencies involved. Thus, all of the low-frequency responses, as well as a subset of the high-frequency responses, were visually monitored graphically, and if any peak responses were found at times more than 50% or so of the way out to 
, then the coefficient of 
for 
was increased accordingly. For higher natural frequencies, that is, 
, it was found that 
generally gave reliable results for high sweep rates (~150-200 Hz/min), while 
gave reliable results for lower sweep rates (~10-20 Hz/minute).
In view of the oscillatory nature of the system, it was important to constrain the maximum integration step size to be at most a small fraction of a cycle. Based on prior experience with similar computations, we chose the maximum step size to be 1/40 of a cycle of the largest frequency of interest, which was the sweep frequency at the 
value described previously. For simplicity, we deliberately chose to constrain the maximum step size based on the largest frequency of interest, encountered at time 
, rather than attempt to change the maximum step size as the frequency changed. For the low sweep frequencies encountered in the early parts of a sweep, this step size was much smaller than 1/40 of a cycle, but this did not create any problems. The numerical integrator used (
) employs an adaptive algorithm that adjusts the step size as needed, subject to any user-prescribed constraints. In addition, we used fifth-order interpolation in the numerical integrator so that the acceleration would be a third-order interpolating function. Finally, in view of the progressive increase of sweep frequency with time, we found it useful to specify a maximum of 100,000,000 integration time steps (considerably more than the integrator’s default value), as in some cases a smaller maximum number of time steps (such as 10,000,000) did not allow the adaptive integrator to reach the global peak response.
It was also required that the closed-form solution be evaluated at very closely spaced time increments in order to reliably find the peak acceleration. This strategy leveraged off of the previously computed numerical solutions, that is, the times at which the numerically obtained peak values occurred, in order to do a very fine-grained search (with the closed-form solution) in the neighborhood of the numerically computed peak value. Although a global list of search points could have been generated in other ways without the use of the numerical solution, using the points generated by the numerical integrator seemed like the most efficient approach. The strategy then was to use the numerically generated estimate of when the peak acceleration occurs and search within plus or minus some number of cycles of this time, at equally spaced increments in the forcing function phase. We found that searching within ±60 cycles with 1,440 phase increments per cycle (i.e. at 0.25° phase increments) yielded reliable results.
9. Peak Response Curves for Linear Sweep
Figure 5 shows the peak response (from the exact solution) normalized by 
, the steady-state resonant response when the excitation frequency is equal to the undamped natural frequency of the system, plotted against the natural frequency of the system for three linear excitation sweep rates, 
, 
and 
Hz/min. The system has a critical damping ratio of 
and its natural frequency was varied from 0.25 Hz to 10 Hz in steps of 0.01 Hz. Each of the (almost 1,000) peak response values on each of these curves was computed via the process for computing peak acceleration (from the exact solution) described in Sections 7 and 8, that is, searching within ±60 cycles of the numerically generated estimate of when the peak acceleration occurs, with 1,440 phase increments per cycle. As can be seen, the attenuation of the peak response relative to the resonant steady-state response is significant for systems with low natural frequencies. As the natural frequency increases, which allows a greater number of response cycles during any given excitation frequency range, the attenuation decreases. These results are consistent with those published by others [6]. What is not consistent is the scalloped behavior of the peak curves at the lower frequencies. This behavior was obtained with both the numerically integrated results and the closed-form solution. Figure 6 shows an expanded close-up view of the lower-frequency range of Figure 5 and was generated by simply changing the horizontal plot range in Figure 5. The details visible in Figure 6 will be discussed in more detail later.
Figure 5. Normalized peak response plotted against natural frequency for several linear excitation sweep rates. Left-to-right curves correspond to top to bottom in key.
Figure 6. Close-up of the low-frequency range of Figure 5. Left-to-right curves correspond to top to bottom in key.
Another observation is that the peak response during a frequency sweep can exceed the steady-state resonant response. This is shown in Figure 7, where the normalized peak responses are shown for two sweep rates (Figure 7 was obtained from Figure 5 by simply adjusting the vertical plot range to focus on the overshoot portion of the response). This might seem counterintuitive, since the frequency of excitation is sweeping through the natural frequency and therefore does not dwell. However, the sweep causes a response that is at the natural frequency of the system and that decays as a function of the system damping. Once the sweep frequency passes the natural frequency, the total response is the response due to the excitation plus the free-decay response of the system at its natural frequency. This is what causes the beating in the response once the sweep frequency passes the natural frequency. The decaying free response plus the transient response to the swept excitation can combine to produce higher peak responses than the resonant response caused by harmonic dwell at the natural frequency. The overshoot observed here is consistent with the overshoot observed by Cronin [3].
Figure 7. Close-up of overshoot phenomenon observed in Figure 5. Left-to-right curves correspond to top to bottom in key.
Figure 8 shows the normalized peak response for various sweep rates plotted against the natural frequency squared divided by the linear sweep rate; this normalization allows comparison to results presented in the literature. The critical damping ratio for this system is 
. The data used in Figure 8 is the same as the data used in Figure 5, only plotted differently. Observe that the curves merge into one, as explained by Hawkes [5].
Figure 8. Normalized peak response for several linear sweep rates plotted against 
, where 
is the natural frequency and 
is the sweep rate.
9.1 Discontinuities in Derivative of Peak Response Curves at Low Frequencies
In Figures 5 through 8, one observes periodic discontinuities in the derivative of the peak response curve. Moreover, the curve does not increase monotonically; sometimes it starts to dip down before hitting a discontinuity in slope and resuming its upward trend. One also observes that at very low frequencies the discontinuities in the derivative are not very regular, but as the natural frequency is gradually increased, they take on a much more regular nature. These discontinuities are best understood in terms of what we will call the “competing peaks” phenomenon, which can be most clearly explained by taking several observations into account:
- The peak response always occurs some time after the sweep frequency reaches the natural frequency of the system.
- As the natural frequency of the system is increased, the time at which the sweep frequency reaches the natural frequency occurs later and later, since for these problems the sweep frequency always started at 1/8 Hz.
- Thus the time at which the peak acceleration occurs can be expected to increase as the natural frequency is increased.
- In the array of plots shown in Figure 9, which show the response time histories for several very closely spaced values of
, one observes that as the time at which the peak acceleration is reached increases, the “dominant” peak (i.e. the largest global peak) is eventually overtaken (from one value of 
to the next, i.e. from one plot to the next) by the secondary peak (i.e. the second-largest global peak), which has been increasing all along. So when this happens, the rate of change of the global peak suddenly changes, since it is now associated with a different peak, and thus there is a discontinuity in the slope of the peak response curve. These peak responses as a function of frequency are summarized in the plot insert at the lower-right corner of Figure 10.
Figure 9. Evolution of peak acceleration as natural frequency is increased (left to right, then top to bottom).
Figure 10. Evolution of peak acceleration as secondary peak overtakes the dominant peak. The first six points come from the preceding plots.
In principle, there are actually three possible types of behavior that can lead to discontinuities in the derivative of the peak response curve, and all can be understood in terms of the preceding logic:
- A decreasing peak is overtaken by an increasing peak (this is the case described in the preceding).
- An increasing peak is overtaken by a more rapidly increasing peak.
- A decreasing peak passes a more slowly decreasing peak so that the more slowly decreasing peak is now the dominant peak (possible in principle, but not observed in this example).
The later the peak (i.e. the larger the natural frequency of the system), the longer the system has to build up to a steady-state-like response, so that successive peak accelerations (corresponding to successively higher natural frequencies) attain higher and higher values, hence the overall upward general trend of the peak response curve. For this same reason, at high sweep frequencies successive peaks in the response versus time curve all have very similar amplitudes, so that when the natural frequency 
is changed slightly and one peak overtakes another, the difference in the rates at which the dominant and secondary peaks are increasing is extremely small and barely noticeable. Thus the peak response curve appears to be smooth at high frequencies.
10. Peak Response Curves for Octave Sweep
As described earlier, in the octave sweep case, it is extremely helpful to first make a change of independent variable by letting 
. The resulting differential equation for 
then has a constant-frequency forcing term in the 
domain (at the expense of coefficients in the equation that are at most quadratic in time). The resulting differential equation for 
, equation (8), was solved both analytically (equation (13)) and numerically, and then transformed back to the time domain.
10.1 Numerical Integration in the 
Domain
The time 
at which the sweep frequency equals the oscillator’s resonant frequency is given by equation (15). However, since the integration is being done in the 
domain, the corresponding expression for the 
value at which the system’s resonant frequency is reached becomes
 |
(16) |
Since in the 
domain the forcing function is a constant-frequency sine wave, we found experimentally that in most cases it was sufficient to integrate to a maximum 
value of 1.5 
, although occasionally it was necessary to go up to 3 or 4 times 
. In some cases it is possible for the value of 
to become less than 1, and so we also imposed a lower bound of 1.05 on the maximum value of 
. We again used fifth-order interpolation for computing derivatives in 
and again allowed the integration to go for a maximum of 100,000,000 time steps: recall that for octave sweep we used the substitution 
, so 
increases exponentially with 
, and thus the number of steps in the 
domain can become much larger than the number of time steps in the 
domain.
10.2 Numerical Optimization to Identify Peak Response
Once the differential equation (for a given set of 
, 
, 
and values) had been solved, the following procedure for finding the peak response was followed:
- Create a list of the
values generated via numerical integration. - Use equation (9) to evaluate
(in the time domain) at each 
value and then from this list select the largest response. - Use the data from steps (1) and (2) to also create an interpolating function for
as a function of 
. - Having found this initial estimate of the peak value, then use the interpolation function returned by
to do a local optimization (via the 
function) around this initial peak, using this peak as a starting point.
10.3 Peak Response Curves for Octave Sweep
Figure 11 shows the normalized peak response to various octave sweep rates. Each of the (almost 1,000) peak response values on each of these curves was computed via the process for computing peak acceleration (from the exact solution) described in Sections 7 and 8, that is, searching within ±60 cycles of the numerically generated estimate of when the peak acceleration occurs, with 1,440 phase increments per cycle. As expected, the slower the sweep rate, the lower the attenuation. In addition, the scalloped behavior in the peak response curves that was observed for the linear sweeps is also present here, although not as pronounced. This is because the octave sweep increases in frequency more rapidly than the linear sweep.
Figure 11. Normalized peak responses with 
and several values of octave sweep rate 
(octaves/minute). At low natural frequencies, 
, the peak response was computed in increments of 0.002 Hz. Left-to-right curves correspond to top to bottom in key.
Figure 12 shows an expanded view of Figure 11 corresponding to the lower frequency systems so that the scalloped behavior can be better seen. Figure 12 was obtained from Figure 11 by simply adjusting the vertical and horizontal plot ranges.
Figure 12. Expanded view of peak response curves with 
at low natural frequencies for various octave sweep rates 
(octaves/minute). Left-to-right curves correspond to top to bottom in key.
Figure 13 shows the results from Figure 11 normalized by the octave sweep rate, as suggested by Hawkes [5]. The data used in Figure 13 is the same as the data used in Figure 11, only plotted differently. As in the case with the linear sweep rate and its normalization factor, the octave sweep rate results also merge into a single curve for systems with the same critical damping ratio.
Figure 13. Normalized peak response curves for 
and various octave sweep rates plotted against 
, where 
is in Hz and 
is in octaves/minute.
Figure 14 shows comparable results to those in Figure 13 for systems with a critical damping ratio of 
.
Figure 14. Normalized peak response curves for 
and various sweep rates with 
, plotted against 
; 
is in Hz and 
is in octaves/minute.
Figures 15 and 16 show the severe attenuation that occurs when the start frequency of the sweep is close to the natural frequency. In both figures, the sweeps were started at 1 Hz. As can be ascertained, the attenuation is significant for systems with natural frequencies close to or below 1 Hz, as would be expected. Hence, the attenuation is not only a function of the natural frequency, damping and sweep rate, but also of the proximity of the start frequency of the sweep to the natural frequency. As with Figure 11, each of the peak response values on each of the curves in Figures 14–16 was computed via the process for computing peak acceleration (from the exact solution) described in Sections 7 and 8.
Figure 15. Normalized peak response curves for octave sweep with 
and 
(instead of 1/8 Hz). Left-to-right curves correspond to top to bottom in key.
Figure 16. Normalized peak response curves for octave sweep with 
and 
(instead of 1/8 Hz). Left-to-right curves correspond to top to bottom in key.
Conclusion
The derivation of closed-form solutions for the responses of single-degree-of-freedom systems subject to linear and octave swept-frequency harmonic excitation was presented. The closed-form solutions were compared to results obtained by direct numerical integration of the equations of motion with excellent agreement obtained. In addition, an in-depth discussion was presented on the numerical difficulties associated with the gamma and error functions of complex arguments that are part of the closed-form solutions, and how these difficulties were overcome by employing exact arithmetic with infinite-precision numbers, that is, rational and/or exact symbolic numbers. This included a study of precision requirements by performing computations with numerical precision exceeding what is available on today’s computers. The closed-form solutions allowed the in-depth study of several interesting phenomena including: (a) computation of the peak response instead of the peak of the envelope function; (b) scalloped behavior of the peak response with frequent discontinuities in the derivative; (c) the significant attenuation of the peak response if the sweep frequency is started at frequencies near or above the natural frequency; and (d) the fact that the swept-excitation response could exceed the steady-state harmonic response when the system is excited at its natural frequency.
Acknowledgments
We are grateful to Luke Titus of Wolfram Research for his valuable suggestions on exact numerical computation.
This work was supported by contract # FA8802-14-C-0001.
References
| [1] | F. M. Lewis, “Vibration during Acceleration through a Critical Speed,” Transactions of the American Society of Mechanical Engineers, 54(1), 1932 pp. 253–261. |
| [2] | R. L. Fearn and K. Millsaps, “Constant Acceleration of an Undamped Simple Vibrator through Resonance,” The Aeronautical Journal, 71(680), August 1967 pp. 567–569. https://doi.org/10.1017/S0001924000055007. |
| [3] | D. L. Cronin, Response of Linear, Viscous Damped Systems to Excitations Having Time-Varying Frequency, Ph.D. thesis, Dynamics Laboratory, California Institute of Technology, Pasadena, California, 1965. https://authors.library.caltech.edu/26518. |
| [4] | R. Gasch, R. Markert and H. Pfutzner, “Acceleration of Unbalanced Flexible Rotors through the Critical Speeds,” Journal of Sound and Vibration, 63(3), 1979 pp. 393–409. https://doi.org/10.1016/0022-460X(79)90682-5. |
| [5] | P. E. Hawkes, “Response of a Single-Degree-of-Freedom System to Exponential Sweep Rates,” Shock, Vibration and Associated Environments, Part II, Bulletin No. 33, February 1964 pp. 296–304. https://apps.dtic.mil/dtic/tr/fulltext/u2/432931.pdf. |
| [6] | J. A. Lollock, “The Effect of Swept Sinusoidal Excitation on the Response of a Single-Degree-of-Freedom Oscillator,” in 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2002, Denver, CO. https://doi.org/10.2514/6.2002-1230. |
| [7] | R. Markert and M. Seidler, “Analytically Based Estimation of the Maximum Amplitude during Passage through Resonance,” International Journal of Solids and Structures, 38(10–13), 2001 pp. 1975–1992. https://doi.org/10.1016/S0020-7683(00)00147-5. |
| [8] | L. Ahlfors, Complex Analysis: An Introduction to the Theory of Analytic Functions of One Complex Variable, New York: McGraw-Hill, 2000. |
| C. C. Reed and A. M. Kabe, “Peak Response of Single-Degree-of-Freedom Systems to Swept-Frequency Excitation,” The Mathematica Journal, 2018. https://doi.org/10.3888/tmj.21-1. |
|
About the Authors
Dr. Chris Reed is a Senior Engineering Specialist in the Structures Department at The Aerospace Corporation. As an applied mathematician, his work has encompassed mechanical vibrations, structural deformation, space-based sensor system performance, satellite system design optimization, flight termination system interference, fluid sloshing, electrostatic discharges, dielectric degradation on satellites and queueing systems. He has two patents and received a Wolfram Innovator award in 2017. His B.S. is from the California Institute of Technology and his M.S. and Ph.D. degrees are from Cornell University.
Dr. Alvar M. Kabe is the Principal Director of the Structural Mechanics Subdivision of The Aerospace Corporation. He has made notable contributions to the state of the art of launch vehicle and spacecraft structural dynamics. He has published numerous papers, is an Associate Fellow of the AIAA, and has received The Aerospace Corporation’s Trustees’ Distinguished Achievement Award and the Aerospace President’s Achievement Award. His B.S., M.S. and Ph.D. degrees are from UCLA.
C. Christopher Reed
Senior Engineering Specialist
Structures Department
M4-912
The Aerospace Corporation
P.O. Box 92957
Los Angeles, CA 90009-2957
cchristopher.reed@aero.org
Alvar M. Kabe
Principal Director
Structural Mechanics Subdivision
M4-899
The Aerospace Corporation
P.O. Box 92957
Los Angeles, CA 90009-2957
alvar.m.kabe@aero.org