Full report

Original page 1

NASA/TM—2016–218227 Extracting Damping Ratio From Dynamic Data and Numerical Solutions M.J. Casiano Marshall Space Flight Center, Huntsville, Alabama September 2016

Original page 2

The NASA STI Program…in Profile Since its founding, NASA has been dedicated to the advancement of aeronautics and space science. The NASA Scientific and Technical Information (STI) Program Office plays a key part in helping NASA maintain this important role. The NASA STI Program Office is operated by Langley Research Center, the lead center for NASA’s scientific and technical information. The NASA STI Program Office provides access to the NASA STI Database, the largest collection of aeronautical and space science STI in the world. The Program Office is also NASA’s institutional mechanism for disseminating the results of its research and development activities. These results are published by NASA in the NASA STI Report Series, which includes the following report types: • TECHNICAL PUBLICATION. Reports of completed research or a major significant phase of research that present the results of NASA programs and include extensive data or theoretical analysis. Includes compilations of significant scientific and technical data and information deemed to be of continuing reference value. NASA’s counterpart of peerreviewed formal professional papers but has less stringent limitations on manuscript length and extent of graphic presentations. • TECHNICAL MEMORANDUM. Scientific and technical findings that are preliminary or of specialized interest, e.g., quick release reports, working papers, and bibliographies that contain minimal annotation. Does not contain extensive analysis. • CONTRACTOR REPORT. Scientific and technical findings by NASA-sponsored contractors and grantees. • CONFERENCE PUBLICATION. Collected papers from scientific and technical conferences, symposia, seminars, or other meetings sponsored or cosponsored by NASA. • SPECIAL PUBLICATION. Scientific, technical, or historical information from NASA programs, projects, and mission, often concerned with subjects having substantial public interest. • TECHNICAL TRANSLATION. English-language translations of foreign scientific and technical material pertinent to NASA’s mission. Specialized services that complement the STI Program Office’s diverse offerings include creating custom thesauri, building customized databases, organizing and publishing research results…even providing videos. For more information about the NASA STI Program Office, see the following: • Access the NASA STI program home page at http://www.sti.nasa.gov • E-mail your question via the Internet to help@sti.nasa.gov • Phone the NASA STI Help Desk at 757 –864–9658 • Write to: NASA STI Information Desk Mail Stop 148 NASA Langley Research Center Hampton, VA 23681–2199, USA

Original page 3

NASA/TM—2016–218227 Extracting Damping Ratio From Dynamic Data and Numerical Solutions M.J. Casiano Marshall Space Flight Center, Huntsville, Alabama National Aeronautics and Space Administration Marshall Space Flight Center • Huntsville, Alabama 35812 September 2016 i

Original page 4

Available from: NASA STI Information Desk Mail Stop 148 NASA Langley Research Center Hampton, VA 23681–2199, USA 757–864–9658 This report is also available in electronic form at http://www.sti.nasa.gov ii

Original page 5

TABLE OF CONTENTS 1. INTRODUCTION ............................................................................................................. 1 1.1 Half-Power Method ...................................................................................................... 3 1.2 Half-Quadratic Gain Method ....................................................................................... 4 1.3 Logarithmic Decrement Method .................................................................................. 4 1.4 Autocorrelation and Power Spectral Density Method ................................................... 8 1.5 Frequency Response Function Method ......................................................................... 8 1.6 Random Decrement Method ........................................................................................ 9 1.7 Numerical Models ........................................................................................................ 10 APPENDIX A — QUALITY FACTOR, HALF-POWER METHOD, AND HALF-QUADRATIC GAIN METHOD .......................................... 14 A.1 Quality Factor .................................................................................................................. 14 A.2 Half-Power Method .......................................................................................................... 15 A.3 Half-Quadratic Gain Method ........................................................................................... 20 APPENDIX B — FITTING A DECAYING EXPONENTIAL TO OSCILLATION MAXIMA .................................................................................................... 22 APPENDIX C — OTHER RELATIONSHIPS ....................................................................... 28 REFERENCES ....................................................................................................................... 29 iii

Original page 6

LIST OF FIGURES 1. Power and quadratic gain function plots using an undamped natural frequency of fn = 1,000 Hz and three different damping ratios: (a) ζ = 0.05, (b) ζ = 0.10, and (c) ζ = 0.20 ............................................................................................... 19 2. Plot showing the expected amplitude level of the normalized quadratic gain at the half-power point ................................................................................................ 19 3. Comparison of amplitude envelope to exponential through local maxima: (a) Free decay response, envelope, and exponential through local maxima and (b) zoomed view ................................................................................................... 23 LIST OF TABLES 1. Relationship for damping ratio ................................................................................... 2 iv

Original page 7

LIST OF ACRONYMS COMSOL COMSOL Multiphysics FFT fast Fourier transform PSD power spectral density SDOF single degree of freedom v

Original page 8

NOMENCLATURE A(t) amplitude envelope a real part of traditional complex eigenvalue, variables in trigonometric identities b imaginary part of traditional complex eigenvalue, variables in trigonometric identities c speed of sound E energy F force F(t) time-dependent force f frequency fd damped natural frequency fi imaginary part of the complex eigenfrequency fl lower frequency at a half-amplitude level fl,H2 lower frequency at the half-quadratic gain level fl,P lower half-power point (lower frequency at the half-power level) fn undamped natural frequency fp peak center frequency of a response fr real part of the complex eigenfrequency fu upper frequency at a half-amplitude level fu,H2 upper frequency at the half-quadratic gain level fu,P upper half-power point (upper frequency at the half-power level) G input power spectral density amplitude Gyy one-sided output power spectral density amplitude H{ } Hilbert transform H(f) frequency respone function k wave number, representative stiffness vi

Original page 9

NOMENCLATURE (Continued) ks representative stiffness m system mass n number of cycles, integer P peak pressure amplitude at a local maxima, pressure amplitude, complex pressure amplitude, power p pressure Q quality factor Ryy autocorrelation function t time X amplitude factor x spatial coordinate, generic dependent variable α spatial absorption coefficient b temporal absorption coefficient d logarithmic decrement ζ damping ratio l complex eigenvalue li imaginary part of the COMSOL complex eigenvalue lr real part of the COMSOL complex eigenvalue 2 σ output spectrum variance y t time duration, autocorrelation time delay td period of a damped oscillation φ phase angle w angular frequency wd angular damped natural frequency wn angular undamped natural frequency wp angular peak center frequency vii

Original page 10

viii

Original page 11

TECHNICAL MEMORANDUM EXTRACTING DAMPING RATIO FROM DYNAMIC DATA AND NUMERICAL SOLUTIONS 1. INTRODUCTION There are many ways to extract damping parameters from data or models. This Technical Memorandum provides a quick reference for some of the more common approaches used in dynamics analysis. Described are six methods of extracting damping from data: the half-power method, logarithmic decrement (decay rate) method, an autocorrelation/power spectral density fitting method, a frequency response fitting method, a random decrement fitting method, and a newly developed half-quadratic gain method. Additionally, state-space models and finite element method modeling tools, such as COMSOL Multiphysics (COMSOL), provide a theoretical damping via complex frequency. Each method has its advantages which are briefly noted. There are also likely many other advanced techniques in extracting damping within the operational modal analysis discipline, where an input excitation is unknown; however, these approaches discussed here are objective, direct, and can be implemented in a consistent manner. Damping ratio, ζ, is characterized since it is commonly used in many disciplines. Though, it should be noted that all linear damping parameters can be related to a single parameter such as damping ratio. There are many ways to describe damping, e.g., damping ratio, quality factor, spatial absorption coefficient, temporal damping coefficient, complex frequency, and many others. By characterizing one parameter, the goal is to have a consistent way to compare damping from data and models. There are many ways to show or derive the damping ratio relationships described in table 1, but the common assumption is that the system response behaves as a single-degree-of-freedom (SDOF) time harmonic oscillator such as described in equation (1): 2 x = F(t) / m , (1) x!!+ 2ζ ω x! + ω n n where x is the dependent variable, wn is the angular undamped natural frequency, F (t) is the timedependent force, and m is a representation of system mass. Many systems can be represented in this form including damped acoustic systems.1 (p. 8), 2 (pp. 22–26), 3 (p. 20) Several other damping parameters are also related to damping ratio so that each method can be consistently compared, whether by extracting damping from test data or from a theoretical model. Table 1 provides a quick reference. 1

Original page 12

Table 1. Relationship for damping ratio. Data Extraction Method Damping Ratio Relationship Notes Half-power method f − f u l ζ ≈ for ζ < 0.05 2 f p Half-quadratic gain method ⎛ ⎛ Convenient estimate from spectrum. Approximate for low damping ratio Simple estimate from spectrum. −1 2 4⎞ Exact for force-excited SDOF 1 f f ⎞ ⎛ f − f ⎞ − ζ = − ⎜ 4 + 4 u l ⎜ 2 ⎜ ⎝ f ⎝ p General logarithmic decrement ln(A A ) method t t 1 2 ζ = 2 (2π f τ) + ln(A d t 1 Power spectral density / Fitting method—see text autocorrelation method Frequency response method Fitting method—see text Random decrement method Fitting method—see text Other Damping Ratio − u l ⎟ response ⎟ ⎜ ⎟ ⎠ ⎝ f ⎠ ⎟ p ⎠ For a linear exponential decay of an amplitude envelope. Exact 2 for unforced SDOF decay A ) t 2 Used with PSD and autocorrelation and allows for function fitting Used on complex FFT and allows for function fitting Estimates a signal proportional to autocorrelation and allows for function fitting Relationships Relationship Notes Complex frequency f i ζ = 2 2 f + f r i Logarithmic decrement δ , where δ = ln(P P ) by selecting successive peaks. ζ = 2 2 (2π ) + δ Temporal absorption coefficient β ζ = 2π f n Spatial absorption coefficient α ⋅c ζ = 2π f n Quality factor 1 ζ = 2Q 2 Output from state-space models and finite element method models Traditional estimate. Evaluated 1 t t n 1 n Exact for unforced SDOF decay Common estimate of physical models—appendix C Classic characterization in acoustics and other branches of physics —appendix C Classic figure of merit—appendix A

Original page 13

Separately, it is noted that a comparison of damping ratio between a state-space model, or COMSOL output, and test data would also serve as a quantitative metric that may be used as a level of validation of linearly stable systems. In evaluating combustion dynamics of a combustor, for example, a damping parameter can be obtained from stationary and linearly stable test data, at the onset of an instability where a clear linear growth is observed, or in the best case, during an unforced free response decay. A linearized model cannot be compared to a nonlinear instability using these metrics as the physics must follow the form of equation (1). Real systems do not behave linearly in and near the ‘unstable’ regime where behavior includes limit cycles, nonlinearly coupled harmonics and modulations, and complex chaotic behavior. In the combustion dynamics example, equation (1) does not hold in this ‘unstable’ regime, i.e., damping parameters from data can only be compared to a linearized model while in a stable regime. 1.1 Half-Power Method The half-power method assumes an SDOF time harmonic oscillator and is best applied to stationary data. The half-power method is commonly used in obtaining estimates of modal damping for structural vibrations.4 (pp. 110–120) Quality factor, Q, is given exactly as equation (2) and is developed in appendix A: Q = 2ζ 1 . (2) Quality factor can be estimated with equation (3),5 (p.16) also developed in appendix A: f p Q ≈ for ζ < 0.05 . (3) f − f u l The half-power method consists of selecting measured frequencies from a spectrum. The measured peak center frequency of a response, fp, is selected assuming it is the frequency with maximum power. This is an approximation because the maximum power for an SDOF occurs at the undamped natural frequency. Frequencies are also selected above the peak, fu , and below the peak, fl , assuming they occur at the half-power level. With the assumption that the square of the SDOF frequency response function (which is proportional to a power spectral density (PSD)) is proportional to power, a PSD can be used to obtain the frequencies at the half-power level where the amplitude is one-half the maximum peak and a fast Fourier transform (FFT) can be used to obtain the frequencies at the half-power level where the amplitude is 1/ 2 the maximum peak. In reality, this is true only as damping approaches zero. The quality factor approximation is only accurate at low damping values, i.e., ζ < 0.05. Combining equations (2) and (3) gives a convenient way of estimating damping ratio shown in equation (4): f − f u l ζ ≈ for ζ < 0.05 . (4) 2 f p 3

Original page 14

The approximations and associated error are discussed further in appendix A, which can be helpful in understanding the acceptability of an estimated damping ratio or quality factor. 1.2 Half-Quadratic Gain Method The half-quadratic gain method is developed in appendix A. It is an exact formula for a forced, SDOF system and is given as equation (5): ⎛ ⎛ 1 f − ζ = − ⎜ 4 + 4 u ⎜ 2 ⎜ ⎝ f ⎝ p −1 2 4⎞ f ⎞ ⎛ f − f ⎞ l − u l ⎟ . ⎟ ⎜ ⎟ (5) ⎠ ⎝ f ⎠ ⎟ p ⎠ Instead of approximating the power curve with the square of the frequency response function (or PSD) as in the half-power method, the half-quadratic gain method uses the square of the frequency response function directly in estimating damping. This provides an exact method of extracting damping ratio from an FFT or PSD since for an SDOF system, an FFT spectrum is proportional to the transfer function gain and a PSD is proportional to the square of the transfer function gain. Unlike in the half-power method, the measured peak center frequency, fp, is equal to the frequency of peak amplitude response for an SDOF and does not necessarily occur at the frequency with maximum power. A PSD can be used to obtain the frequencies at the half-quadratic gain level where the amplitudes are one-half the maximum peak and an FFT can be used to obtain the frequencies at the half-quadratic gain level where the amplitudes are 1 2 the maximum peak. The frequencies at the half-quadratic gain level are exactly the frequency selected above the peak, fu, and below the peak, fl , of the measured spectrum. 1.3 Logarithmic Decrement Method The logarithmic decrement method assumes an SDOF time harmonic oscillator and is best applied to underdamped exponentially decaying data. It is commonly used in obtaining modal damping for structural vibrations.4 (pp. 87–98), 6 (pp. 52,53) However, given an impulse response from an aperiodic transient, such as a combustion stability bomb in the combustion dynamics example, a similar approach can be used to estimate the decay rate at a given frequency, assuming there is an acoustically exponential decay of the amplitude envelope. First, a general approach will be described followed by the classic logarithmic decrement. For an underdamped system, 0 < ζ < 1, the particular solution (no force) of equation (1) is uation (6):7 (p. 245) the underdamped response and given by eq p(0) −ζ ω t −1 p(t) = ⋅X ⋅e 2 1−ζ n cos(ω t − tan (φ)) . (6) d Pressure, p, is the dependent variable, wd, is the angular damped natural frequency, and p(0) is the initial pressure amplitude condition. The factor X and the phase angle φ are discussed in appendix B. 4

Original page 15

The angular damped natural frequency can be written as equation (7) for an SDOF. This is the frequency that is observed during an SDOF free decay as indicated by the cosine argument in equation (6): ω = ω d n 2 1−ζ . (7) equation (6) is given as equation (8):7 (p. 249) The decaying amplitude envelope of p(0)⋅X −ζ ω t A(t) = 1−ζ e n , (8) 2 where A(t) is the amplitude envelope. Analytically, using the Hilbert transform of the pressure, H{p(t)}, it can also be shown that the amplitude envelope follows the form as expressed in equation (9):8 { )} = e n . (9) A(t) = p(t) + iH p(t p(0)⋅X −ζ ω t 2 1−ζ It is common to apply the Hilbert transform numerically to a data set to remove the oscillatory components and reduce a time history to an amplitude envelope. The decaying exponent can be defined using a parameter called the temporal absorption coefficient, also referred to as the decay rate, shown in equation (10): ζω . (10) β ≡ n A ratio of two amplitudes on the decaying exponential curve can be written using equation (8) and is shown as equation (11) where t1 and t2 are subsequent times at any given amplitude envelope value: A(t ) −β(t −t ) 2 = e A(t ) 1 2 1 . (11) Equations (12) and (13) result when using A t( ) = A , and A t( ) = A , where A and A are measured amplitude envelope values at t1 and t2: A = A e t 2 and 1 β = ln τ 1 t 2 t t t 1 2 1 2 −βτ t (12) 1 ⎛ A ⎞ t 1 ⎜ ⎟ . (13) ⎜ A ⎟ ⎝ t ⎠ 2 5

Original page 16

Also, t = t2 – t1 where the parameter t is the measured duration between the time of the measured amplitude envelope values. The response curve, equation (6), is tangent to the amplitude envelope curve, equation (8), and the tangents do not occur at the oscillation peaks, i.e., the local oscillation maxima are not values that occur on the amplitude envelope curve.6 (pp. 52,53), 7 (p. 249) However, it is important to note that while values from the exact amplitude envelope are needed to produce an exponential fit to the entire equation (8), the damping ratio can still be determined exactly using the local oscillation maxima to produce an exponential fit since the same decay rate is retained. Only the amplitude of equation (8) would be incorrect with a fit to the local oscillation maxima; however, an amplitude correction can be applied. This is discussed further in appendix B. Equations (2) and (10) can be used to express the temporal absorption coefficient in terms of damping ratio, shown in equation (14): ζ = 2π f β . (14) n However, the observed oscillations occur at the damped natural frequency for a decaying transient, so equation (7) can be used to put the damping ratio in terms of measurable quantities: ζ = (2π f ) + β β . (15) 2 2 d Simplifying further in terms of measurable parameters using equation (13) gives equation (16): (A A ) ln ζ = 2 t t 1 2 . (16) 2 (2π f τ) + ln(A A ) d t t 1 2 This is exact for a free decay and can be applied over any duration given values for the amplitude envelope. For the linear period at the onset of a growing signal the formulation can only be used as an approximation of a growth rate since the previous derivations are based on an ‘unforced’ decay or free damped vibration. In theory, while an actual envelope needs to be estimated for the general logarithmic decrement approach described in equation (16), as noted previously, the local oscillation maxima values can also be used to obtain the damping ratio. Using successive oscillation peaks to calculate damping is an approach known classically as the logarithmic decrement method. It is shown in appendix B that this procedure also provides an exact value for damping ratio and that local oscillation maxima values follow the same decay rate as the amplitude envelope. The logarithmic decrement, d, is given as equation (17) where the sinusoidal peak amplitudes, P and P , are chosen between successive peaks:6 (pp. 52,53) 6 t t 1 2

Original page 17

(P P ) . (17) δ ≡ ln t t lized to n cycles as equation (18):6 (pp. 55,56) The logarithmic decrement can be genera 1 1 2 δ = ln(P P ) , (18) n t t 1 n where Pt represents the sinusoidal peak amplitude value after n cycles. Combining equations (6) n and (17), where P(t1) = Pt and P(t2) = Pt results in equation (19) and is further discussed 1 2 in appendix B: ⎛ −ζ ω t ⎞ e n 1 δ = ln⎜ ⎟ = −ζ ω t ⎝ e n 2⎠ ζ ω (t − t ) . (19) n 2 1 This can then be simplified in terms of damping ratio shown in equation (20): ζ = 2π f τ δ , (20) n d where td is the period of a damped oscillation since this is the time duration between peaks of free damped vibrations. The period can simply be calculated from equation (7) and is written as equation (21): τ = d f n 1 . (21) 2 1−ζ Combining equations (20) and (21) and solving for damping ratio gives equation (22): ζ = (2π ) + δ δ . (22) 2 2 Simplifying further in terms of measurable parameters, using the multiple-cycle form of logarithmic decrement equation (18), gives equation (23): ln(P P ) t t ζ = 2 1 n . 2 (23) (2πn) + ln(P P ) t t 1 n 7

Original page 18

Since the exponential decay rate through the local maxima is the same as the decay rate of the amplitude envelope, it can be observed that equation (23) is equivalent to the general logarithmic decrement equation (16) applied over the duration of n cycles. 1.4 Autocorrelation and Power Spectral Density Method This method assumes an SDOF time harmonic oscillator that has a constant power spectrum input (white noise) and is applicable for stationary data. For the combustion dynamics example, constant power spectrum input is approximate for combustion noise across a narrow bandwidth of interest. In most applications, the excitation is sufficiently broad so as not to introduce significant errors. The PSD and autocorrelation functions are known functions for an SDOF mechanical system with a stationary process.9 (pp. 124–126) These are given as equations (24) and (25), respectively: G (τ ) = yy ⎡1− f f 2 G / k s 2 (24) 2⎤ 2 ) + ⎡2ζ f f ⎤ ⎢ ( n ⎥ ⎣ n ⎦ ⎣ and −2π f ζ τ ⎛ πG f e n n ⎦ ⎞ ζ R (τ ) = ⎜ cos 2(π f τ) + sin 2(π f τ )⎟ , (25) yy 2 4ζ k ⎜ s ⎝ d 2 d 1−ζ ⎟ ⎠ where Gyy is the one-sided output PSD, G is the input PSD amplitude (constant in this case with white noise), ks is the representative stiffness, f is the frequency, fn is the undamped natural frequency, fd is the damped natural frequency, Ryy is the autocorrelation function, and t is the autocorrelation time delay. The measured data can be reduced using the PSD and autocorrelation functions and then fit to these forms to extract a damping ratio. To aid in fitting, equation (26) gives the variance of the output spectrum: 2 π G f n R (0) = σ = . (26) yy y 2 4ζ k s This is discussed further in references 9 (pp. 124 – 126) and 10. 1.5 Frequency Response Function Method This method is similar to section 1.4 but directly fits data to a second-order system response. The frequency response function gain, maximum gain, and phase for a second-order system with harmonic excitation are known functions and given as equations (27) through (29), respectively.9 (pp. 18–21) Normally, these would be fit using linear cross spectrum (transfer function) given an input and output measurement. This model assumes an SDOF time harmonic oscillator and is best applied to stationary data; however, in general, a frequency sweep input excitation can be 8

Original page 19

performed at sweep rates slow enough to capture the full response. Assuming the input is harmonic over the spectrum, or in the case of white noise with superposition of harmonic components at various frequencies, a damping ratio approximation can be obtained. H( f ) = ⎡1− f f F k , (27) 2 2⎤ 2 ) + ⎡2ζ f f ⎤ ⎢ ( n ⎥ ⎣ n ⎦ ⎣ H = max 2ζ 1−ζ and ⎛ −1 φ( f ) = tan ⎜ ⎜ 1−( f f ) ⎝ ⎦ F k , (28) 2 ⎞ 2ζ f f n ⎟ . (29) 2 ⎟ n ⎠ It is noted that some software tools may plot the fast Fourier transform (FFT) phase using the twoargument arctangent which covers the entire circular range and equation (29) would be written as −1 2 φ( f ) = tan2 (2ζ f f , 1−( f f ) ) . Measured data can be reduced using the complex FFT and n n fitting the data to the FFT phase and gain functions to extract a damping ratio. Its advantage is the ability to optimize an estimate by fitting both gain and phase; however, it is also possible to fit the form using just the gain or just the phase. 1.6 Random Decrement Method The concept of random decrement is based on averaging enough samples of the same data set to remove the random part and retain the deterministic part.11 (p. 20) The deterministic part for an SDOF response will be exactly proportional to the autocorrelation function under stationary and Gaussian white noise conditions.12 Under these conditions, the autocorrelation function is also exactly proportional to the free decay response of the SDOF system. So the recovered deterministic part, referred to as the Randomdec signature, can be used to interrogate the free decay of an SDOF system. An amplitude value is set at a level where the data fluctuations regularly exceed, filtered for a specific mode. The filter signal captures the response of just the single mode. For each instant of time the data set reaches this set amplitude level, a new time history is created by shifting the data back to the initial time that the data first reached the set amplitude level. All of the time history data sets are then averaged together. This signal is the Randomdec signature and is simply a trace formed by a waveform averaging of a number of specially selected segments from an observed time history.12 9

Original page 20

Damping ratio can be found by fitting the Randomdec signature using a Gaussian white noise SDOF relationship such as equation (25) or applying a logarithmic decrement approach as in equations (16) or (23). A simple approach is to fit the amplitude envelope curve of the Randomdec data to the exponential form of equation (8). 1.7 Numerical Models Eigenanalysis, or modal analysis, can be used to extract dynamic properties of systems including damping ratio. This can be done traditionally using state-space modeling or for more complex systems using the finite element method. First, the traditional state-space representative approach is described and then an approach using COMSOL Multiphysics is described. COMSOL Multiphysics is a finite element analysis, solver, and simulation software package for various physics and engineering applications, especially coupled phenomena, or multiphysics. Care must be made when examining the solution in COMSOL. In the COMSOL eigen-solver, COMSOL does not use the traditional definition of an eigenvalue. 1.7.1 Traditional Approach system is assumed to behave as the model in equation (1)1 (p. 8), A damped acoustic 2 (pp. 22–26), 3 (p.20) and a pressure solution is described by equation (30), where P is the complex pressure amplitude:2 (pp. 22–26), 5 (p. 212) Pe . (30) p = iω t The traditional state-space system defines the eigenvalue as λ =ˆ iωˆ and so, the frequency, in terms of the eigenvalue, is expressed as f = −λˆi 2π. The complex eigenvalue is defined as equation (31): λ λ = a + bi , (31) trad where a and b are the real and imaginary parts of the traditional complex eigenvalue, ltrad. For a second-order system, the eigenvalues (or system poles) are found such that the real and imaginary parts are a = – ζ wn and b = ±wd. The angular natural frequency and damping ratio can then be written as equations (32) and (33), respectively: 2 2 ω = a + b (32) n and −a ζ = 2 2 a + b . (33) Equations (32) and (33) can also be observed clearly by examining an underdamped second-order system in the complex plane. 10

Original page 21

The angular damped natural frequency is given as equation (34): ω = ω d n 2 1−ζ = b (34) and the angular damped resonant frequency, or angular peak center frequency, wp, is given as9 (pp. 18–21) equation (35): 2 2 2 ω = ω 1− 2ζ = b − a . (35) p n ition of the eigenvalue, λ =ˆ i2π fˆ, and the definition from equa- Using the traditional defin tion (31), the complex eigenfrequencies can easily be determined as equation (36): b a f = f + f i = − i , (36) trad r i 2π 2π where fr and fi are the real and imaginary part of the complex eigenfrequency. So equation (36) gives a = –2 p fi and b = 2 p fr . These can be substituted into equations (32)–(35) which gives the natural frequency, damping ratio, damped natural frequency, and peak frequency in terms of the complex eigenfrequency shown in equations (37)–(40), respectively; 2 2 f = f + f , (37) n r i ζ = 2 2 f + f r i f = f , (39) d and f i , (38) r 2 2 f = f − f ,. (40) p r i 1.7.2 COMSOL Approach Opposed to the traditional definition, COMSOL defines the eigenvalues as λ = −ˆiωˆ for a pressure solution of the form in equation (30) in the acoustic eigensolver and so, the frequency, in terms of the eigenvalue, is expressed as f = λˆi 2π . This is a subtle difference from the tradi- λ tional definition, but is important in the interpretation of the COMSOL results. 11

Original page 22

The COMSOL complex eigenvalue is described in equation (41): λ COMSOL Because the COMSOL defined notation is λ COMSOL trad = λ + λ i . (41) r i = −λ , equation (31) gives the COMSOL eignevalues real and imaginary parts as lr= –a and li= –b. Using these relationships in equations (32) and (33) gives the natural frequency and damping ratio in the COMSOL framework as equations (42) and (43), respectively: 2 2 ω = λ + λ (42) n and ζ = λ + λ r i r i λ r . (43) 2 2 Substituting these relationships into equations (34) and (35) as well gives the angular damped natural frequency and angular damped resonant frequency as equations (44) and (45): 2 ω = ω 1−ζ = λ (44) d n and 2 2 2 ω = ω 1−2ζ p n i = λ − λ . (45) i r Using the COMSOL defined eigenvalue, λ = −ˆi2π f, complex eigenfrequenies can easily be determined: f i = − i + r i . (46) f = f + COMSOL r ˆ and equation (41), the COMSOL λ λ i 2π 2π So equation (46) gives li = –2 p fr and lr = 2 p fi . These can be substituted into equations (42)–(45). This gives the natural frequency, damping ratio, damped natural frequency, and peak center frequency, in terms of the complex eigenfrequency shown in equations (47)–(50), respectively: f = f n r i ζ = 2 2 f + f r 12 2 2 + f , (47) f i , (48) i

Original page 23

f = f , (49) d and f = f p r i r 2 2 − f . (50) While the actual eigenvalues in COMSOL are negative of the traditional eigenvalues, the key parameters, undamped natural frequency, damping ratio, damped natural frequency, and peak frequency, are all estimated in the same manner in COMSOL as in the traditional approach. In other words, equations using the traditional approach, equations (37)–(40), are identical to equations using the COMSOL approach, equations (47)–(50). The damping ratio can be calculated from eigenvalues using equation (33) for traditional state-space models and from equation (43) for COMSOL. The damping ratio can also be calculated from complex frequency; equation (38) for a state-space model is identical to equation (48) for COMSOL. 13

Original page 24

APPENDIX A — QUALITY FACTOR, HALF-POWER METHOD, AND HALF-QUADRATIC GAIN METHOD A.1 Quality Factor Quality factor is defined for a harmonic excitation at the undamped natural frequency , divided by the total energy lost, E , in one cycle:13 (pp. 3–13) as the total energy stored, Estored E Q ≡ 2π lost stored . (51) E lost As described in reference 13, the total energy stored for the SDOF mechanical system can be written as equation (52) using the response functions described in section 1.5: 1 2 E = kx stored max 2 k 2 = H( f ) , (52) 2 where xmax is the maximum value of the dependent variable, e.g., position or pressure, and the energy dissipated over one cycle is equivalent to the work done by the applied force over one cycle: dx = F(t)x!(t) dt . (53) E = ∫ F(t) lost 1 cycle 2π /ω ∫ 0 For an applied harmonic force of F sin(w t), the SDOF solution can be found and the terms in equation (53) can then be written as equation (54): 2π /ω E = ∫ F sin(ωt) ( H( f ) ω cos(ωt − φ ( f )))dt . (54) lost 0 Using equation (29), this is simplified to equation (55): ⎛ f 2 E = 2π ζ k ⎜ lost ⎝ f ⎠ 14 ⎞ ⎟ H( f ) . (55) n

Original page 25

Substituting equation (52) and equation (55) into equation (51) gives equation (56): Q = 2ζ ( f f ) 1 . (56) n With maximum power at the undamped natural frequency, the quality factor is given by equation (57): Q = 2ζ 1 . (57) A.2 Half-Power Method The time-averaged power over one period, given by equation (58), is similar to the energy dissipated over one cycle described by equation (53): 2π /ω ω P = ∫ 2π 0 F(t)x!(t)dt . (58) The average power over one cycle is then given as equation (59): ⎛ P = 2π fζ k ⎜ ⎝ f ⎠ f ⎞ 2 ⎟ H( f ) . (59) n Immediately, it is seen that power is not a scalar multiple of gain squared since there is also a functional dependence on frequency and damping ratio. Expanding power by incorporating equation (27) gives equation (60): 2 2π fζ ( f f )F P = n . 2 ⎞ (60) ⎛ ⎡ 2⎤ 2 k 1−( f f ) ⎜ ⎢⎣ n ⎝ + ⎡2ζ f f ⎤ ⎥⎦ ⎣ n ⎦ ⎟ ⎠ It can be shown by differentiating equation (60) that the maxima occurs at equation (61). Therefore, the maximum power—no matter the damping ratio—occurs at the undamped natural frequency5 (pp. 14,15) f = f (61) n . 15

Original page 26

The maximum power is found as equation (62) by substituting equation (61) into equation (60): π f F P = max 2 n . (62) 2kζ Equations (60) and (62) can be used to solve for frequencies at half the maximum power. With P = Pmax/2, the upper and lower frequencies are found by solving for damping ratio and given as equations (63) and (64), respectively: 2 f =(ζ + ζ +1) f (63) u,P and n 2 f = −(ζ + ζ +1) f , (64) l,P n where fu,P and fl,P are the upper and lower frequencies at the half-power level on the SDOF power curve. These are also referred to as the upper and lower half-power point, respectively. Subtracting equation (64) from equation (63) gives equation (65): f − f = 2ζ f . (65) u,P l,P n Equation (66) is found by rearranging equation (65): f u,P ζ = 2 f − f l,P . (66) n Therefore, it is verified that, when using the power curve, the half-power method is exact. Given the undamped natural frequency and the power response curve, a damping ratio can then be estimated exactly. However, the undamped natural frequency and the power response are not usually known, especially when extracting information from measured data. In the half-power method, an approximation is made that the square of the frequency response curve is the same functional form as the power curve. With this approximation, the frequency at the peak amplitude of a spectrum is assumed to be equivalent to the undamped natural frequency. Additionally, the frequency at the half-power level on the power response curve is assumed to occur at half the quadratic gain level on the frequency response function curve. With these approximations, the half power method can be written approximately as f − f u l ζ ≈ for ζ < 0.05 , (67) 2 f p 16

Original page 27

where the frequencies can be extracted from the square of the frequency response function as described in section 1.1. Quality factor can also be approximated from equations (57) and (67) and written as equation (68): f p for ζ < 0.05 . (68) Q ≈ f − f u l A.2.1 Half-Power Method Error The error associated with the half-power method assumptions are eliminated as damping approaches zero. This error can be visualized and estimated by comparing the power curve directly to the square of the frequency response curve. First, normalizing the power curve by the maximum power, using equations (60) and (62) gives equation (69): P P = = norm P 2 (2ζ f f ) n . (69) 2 2 2 max 1− f f + 2ζ f f ( ( n) ) ( n) The frequency response function is previously given by equation (27) and rewritten here as equation (70). It can be shown by differentiating equation (70) that the maxima occurs at equation (71), which was also mentioned as equation (35). The maximum value was previously given as equation (28) and is rewritten here as equation (72): H( f ) = ⎡1− f f ⎢ ( n ⎣ F / k , (70) 2 2⎤ 2 ) + ⎡2ζ f f ⎤ ⎥⎦ ⎣ n ⎦ f = 1− 2 ζ f ≡ f , (71) and n p F / k H( f ) = . (72) p max 2 2ζ 1−ζ 17

Original page 28

Next, normalizing the frequency response curve by the maximum power gives equation (73): H( f ) H( f ) = = 2ζ norm H( f ) p max 2 1−ζ . 2 (73) 2 2 (1−( f f ) ) +(2ζ f f ) n n Then, squaring equation (73) gives equation (74), the normalized form of the quadratic gain: 2 H( f ) = norm 2 4 4(ζ −ζ ) . 2 (74) 2 2 (1−( f f ) ) +(2ζ f f ) n n Comparing equation (69) to equation (74) directly shows that there is a notable difference. The error and fractional error can be given as equations (75) and (76), respectively; and represents the error associated with applying the half-power method to a PSD: 2 ε = P − H( f ) = norm norm (1−( f f ) ) +(2ζ f f ) and 2 P − H( f ) norm norm ε = frac P norm 2 2 2 (2ζ f f ) + 4ζ (ζ −1) n 2 (75) 2 2 n n 2 ⎛ f ⎞ n = 1−(1−ζ )⎜ ⎟ . (76) ⎝ f ⎠ Figure 1 shows three comparative plots of the normalized power curve, equation (69), and the normalized square of the frequency response function (normalized quadratic gain), equation (74), for undamped natural frequency of fn = 1,000 Hz and damping ratio of ζ = 0.05, ζ = 0.10, and ζ = 0.20. The intersection of the curve with the dashed line at an amplitude level of 0.5 represents the half-power point or the half-quadratic gain level. The separation of these points represents a margin of error. Another way to visualize the error is to examine what is the amplitude of the square of the frequency response function (quadratic gain) when the frequency is at the half-power level based on the power curve; i.e., using equations (63) and (64) in equation (74). This frequency always occurs at half the power no matter what the damping. The error in using the normalized square of the frequency response function can be visualized by observing the deviation of the value from one-half. Figure 2 shows a representation of error margin associated with using the upper and lower frequency half-power point. It is seen that as the damping value approaches zero, it becomes acceptable to use the frequency response function at the half-quadratic gain level in the half-power 18

Original page 29

method. But for values where ζ > 0.05, the amplitude level on the quadratic gain curve at the half-power point occur at values that diverge significantly from one-half. This is clearly observed in figure 2. Since the gain-squared response curve shifts left from the power curve, as observed in figure 1, the amplitude level at the upper frequency decreases and the amplitude level at the lower frequency increases from the one-half-power amplitude level at the half-power point. 1.5 1.5 Power =0.05 Power 1.5 =0.10 Power =0.20 Gain-Squared =1,000 Hz Gain-Squared =1,000 Hz Gain-Squared =1,000 Hz 1 1 0.5 0.5 1 0.5 Normalized Amplitude Normalized Amplitude Normalized Amplitude 0 0 0 700 800 900 1,000 1,100 1,200 1,300 700 800 900 1,000 1,100 1,200 1,300 700 800 900 1,000 1,100 1,200 1,300 (a) Frequency (Hz) (b) Frequency (Hz) Frequency (Hz) (c) Figure 1. Power and quadratic gain function plots using an undamped natural frequency of fn = 1,000 Hz and three different damping ratios: (a) ζ = 0.05, (b) ζ = 0.10, and (c) ζ = 0.20. 1 Upper Frequency F1_1623 Lower Frequency 0.8 0.6 0.4 0.2 Normalized Gain-Squared Amplitude Level 0 0 0.1 0.2 0.3 0.4 0.5 ζ (Damping Ratio) Figure 2. Plot showing the expected amplitude level of the normalized quadratic gain at the half-power point. F2_1623 19

Original page 30

A.3 Half-Quadratic Gain Method The half-power method is a convenient method to estimate damping for low damping ratio. However, an exact method can be developed so the spectra can be used directly to estimate damping. Instead of assuming that the square of the frequency response function is proportional to the power, the square of the frequency response function can be used directly to derive a formula to estimate damping. The half-quadratic gain method considers that the system behaves as an SDOF described by equation (1). Since a PSD is directly proportional to the square of the frequency response function, it can be used at a specified level to estimate damping exactly. The formulation is based on using the frequencies that are associated with half the peak amplitude of the square of the frequency response function. Solving equation (77) using equations (70) and (72) gives the frequencies at half the gain squared value: 2 H( f ) p max 2 2 The upper and lower frequencies, f 2 and f 2 u,H l,H = H( f ) . (77) , at the half-gain-squared (half-quadratic gain) level for an SDOF are found to be equations (78) and (79) in terms of the undamped natural frequency: 2 2 f 2 = 1−2ζ + 2ζ 1−ζ f (78) u,H and 2 n 2 f 2 = 1−2ζ − 2ζ 1−ζ f . (79) l,H n The frequency at the peak amplitude of an SDOF is known in terms of the undamped natural frequency from equation (71) and can be substituted into equations (78) and (79). This results in equations (80) and (81): f 2 = 1+ u,H and 2ζ 1−ζ f 2 = 1− l,H 20 2 2ζ 1−ζ f (80) 2 p 1−2ζ 2 f . (81) 2 p 1−2ζ

Original page 31

Although the expressions are much more complex, a similar approach as that used in the half-power method to obtain equation (65), is used to obtain a half-quadratic gain formula. Subtracting equation (81) from equation (80) gives equation (82): ⎛ 2 2 ⎞ 2ζ 1−ζ ζ 1−ζ f − f = ⎜ 1+ 2 2 2 − 1− ⎟ f . (82) 2 2 p u,H l,H ⎜ 1−2ζ 1−2ζ ⎝ ⎟ ⎠ Solving this for damping ratio gives equation (83), and since a PSD of an SDOF is a scalar multiple of the frequency response function squared, the following can be substituted: f 2 = f and f 2 = f : u,H u l,H l ⎛ ⎛ −1 2 4⎞ 1 ( f f )⎞ ⎛( f − f )⎞ − u ζ = − ⎜ 4 + 4⎜ 2 ⎜ ⎝ f ⎝ p l u l ⎟ − ⎜ ⎟ ⎟ . (83) ⎠ ⎝ f ⎠ ⎟ p ⎠ 21

Original page 32

APPENDIX B — FITTING A DECAYING EXPONENTIAL TO OSCILLATION MAXIMA Appendix B shows that the decay rate of a decaying envelope of an SDOF underdamped response is identical to the decay rate of a curve through the SDOF underdamped response local maxima. This allows for determination of damping ratio based on an exponential fit through the local maxima and also shows that the logarithmic decrement method is exact. For an underdamped system, 0 < ζ < 1, the particular solution of equation (1) with pressure nt variable is given by equation (84):7 (p. 245) as the depende ⎡⎛ −ζ ω t ζ 1 n ⎢ p(t) = e ⎜ p(0) + ⎢⎜ 2 ω ⎣⎝ 1−ζ d ⎞ ⎤ p!(0)⎟ sin(ω t) + p(0)cos(ω t)⎥ . (84) d d ⎟ ⎥ ⎠ ⎦ The initial conditions are given as the initial pressure amplitude, p(0), and initial rate of change of pressure, p!(0). Equation (84) can be simplified using a linear combination of sinusoids, 2 2 –1 a + b )), where tan2–1 is the two-argument arctangent and acos(x) + bsin(x) = cos(x – tan2 (b,a covers the entire circular range. By factoring equation (84) in such a way to obtain equation (85), the signum function of the initial pressure amplitude emerges: ⎛ p(0)⋅sgn( p(0)) −ζ ω t n p(t) = ⋅X ⋅e cos⎜ω t 2 d 1−ζ ⎜ ⎝ ⎛ ⎞ ⎞ −1 ζ 1 − tan2 ⎜ p(0) + p!(0), p(0)⎟ ⎟ . (85) ⎜ 2 ω ⎟ ⎟ ⎝ 1−ζ d ⎠ ⎠ However, equation (85) can be further reduced to equation (86) since, in this form, a sign change is nullified by the difference in phase between the two-argument arctangent and the standard arctangent. Thus, the two-argument arctangent reverts back to the standard arctangent upon factoring p(0) out of the 2 2 term. The initial amplitude, p(0), was originally factured out of X to retain a + b p(0) as an independent factor in equation (85). The simplified form is given as equations (86)–(88). Equation (86) was previously written as equation (6): p(0) −ζ ω t −1 p(t) = ⋅X ⋅e 2 1−ζ n cos(ω t − tan (φ)) , (86) d 2 2 p!(0) 2 p(0) X = 1+ 2ζ 1−ζ (0)ω p(0) ω p and 22 ! +(1−ζ ) , (87) 2 2 d d

Original page 33

ζ φ = 2 1−ζ p!(0) + . (88) p(0)ω d The amplitude envelope of equation (86) was given previously in equation (8). For a zero initial rate of change of pressure, commonly examined in most textbooks, equation (86) simply becomes equation (89),7 (p. 245) and then the initial pressure at t = 0 is exactly a local maximum for positive p(0): p(0) −ζ ω t ⎛ −1 2 ⎞ p(t) = ⋅e n cos ω t − tan ζ 1−ζ (89) ( ) . 2 ⎝ d ⎠ 1−ζ It is noted previously that the envelope of a damped free response is given as an exponential decay. As shown in figure 3, the amplitude envelope (red curve), described in equation (8), is actually tangent to the equation time history (blue curve) at locations near the oscillation maxima. It can be shown analytically that the decay rate is identical for the decaying amplitude envelope (red curve) and for an exponential decay through the response maxima (black dashed curve). There is only an amplitude correction needed if the exponential decay envelope is desired. However, since the decay rates are identical, either the envelope or an exponential decay through the maxima can be used to give the damping ratio exactly. The local oscillation maxima of the decaying response function, equation (86), can be found by calculating the derivative and solving the inequality equation (90). This will explicitly give a range where the pressure is increasing beginning at a local minima and ending at a local maxima. ∂ p(t) ∂t 1 Free Decay Response Envelope Exponential Exponential Through Maxima 0.5 0 Amplitude –0.5 ƒn=15 Hz ζ=0.10 p(0)=1.0 –1 0 0.1 0.2 0.3 0.4 0.5 (a) Time (s) (b) > 0. (90) 0.545 Free Decay Response 0.54 Envelope Exponential Exponential Through Maxima 0.535 0.53 Amplitude 0.525 0.52 0.515 0.065 0.066 0.067 0.068 0.069 0.07 0.071 Time (s) Figure 3. Comparison of amplitude envelope to exponential through local maxima: (a) Free decay response, envelope, and exponential through local maxima and (b) zoomed view. 23

Original page 34

Simplifying the derivative, knowing that X ≥ 0 for 0 ≤ ζ ≤ 1, gives the following range of t for n = 0,1,2,... ⎧ ⎧ p(0) > 0 and p! (0) ≥ 0 ⎪ 1 −1 (2n −1)π 1 −1 2n ⎪ ⎪ tan (Ψ ) + < t < tan (Ψ ) + ω ! ! 2 ω ω ω ⎪ d d d ⎪ ⎪ ⎪ π , ⎨ d ⎪⎩ p(0) < 0, p (0) > 0, and p(0) / p (0) ≥ −ζ 1−ζ ω d ⎧ p (0) < 0 and p! (0) ≤ 0 ⎪ ∂ p(t) ⎪ 1 −1 2nπ 1 −1 (2n +1) ⎪ 2 > 0 when ⎨ tan (Ψ ) + < t < tan (Ψ ) + ∂t ⎪ ω ω ω ω d d d d ⎪ ⎪ ⎪ π , ⎨ p(0) < 0, p! (0) > 0, and p(0) / p! (0) < −ζ 1−ζ ω d ⎪ 2 ⎪ p(0) > 0, p! (0) < 0, and p(0) / p! (0) > −ζ 1−ζ ω ⎩ d 1 −1 (2n +1)π 1 −1 (2n + 2) 2 ⎪ tan (Ψ ) + < t < tan (Ψ ) + ⎪ ω ω ω d d d ⎪ ⎩ π , { p(0) > 0, p! (0) < 0, and p(0) / p! (0) ≤ −ζ 1−ζ ω . d ω d (91) The phase angle, Ψ, in equation (91) is obtained when simplifying the derivative using the arctangent difference trigonometric identity: tan–1(a) – tan–1(b) = tan2–1(a – b, 1 + a · b) and is given as equation (92): 2 φ 1−ζ −ζ Ψ = = 2 1−ζ + φζ p!(0) 2 p!(0)⋅(1−ζ ) . (92) 2 ζ 1−ζ + p(0)ω d The standard arctangent is used in equation (91), as appropriate conditions and periodicity are incorporated throughout. From the ranges of an increasing response given in equation (91), the local oscillation maxima can be determined to occur at tmax,n for n = 0,1,2,... ⎧ ⎧ p(0) > ⎪ 1 −1 2nπ ⎪ ⎪ tan (Ψ ) + , ⎨ 0 and p!(0) ≥ 0 2 ω ω ⎪ p(0) < 0, p(0) > 0, and p(0) / p(0) ≥ −ζ 1−ζ ω ⎪ d d ⎩ ⎪ ! ! d ⎪ ⎧ p(0) < 0 and p(0) ≤ 0 ⎪ ⎪ ⎪ 1 −1 (2n +1)π ⎪ t = ⎨ tan (Ψ ) + , ⎨ p(0) < max,n ⎪ ω ω ⎪ d d ⎪ ⎪ p(0) > ⎩ ⎪ ⎪ 1 (2n + 2)π −1 ! 2 0, p!(0) > 0, and p(0) / p!(0) < −ζ 1−ζ ω d 2 0, p!(0) < 0, and p(0) / p!(0) > −ζ 1−ζ ω d 2 { 0, p!(0) < 0, and p(0) / p!(0) ≤ −ζ 1−ζ ω . (93) ⎪ tan (Ψ ) + , p(0) > ⎪ ω ω d d ⎪ ⎩ 24 d

Original page 35

Substituting equation (93) into equation (86) and simplifying gives the values at the maxima for n = 0,1,2,... ⎧ ζ −1 ⎪ − (2⋅n⋅π +tan (Ψ )) p(0) 2 ⎪ 1−ζ −1 ⎧ p(0) > 0 and p! (0) ≥ 0 −1 ⎪ ⋅X ⋅e cos 2(⋅n⋅π +tan (Ψ ) −tan (φ)) , ⎨ ⎪ 2 1−ζ ⎪ ⎪ ⎪ ζ −1 ⎪ − ((2⋅n+1)⋅π +tan (Ψ )) ⎪ p(0) 2 1 −1 ⎪ 2 1−ζ − p t( ) = ⎨ ⋅X ⋅e cos((2⋅n+1)⋅π +tan max, n 2 ⎪ 1−ζ ⎪ ⎪ ⎪ ζ −1 ⎪ − ((2⋅n+2)⋅π +tan (Ψ )) p(0) 2 2 ⎪ p(0) < 0, p! (0) > 0, and p(0) / p! (0) ≥ −ζ 1−ζ ω ⎩ d ⎧ p (0) < 0 and p! (0) ≤ 0 ⎪ (Ψ ) −tan (φ)) , ⎨ p(0) < 0, p! (0) > 0, and p(0) / p! (0) < −ζ 1−ζ ω d ⎪ 2 ⎪ p(0) > 0, p! (0) < 0, and p(0) / p! (0) > −ζ 1−ζ ω ⎩ d ⎪ 1−ζ −1 −1 2 ⋅X ⋅e cos((2⋅n+ 2)⋅π +tan ⎪ 2 1−ζ ⎪⎩ (Ψ ) −tan (φ)) , { p(0) > 0, p! (0) < 0, and p(0) / p! (0) ≤ −ζ 1−ζ ω . d (94) The arctangent functions are not combined into the two-argument arctangent in equation (94); −1 −1 −1 however, the expression tan2 (Ψ − φ, 1+ φΨ ) = tan (Ψ ) − tan (φ) arctangent terms. can be used to replace the The value of the amplitude envelope at this time is different from the value at the response function at the maxima and found by substituting equation (93) into equation (8) for n = 0,1,2,... ⎧ ζ −1 ⎪ − (2⋅n⋅π +tan (Ψ )) p(0) 2 ⎪ 1−ζ ⋅X ⋅e ⎪ 2 1−ζ ⎪ ⎪ ⎪ ζ −1 ⎧ p(0) > 0 and p! (0) ≥ 0 ⎪ , ⎨ 2 ⎪ p(0) < 0, p! (0) > 0, and p(0) / p! (0) ≥ −ζ 1−ζ ω ⎩ d ⎧ p (0) < 0 and p! (0) ≤ 0 ⎪ − ((2⋅n+1)⋅π +tan (Ψ )) ⎪ ⎪ p(0) 2 1−ζ A t( ) = p t( ) = ⎨ ⋅X ⋅e max, n max, n 2 ⎪ 1−ζ ⎪ ⎪ ⎪ ζ −1 ⎪ − ((2⋅n+2)⋅π +tan (Ψ )) p(0) 2 1−ζ ⎪ ⋅X ⋅e ⎪ 2 1−ζ ⎪⎩ ⎪ 2 , ⎨ p(0) < 0, p! (0) > 0, and p(0) / p! (0) < −ζ 1−ζ ω d ⎪ 2 ⎪ p(0) > 0, p! (0) < 0, and p(0) / p! (0) > −ζ 1−ζ ω ⎩ d 2 , { p(0) > 0, p! (0) < 0, and p(0) / p! (0) ≤ −ζ 1−ζ ω . d (95) 25

Original page 36

There are several ways to show that both the decaying response amplitude envelope and the decaying response through the oscillation maxima follow the same decay. The simplest approach is to examine the ratio of the amplitude envelope and the decaying response at all the times the decaying response has a maximum. This ratio is equal to equation (96) for n = 0,1,2,... ⎧ ⎪ −1 −1 ⎪ sec 2( ⋅ n ⋅π + tan (Ψ ) − tan (φ)) , ⎪ ⎪ ⎪ A t ⎪ ( max, n) ⎪ −1 −1 = ⎨ sec((2 ⋅ n +1)⋅π + tan (Ψ ) − tan (φ)) , p t( ) ⎪ max, n ⎪ ⎪ ⎪ −1 −1 ⎪ sec((2 ⋅ n + 2)⋅π + tan (Ψ ) − tan (φ)) , ⎪ ⎪ ⎩ ⎧ p(0) > 0 and p!(0) ≥ 0 ⎪ ⎨ 2 ⎪ p(0) < 0, p!(0) > 0, and p(0) / p!(0) ≥ −ζ 1−ζ ω ⎩ d ⎧ p(0) < 0 and p!(0) ≤ 0 ⎪ ⎪ 2 ⎨ p(0) < 0, p!(0) > 0, and p(0) / p!(0) < −ζ 1−ζ ω d ⎪ 2 ⎪ p(0) > 0, p!(0) < 0, and p(0) / p!(0) > −ζ 1−ζ ω ⎩ d 2 { p(0) > 0, p!(0) < 0, and p(0) / p!(0) ≤ −ζ 1−ζ ω . d (96) For all the maxima, i.e., all the integer values of n, with some manipulation, this simplifies to equation (97): A t( ) max, n = . (97) p t( ) max, n 1 2 1−ζ The ratio is a constant for a given damping ratio and not dependent on frequency. This confirms that the exponential decay rate for the amplitude envelope is identical to the decay rate of the response maxima. Equation (97) also provides an amplitude correction factor to obtain the amplitude envelope if a fitting expression is found based on the local maxima. An exponential fitting exercise comparing the amplitude envelope and the exponential decay through the oscillation maxima also verifies this solution. 26

Original page 37

For completion, the time of tangent intersection is found by equating equations (8) and (86) and solving for the time for n = 0,1,2,... ⎧ ⎧ (0) > 0 and p!(0) ≥ 0 p ⎪ 1 −1 2nπ ⎪ ⎪ tan (φ) + , ⎨ 2 ω ω ⎪ p(0) < 0, p(0) > 0, and p(0) / p(0) ≥ −ζ 1−ζ ω ⎪ d d ⎩ ⎪ ! ! d ⎪ ⎧ p(0) < 0 and p(0) ≤ 0 ⎪ ⎪ ⎪ 1 −1 (2n +1)π ⎪ ! 2 ⎨ ⎨ (0) < 0, p!(0) > 0, and p(0) / p!(0) < −ζ 1−ζ ω t = tan (φ) + , p intersection, n ⎪ ω ω ⎪ d d (0) > 0, p!(0) < 0, and p(0) / p!(0) > −ζ 1−ζ ω ⎪ ⎪ p ⎩ ⎪ ⎪ 1 −1 (2n + 2)π d 2 d 2 ⎪ tan (φ) + , p(0) > 0, p(0) < 0, and p(0) / p(0) ≤ −ζ 1−ζ ω . { ⎪ ω ω d d ⎪ ⎩ ! ! d (98) Substituting this into equations (8) or (86) and simplifying gives the values at the intersection for n = 0,1,2,... ⎧ ζ −1 ⎪ − (2⋅n⋅π +tan (φ)) p(0) 2 ⎪ 1−ζ ⋅X ⋅e , ⎪ 2 1−ζ ⎪ ⎪ ⎪ ζ −1 ⎪ − ((2⋅n+1)⋅π +tan (φ)) ⎪ p(0) 2 1−ζ ⎧ p(0) > 0 and p!(0) ≥ 0 ⎪ ⎨ 2 ⎪ p(0) < 0, p!(0) > 0, and p(0) / p!(0) ≥ −ζ 1−ζ ω ⎩ d ⎧ p(0) < 0 and p!(0) ≤ 0 ⎪ ⎪ 2 ( ) = ⎨ , ⎨ p(0) < 0, p!(0) > 0, and p(0) / p!(0) < −ζ 1−ζ ω A t ⋅X ⋅e intersection, n 2 ⎪ 1−ζ ⎪ ⎪ ⎪ ζ −1 ⎪ − ((2⋅n+2)⋅π +tan (φ)) p(0) 2 1−ζ d ⎪ 2 ⎪ p(0) > 0, p!(0) < 0, and p(0) / p!(0) > −ζ 1−ζ ω ⎩ d 2 ⎪ ⋅X ⋅e , p(0) > 0, p(0) < 0, and p(0) / p(0) ≤ −ζ 1−ζ ω . ⎪ 2 1−ζ ⎪⎩ { ! ! d (99) 27

Original page 38

APPENDIX C — OTHER RELATIONSHIPS A damped time harmonic pressure wave can be written as equation (100):2 (p. 17), 5 (p. 212) iω t −i kx −αx iω t −i kx −αct iω t −i kx −βt p(x,t) = P e( e )e = P e( e )e = P e( e )e , (100) where c is the sound speed, P is the pressure amplitude, p is the pressure, x is the spatial coordinate, t is the time coordinate, k is the wave number, and w is the angular frequency. The temporal absorption coefficient in equation (100), b, is defined in equation (10).5 (pp. 8–11,17) This parameter is commonly used in physical models including acoustics and combustion stability for example. It can also be represented as the spatial absorption coefficient by equation (101):3 (p. 299),5 (p. 217) α ⋅c . (101) β = The temporal and spatial absorption coefficient can be written in terms of quality factor using equations (10), (57), and (101):5 (pp. 8–11,17) π f n Q = β π f = n . (102) α ⋅c In terms of damping ratio, using equations (57) and (102): β ζ = 2π f n α ⋅c = . (103) 2π f n These are listed in table 1 and can be compared to experimentally extracted values. 28

Original page 39

REFERENCES 1. Casiano, M.J.; and Zoladz, T.F.: “Theoretical Consolidation of Acoustic Dissipation,” NASA/ TM—2012–217455, NASA Marshall Space Flight Center, Huntsville, AL, pp. 8, March 2012. 2. Munjal, M.L.: Acoustics of Ducts and Mufflers, John Wiley & Sons. Inc., New York, pp. 17, 22–26, 1987. 3. Blackstock, D.T.: Fundamentals of Physical Acoustics, John Wiley & Sons, Inc., New York, pp. 20, 299, 2000. 4. Meirovitch, L.: Fundamentals of Vibrations, Waveland Press, Inc., Long Grove, IL, pp. 87–98, 110–120, 2010. 5. Kinsler, L.E.; Frey, A.R.; Coppens, A.B.; and Sanders, J.V.: Fundamentals of Acoustics, John Wiley & Sons, Inc., New York, pp. 8–11, 16, 17, 212, 217, 2000. 6. Thomson, W.T.: Mechanical Vibrations, Second Edition, Prentice-Hall, Inc., Englewood Cliffs, NJ, pp. 52, 53, 55, 56, 124–126, 1956. 7. Ogata, K.: System Dynamics, 3rd ed. Prentice-Hall, Inc., Upper Saddle River, NJ, pp. 245, 249, 1998. 8. Thrane, N.; Wismer, J.; Konstantin-Hansen, J.H.; and Gade, S.: Application Note - Practical use of the “Hilbert transform,” Brüel & Kjær, Denmark, 1995. 9. Bendat, J.S.; and Piersol, A.G.: Engineering Applications of Correlation and Spectral Analysis, 2nd ed., John Wiley & Sons, Inc., New York, pp. 18–21, 124–126, 1993. 10. Kenny, R.J.; Lee, E.; Hulka, J.R.; and Casiano, M.J.: “Signal Processing Methods for Liquid Rocket Engine Combustion Spontaneous Stability and Rough Combustion Assessments,” 60th JANNAF Propulsion Meeting and 7th Liquid Propulsion Subcommittee Meeting, Colorado Springs, CO, December 5–9, 2013. 11. Ibrahim, S.R.: “Modal Identification of Structures from the Responses and Random Decrement Signatures,” NASA-CR-155321, NASA Langley Research Center, Hampton, VA, p. 20, 1977. 12. Vandiver, J.K., Dunwoody, A. B., Campbell, R. B., and Cook, M. F.: “A Mathematical Basis for the Random Decrement Vibration Signature Analysis Technique,” Journal of Mechanical Design, Vol. 104, pp. 307–313, April 1982. 13. Brand, O.; Dufour, I.; Heinrich, S.M.; and Josse, F. (Eds.): “Resonant MEMS: Fundamentals, Implementation, and Application,” Advanced Micro & Nanosystems, Wiley-VCH Verlag & Co. KGaA, Weinheim, Germany, pp. 3–13, 2015. 29

Original page 40

REPORT DOCUMENTATION PAGE Form Approved OMB No. 0704-0188 The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operation and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To) 01–09–2016 Technical Memorandum 4. TITLE AND SUBTITLE Extracting Damping Ratio From Dynamic Data and Numerical Solutions 6. AUTHOR(S) M.J. Casiano 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) George C. Marshall Space Flight Center Huntsville, AL 35812 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) National Aeronautics and Space Administration Washington, DC 20546–0001 12. DISTRIBUTION/AVAILABILITY STATEMENT Unclassified-Unlimited Subject Category 71 5a. CONTRACT NUMBER 5b. GRANT NUMBER 5c. PROGRAM ELEMENT NUMBER 5d. PROJECT NUMBER 5e. TASK NUMBER 5f. WORK UNIT NUMBER 8. PERFORMING ORGANIZATION REPORT NUMBER M–1418 10. SPONSORING/MONITOR’S ACRONYM(S) NASA 11. SPONSORING/MONITORING REPORT NUMBER NASA/TM—2016–218227 Availability: NASA STI Information Desk (757–864–9658) 13. SUPPLEMENTARY NOTES Prepared by the Propulsion, Structural, Thermal and Fluid Analysis Division, Engineering Directorate 14. ABSTRACT Damping can be extracted from test data or be obtained through theoretical modeling. Six methods of extracting damping from data are described: (1) the half-power method, (2) logarithmic decrement method, (3) an autocorrelation/power spectral density fitting method, (4) a frequency response fitting method, (5) a random decrement fitting method, and (6) a newly developed half-quadratic gain method. Additionally described are state-space models (and finite element method models) which provide a theoretical damping via complex frequency. Relationships are also provided for several other damping parameters in terms of a single damping parameter which can aid in comparing data to theoretical models of a system. 15. SUBJECT TERMS extracting damping, extracting damping ratio, damping, damping ratio, damping rate, combustion stability, combustion stability model validation 16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF ABSTRACT 18. NUMBER OF 19a. NAME OF RESPONSIBLE PERSON a. REPORT b. ABSTRACT c. THIS PAGE U U U UU 30 PAGES STI Help Desk at email: help@sti.nasa.gov 42 19b. TELEPHONE NUMBER (Include area code) STI Help Desk at: 757–864–9658 Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std. Z39-18

Original page 41

Original page 42

National Aeronautics and Space Administration IS02 George C. Marshall Space Flight Center Huntsville, Alabama 35812