Center Impedance Method for Estimating Complex Modulus with Wide Frequency Range and Large Loss Factors

A simple method for determining viscoelasticity over a wide frequency range using the frequency response function (FRF) mobility obtained by the center impedance method is presented. As user data comprise the FRF between the velocity of the excitation rod and excitation force, it is challenging to separate the signal and noise. Our proposed method is based on the FRF obtained from the analytical solution of the equation of motion of the viscoelastic beam and relationship between the complex wavenumber (real wavenumber and attenuation constant) of flexural wave and viscoelasticity. Furthermore, a large loss factor can be handled over a wide frequency range without using the half-power bandwidth. In this study, actual FRF mobility data containing noise were processed using preprocessing, inverse calculation, and postprocessing. Preprocessing removed low-coherence data, compensates for the effects of instrument gain, and transformed the FRF into its dimensionless equivalent. Then, inverse calculations were used to solve the mobility equation and determine the complex wavenumber. In postprocessing, the complex wavenumber obtained by the inverse calculation was curve fitted using functions with mechanical significance. Consequently, the storage modulus based on the curve-fitted complex wavenumber was a monotonically increasing frequency function. The loss factor had a smooth frequency dependence such that it has the maximum value at a single frequency. The proposed method can be applied to composite materials, where the application of time-temperature superposition is challenging. We utilized the measured FRF mobility data obtained over a duration of several seconds, and this method can also be applied to materials with large loss factors of 1 or more.


Introduction
Acoustic radiations from the surface of plates at the floor of a car decrease with its bending vibrations, which is effective for reducing noise from interiors of automobiles, railway trains, and aircrafts. In addition, vibration damping is necessary for passenger comfort. Recently, polymeric composites have been widely used for this purpose, and therefore, vibration damping properties of various plate materials must be evaluated. To this end, tests are performed at the point of beam vibrations, and the storage modulus and loss factor of the specimen are estimated from the responses.
e Oberst beam method is used to measure the responses of composite beams with the material attached to an aluminum or steel substrate for damping materials with low storage moduli [1,2]. According to ASTM E756, the Oberst beam is magnetically vibrated, and the beam responses are measured using a noncontact magnetic transducer. However, the beam center can be vibrated without using a magnetic substrate, and this method is called the center impedance or central excitation method as specified in JIS K 7391 [2]. Although this is increasingly used in Japan than elsewhere, its usage worldwide is expanding because of the relatively convenient installation of the equipment [3]. In the center impedance method, the vibration damping property of the specimen is evaluated using the transfer function between the specimen displacement (velocity or acceleration) and vibration force applied to the specimen center. A piezo element and force gauge are constructed into the impedance head attached to the tip of the vibration rod to measure the vibration force and acceleration.
In the extant measurement standards [1,2], Young's modulus for a specimen regarded as a purely elastic material is estimated using the resonant frequency. However, that of a viscoelastic material varies because of its large loss factor.
ese standards indicate that the half-power bandwidth can be used to estimate this factor near the resonant frequency when assuming structural loss factors to be low. When the half-power bandwidth is used, the system comprising viscoelastic beams is approximated by an equivalent lumped constant system consisting of a mass point, spring, and dashpot. Here, it is assumed that the square of the damping ratio is much less than 1 and the structural loss factor is twice the damping ratio. However, in the half-power bandwidth method, a large damping ratio cannot be handled because the second and higher order terms of the damping ratio are ignored. Papagiannopoulos and Hatzigeorgiou [4] increased the upper limit of the damping ratio to ensure that the halfpower bandwidth can be applied by incorporating higher order terms; however, only the loss factor near the resonant frequency could be obtained, which is a limitation. Wojtowicki et al. [5] noted that the frequency response function (FRF) of the displacement amplitude expressed as a function of a complex wavenumber contains all requisite information regarding viscoelasticity, and their method does not require the mode-dependent half-power bandwidth. e equivalent viscoelasticity of the Oberst beam was associated with the complex modulus of the damping layer [5]; however, the complex wavenumber and viscoelasticity were not explicitly related and numerical processing using the nonlinear leastsquares method was required.
For determining the viscoelasticity for a beam specimen conventionally, Young's modulus and loss factor of the composite beam are estimated at a certain resonant frequency. In the Oberst beam method, the FRF is measured at a constant temperature in an environmental chamber and the viscoelasticity of the damping layer is determined in several incremental steps for each measurement under a constant temperature. Furthermore, the frequency dependence of the viscoelasticity is estimated at an ambient temperature of 20°C by applying a time-temperature superposition (TTS) to the results [5]. However, applying TTS to composite materials having complex configurations is challenging. As an alternative, the FRF obtained from the analytical solution is employed for the equation of motion of the viscoelastic beam over a wide frequency range. is could be applicable to the FRF used by Wojtowicki et al. [5]. However, they applied TTS to the viscoelasticity of polyvinyl chloride (PVC) obtained from the FRF of the Oberst beam near the resonance point for determining the frequency dependence of viscoelasticity at ambient temperature. e material viscoelasticity has been identified previously using the relationship between the complex wavenumber and viscoelasticity in applications other than vibration damping tests. Benatar et al. [6] derived an analytical solution of the Pochhammer frequency equation for complex wavenumbers in a split Hopkinson pressure bar for wavelengths comparable to its diameter of the bar. ey identified viscoelasticity by analytically solving the Pochhammer frequency equation for small damping, while Zvietcovich et al. [7] determined it by optically measuring the responses of acoustically vibrating living tissues such as the cornea using Gaussian beam pulses. Owing to the complexity of the governing equations of wave propagation in these studies, the relationship between complex wavenumbers and frequencies was simplified using linearization. Conversely, the solution must be strictly managed to express the analytical solution of the flexural wave in the thin beam specimen in a simple form. We demonstrated that, for thin viscoelastic rods and beams (plates), a simple relationship between the viscoelasticity and complex wavenumber exists [8,9]. Upon combining these with an analytical solution for the FRF (mobility (V/F)) between the velocity and excitation force, the viscoelasticity of the material can be determined over a wide frequency range from FRF data using complex wavenumbers [9]. e remainder of this paper is organized as follows. In Section 2, we outline the measurements using the center impedance method. In Section 3, the complex wavenumber for the flexural wave in the viscoelastic beam is expressed as an explicit function of viscoelasticity; the inverse formula for determining the viscoelasticity for a given complex wavenumber (real wavenumber and attenuation constant) is displayed along with the FRF mobility expressed as an analytical function of the complex wavenumber. Consequently, the complex wavenumber is obtained from the FRF by solving this complex function in reverse. e viscoelasticity can be determined at an arbitrary frequency in the measurement range via the complex wavenumber obtained from the measured FRF mobility data and the above analytical treatment [8,9]. e complex wavenumber and viscoelasticity (storage modulus and loss factor) were determined for arbitrary frequencies by inverse calculations using dimensionless mobility data obtained from preprocessing. e obtained viscoelasticity has a one-to-one correspondence with the measured FRF mobility; however, it considerably fluctuates with frequency because of the noise in the measured data. us, noise from the FRF was removed by fitting the relationship between the complex wavenumber and frequency with functions of mechanical significance wherein the real wavenumbers were treated as monotonously increasing frequency functions. e ratio of attenuation constant to real wavenumber was expressed as a function of the loss factor and fitted as that of the frequency using a regression curve proportional to the density function of the gamma distribution. A relationship similar to this curve can be observed in the viscoelastic model of Córtes and Elejabarrieta [10]. Finally, the FRF mobility was calculated using the curve-fitted complex wavenumber as compared with the measured results. e application of the inverse complex wavenumber formula on the curvefitted complex wavenumber transformed the viscoelasticity (storage modulus and loss factor) to a smooth frequency function.

Experimental System
A schematic of the measurement system configuration is presented in Figure 1; the beam center is vibrated with white noise. As per convention, the measurement system was used to evaluate the loss factor using the half-power bandwidth method and estimate the elastic modulus based on the resonant frequency. However, FRF was used differently. Although the FRF mobility (V/F) was calculated based on the measured acceleration and excitation force of the excitation rod, it was expressed as an analytical function of a complex wavenumber using the beam equation of motion, as discussed in Section 3.

Identification of Viscoelasticity by Analytical Solutions
is section briefly summarizes the descriptions [8,9] regarding the governing equations and analytical expression of FRF mobility.

Inverse
Equations. e following inverse equations were obtained from equations (4) and (5) [9]:  Using the inverse equation (equation (6a)), the storage modulus E and loss factor η were obtained for a given complex wavenumber; equation (6b) denotes conditions for the wave solution (Appendix A).

FRF Mobility at Excitation Point.
Using the solution of the equation of beam motion, the dimensionless FRF mobility W can be expressed as [8,9] where M * ob denotes the dimensional FRF mobility (velocity amplitude at the excitation point)/(amplitude of the excitation force) (ms − 1 N − 1 ), F x�L/2 � (EI) * (z 3 y d /zx 3 ) x�L/2 represents the shearing force on the rod (N), L denotes the beam length (m), and y d0 denotes the amplitude of rod displacement (m) [9].
W defined using equation (7a) can be expressed as a function of the dimensionless complex wavenumber Z [9] as follows.

Dimensionless Mobility Equation
Equations (8b)-(8d) are the dimensionless complex wavenumber, real wavenumber, and attenuation constant ×(− 1), respectively. Equation (8e) represents the constraints, where the first expression is an additional condition for excluding evanescent waves. Figure 2 portrays the behavior of the complex function W with Y as a parameter, where the following properties of W are portrayed as well [9].

Properties of Dimensionless Mobility (W)
(a) If viscoelasticity is constant, the dimensionless real wavenumber X increases with frequency. (b) e X-axis is delimited with resonance points at intervals of approximately π. (c) e phase angle of W ranges from 0 to π (180°). (d) e relationship between the phase angle of dimensional mobility M * ob and W is given by arg(M * ob ) � arg(W) + π/2. (e) e phase angle increases from a frequency slightly lower than the antiresonance point and starts to decline from that slightly lower than the resonance point. Although these alterations are stepped at small attenuation, they become gentle as the attenuation increases.
When the FRF mobility is measured, the real wavenumber β R and attenuation constant α are obtained for each frequency via the inverse calculation of equation (8a) for Z. e ramp response method developed in a previous study was used for the inverse calculation [9] (Appendix B). In general, properties (b) and (e) are vital for examining the dimensionless mobility data. e real part of Z must be between the dimensionless real wavenumbers for the n R− 1 th and n R th resonance points. To ensure this, the information regarding the closest resonance point on the high-frequency side of the measurement frequency f designated as the n R th resonance point is used as the table value. Furthermore, viscoelasticity (E, η) was determined for each frequency by applying the inverse formula to the complex wavenumber data.

Resonance.
e dimensionless real wavenumber of the n R th resonance point is the n R th root of tan X R + tanh X R � 0.
3.2.4. Antiresonance. e dimensionless real wavenumber of the n th antiresonance point represents the n th root of the following equation. In general, almost no vibration is seen near the exciting rod; however, the free end vibrates vigorously at the antiresonant frequency.
We calculated the dimensionless mobility W of the Oberst beam ( Figure 3) based on the viscoelasticity of PVC measured by Wojtowicki et al. [5]. ereafter, we regarded this as virtual measurement data and performed inverse calculations [9]. e obtained complex wavenumber was applied to equations (6a)-(6b) for obtaining the storage modulus and structural loss factor of the composite beam. As presented in Figure 4, the virtual measurement data ("measured") for viscoelasticity were reproduced via inverse calculations ("estimated"). In their study [5], a 6 mm-thick and 400 mm-long PVC-based material was glued onto a 1.6 mm-thick aluminum substrate.

Results and Discussion: Actual Measurements
and Data Processing e measured FRF mobility (V/F) data contained noise and were not smooth, as illustrated in Figure 3. To develop a denoising method as discussed, inverse calculation was applied to these data, which resulted in unrealistic viscoelasticities with large fluctuations of frequency.   Shock and Vibration

Data Processing Flow and Noise
Removal. Figure 5 depicts the data processing flow that comprises (a) preprocessing, (b) inverse calculation, and (c) postprocessing.
In the first step, measured data were normalized and converted to dimensionless mobility values W. Second, the complex wavenumber was obtained using equations (8a)-(8e) with respect to Z. Subsequently, the viscoelasticity was identified from the complex wavenumber using equations (6a) and (6b). If W contained noise, the relationship between the estimated viscoelasticity and frequency becomes unsmooth and unrealistic. Finally, the relationship between the complex wavenumber and frequency was curve fitted with functions of mechanical significance instead of directly eliminating noise from the viscoelastic data. Moreover, equations (8a)-(8e) were applied to the curve-fitted complex wavenumber to obtain the dimensionless mobility W fit as a smooth frequency function. However, the smooth relationship between viscoelasticity (E fit , η fit ) and frequency can be found by applying equations (6a) and (6b) to the curvefitted complex wavenumber. (1). Data points with a coherence of ≤0.99 were first removed from the measured FRF mobility data. In addition, the term "coherence" is defined and summarized in the following section. Furthermore, the phase angle range of the mobility M * ob was normalized within 90°-270° [9]. e preprocessed FRF mobility data are presented in Figure 6, wherein the phase angle data include spike-like noise represented by dashed circles. Similar to the properties of W, the phase angle increased from a frequency slightly under the antiresonance point and started declining before the resonance point. A few phase angle data that cannot be normalized within the range of 90°-270°are present, which were subjected to the following equations:

Preprocessing
In these cases, ε L , ε H > 0 represent small positive numbers and phase Mob denotes the phase angle of dimensional mobility (°) .

Definition of Coherence.
e term "coherence" refers to the output value of the measurement system and is defined as where f denotes the frequency (Hz) and S VV (f) and S FF (f) represent the power spectra of the excitation point velocity and excitation force, respectively. Furthermore, S VF (f) represents the cross spectrum between the excitation point velocity and excitation force. e units of (S VV (f) S FF (f)) 1/2 and S VF (f) are identical to each other such that coherence is a dimensionless number with values in the range of 0-1. Coherence has a meaning similar to the correlation coefficient in the present context. For a small coherence, the velocity of the excitation point is considered to be noise independent of the excitation force.

Preprocessing (2): Gain Correction.
e FRF mobility data are influenced by the gain of the measurement system, which is estimated in equation (12). Parameters f 0 and k G in equation (12) were determined by employing the leastsquares method using the analytical solution and measured mobility data for an aluminum beam.
where f denotes the frequency (Hz), f 0 (>0) is a constant (Hz), and k G (>0) is a constant (-). Using this gain, the measured dimensional mobility M * obM at f was corrected using   Shock and Vibration where equation (13b) represents the absolute mobility in decibels.

Determination of Gain Parameters by Measuring an
Aluminum Beam. Assuming that the gain of the measurement system is independent of the beam material, the parameters were estimated using the measured results of the aluminum beam. Consequently, f 0 � 10 (Hz) and k G � 0.15 (− ) were obtained. e impact of gain correction on the absolute value of the dimensionless mobility data of the aluminum beam is presented in Figure 7. For appropriate gain correction, a dimensional mobility of 0 dB corresponds to 1 ms − 1 N − 1 . In the least-squares method, the residual S res was minimized: where u � R e W, v � I m W, and C, M, and n represent calculated, measured, and n-th frequency, respectively. (g) Calculation of viscoelasticity (E, η) from the complex wavenumber using inverse formula of the complex wave number formula.
Postprocessing (Removing noise from the complex wavenumber data) (h) Curve fitting of the relationship between complex wavenumber β* = β Rand frequency with smooth functions having mechanical meaning (i) Calculation of FRF mobility from the curve-fitted complex wave number data (j) Estimation of viscoelasticity (E, η) from the curve-fitted complex wavenumber data

Preprocessing (3): Rendering FRF Mobility Dimensionless.
As mobilities after preprocessing steps (1) and (2) are displayed in decibels, the absolute value |M * ob | (ms − 1 N − 1 ), real part R e (M * ob ), and imaginary part I m (M * ob ) of the dimensional mobility can be described using where phaseMob denotes the phase angle of the dimensional mobility (°) .

Dimensionless Mobility.
e measured data M * ob of the FRF mobility can be converted to dimensionless mobility W using where ρbhL denotes the mass of the beam specimen. e phase angle arg(W) of the dimensionless mobility portrayed in Figure 8 ranges from 0 to π.

Inverse
Calculation. e application of inverse calculations expressed in equations (8a)-(8e) to W in Figure 8 yielded the relationship between the dimensionless complex wavenumber and frequency, as indicated in Figure 9. e thick lines represent the dimensionless real wavenumber X, while the thin lines represent the dimensionless attenuation constant (− Y). A drastically steep rise can be observed in the dimensionless real wavenumber X near 5.5 kHz in Figure 9(a) (ABS) and 6 kHz in Figure 9(b) (PP). However, the dimensionless attenuation constant alternated between small and large values. e storage modulus and loss factor were obtained by applying equations (6a) and (6b) to the real wavenumber β R � 2X/L and attenuation constant α � − 2Y/L obtained from the results in Figure 9. ese results are presented in Figure 10, where the elastic modulus of an elastic body with respect to the resonant frequency is "purely elastic." In addition, both the storage modulus and loss factor fluctuated significantly with the frequency, which is an anomaly for viscoelastic data. Moreover, the loss factor near the resonant frequency obtained using the halfpower bandwidth method was calculated using the measurement system and presented as the "half-power bandwidth."

Postprocessing: Noise Removal by Fitting Complex
Wavenumber vs. Frequency Relationship to a eoretical Curve. Instead of directly eliminating the noise from the viscoelasticity vs. frequency data, we developed a method to fit a curve for the relationship between the complex wavenumber and frequency using smooth functions of intermediate data.

Real Wavenumber vs. Frequency Relationship.
e generalized Maxwell model is applicable to storage modulus E of materials, which is a monotonically increasing function of frequency. As ω ⟶ ∞ at high frequencies, E approaches the instantaneous modulus E inst , and as ω ⟶ 0 at low frequencies, it approaches the long-term modulus E ∞ . If η has a small frequency dependence, equation (4) can be considered for the following relationship between the real wavenumber β R and radial frequency ω.

Shock and Vibration
erefore, if the relationship of l n β R vs. l n ω is plotted, the slope at low frequencies is <1/2 and approaches 1/2 at high frequencies. us, we can determine a curve-fitting equation such that the slope satisfies equations (18a) and (18b).

Equation to Fit the Relationship between Real Wavenumber vs. Frequency
Slope from Reference Point. e plot of l n β R vs. l n ω has many irregularities. Instead of the gradient z(l n β R )/z(l n ω) at each frequency, we define the slope k f (ω) of the straight line connecting the points (l n ω, l n β R ) and (l n ω 0 , l n β 0 ) from the reference point using the following equations: where where ω 0 and β 0 denote the reference radial frequency (s − 1 ) and reference wavenumber (m − 1 ), respectively; furthermore, β R � β 0 for ω � ω 0 .
where E 0 denotes the storage modulus at ω 0 . From equations (18a) and (18b), the slope k f (ω) satisfies as ω ⟶ ∞. (20) Furthermore, ω 0 was (A) Located near the center of the frequency range of interest and (B) Situated far from the peaks and valleys due to resonance points, antiresonance points, and noise of the dimensionless mobility W e slope k f (ω) does not always satisfy equation (20) due to variations in the real wavenumber β R from the inverse calculation on the FRF data. erefore, the following conditions were used to satisfy equation (20): ⟵ k min k min denotes a small positive number . (21)

Curve-Fit
Formula for Slope k f (ω). k f (ω) obtained using the sigmoid function, a smooth function with a value between 0 and 1, was adopted as the curve-fitting equation satisfying equation (20). where where x F d denotes the difference in the logarithm of frequency (logarithm of frequency ratio) (-), ω I denotes the radial frequency of the inflection point (s − 1 ), and y s (x F d ) represents the sigmoid function of x F d .
x F � l n ω, e method for identifying ω I is described in Appendix C. Evaluating the sigmoid function is not necessary once ω I has been calculated. Moreover, equations (22a)-(22d) can be simplified as Applying the slope of equation (23) to equation (19a), the real wavenumber β R becomes a smooth function of the frequency. Figure 11 presents examples depicting the relationship between logarithms of the dimensionless real wavenumber and frequency. e real wavenumber β R was converted to a dimensionless real wavenumber X � β R L/2. In Figure 11, the points at which the thick line deviated from the regression curve were considered to be affected by noise. Vertical and horizontal axes of Figure 11 were rewritten in Figure 12 in linear forms. Here, for both ABS and PP, the thick and dashed lines were in good agreement for frequencies of 4 kHz or less.

Equation to Fit Relationship between Attenuation
Constant and Frequency. In Figure 10, the loss factor obtained by the inverse calculation from FRF data fluctuates in the range from 0 to approximately 0.3. However, it must satisfy conditions (A), (B), and (C). In addition, the ratio α/β R of the attenuation constant to real wavenumber was fitted to a theoretical curve based on these conditions. ereafter, E and η were calculated from the values of β R and α using equations (6a) and (6b).
(A) Relationship between Ratio α/β R and Loss Factor η. From the complex wavenumber equations (4) and (5), α/β R can be expressed as a function of only the loss factor η, which is expressed as

(B) Frequency at Which Ratio α/β R and Loss Factor η Are
Maximized. e following relationship holds from equations (6a) and (6b): where Equations (25a) and (25b) maximize (α/β R ) at the frequency corresponding to the maximization of the loss factor η.
(C) For Many Viscoelastic Materials, the Loss Factor η Has the Maximum Value at a Single Frequency. From properties (A)-(C), the relationship of the ratio (α/β R ) � (− Y/X ) vs. ω should be fitted with an appropriate curve to consider the maximum value at a particular ω. In this case, the loss factor in Figure 10 was replaced by a smooth curve with only one maximal peak. Figure 13 shows examples of the relationship of (α/β R ) and frequency for both materials, which shows a large fluctuation of (α/β R ). erefore, we considered the horizontal axis (frequency (Hz)) as a random variable, and the value on the vertical axis was considered as the number of counts in the histogram. us, the curve fit of the data in Figure 13 can be replaced with a case of an appropriate probability distribution applied to the experimental data represented by a histogram. erefore, a function proportional to the density function of the gamma distribution, a continuous distribution similar to the Poisson distribution, was used as the curve-fit function. e proposed equation for the regression curve of (α/β R ) is presented in equations (26a)-(26e), wherein the normalization by the step size r 0 divided into 16 intervals between the minimum and maximum (α/β R ) was adopted in equations (26a)-(26e). e plots obtained using this normalization are presented in Figure 14.

Regression Curve Proportional to Gamma Distribution.
where where g(f) denotes the probability density of the gamma distribution and Γ(p) denotes the gamma function.
where r min and r max represent the minimum and maximum values of α/β R in the frequency range, respectively; r 0 denotes the step width (-) and f � ω/(2π) denotes the frequency (Hz); σ( > 0) and p( > 0) (− ) are the parameters in the gamma distribution; and n 0Rf represents the coefficient for the height of the distribution ((H Z ) − p ), as defined in Considering the thin line in Figure 14 as a histogram with frequency f as a random variable, the sample mean m f and sample standard deviation s f with respect to frequency f can be defined. Parameters p and σ of the gamma distribution in equation (26b) can be expressed by equations related to mean m f and standard deviation s f . e coefficient n 0Rf was determined using equation (27) such that the integral value of the curve-fit equation and plot data is equal.  e sample mean m f and sample standard deviation s f were calculated using equation (29b). In equation (29b), the region where the convergence points of the inverse calculation (markers in Figure 14) are densely present has large weights.

Shock and Vibration
Here, N is similar to the total count but is generally a noninteger real number. e inverse calculation of equation (8a) was performed for 635 frequency points from 30 to 6400 Hz with 10 Hz increments. e number of points where equation (8a) converged was 490 for ABS and 543 for PP, which are represented by markers in Figures 13 and 14. Plot point numbers k are assigned to these convergence points in the ascending order of frequency.

Curve-Fit for Ratio α/β R vs. Frequency.
e thin lines in Figure 14 show the normalized plots of ratio α/β R for ABS and PP. e thick lines are regression curves proportional to the gamma distributions, and the curve-fit parameters are shown in Figure 14. e gamma distribution parameters p are 3.4 and 4.1 for ABS and PP, respectively. e anomalous jumps are circled in Figure 14(a), which appear near the spike-like noises of the phase angles in Figures 6 and 8. Although the thick lines did not apparently correspond with the thin lines well, the curve-fit result was evaluated by the degree to which the original dimensionless mobility could be expressed. e reason why the thin lines (value by direct inverse calculation from FRF data) vary is discussed in Appendix D.

Viscoelastic Model with Frequency Dependence Whose Loss Factor Is Proportional to the Gamma Distribution.
Córtes and Elejabarrieta [10] derived a viscoelastic model represented by equations (30a) and (30b) via a theoretical study of the relaxation function. ey experimentally validated this model using a material in which concrete was impregnated with epoxy resin. Equation (30b) is the case where p � 2 in equations (26a)-(26e). e constant η 0 is 2 × 10 − 2 , and therefore, the approximation of equation (24b) holds for such a small loss factor. Equation (30b) is a special case of our curve-fit equations (26a)-(26e). Shock and Vibration 13 where E ′ (ω) denotes the storage modulus (Pa) and η(ω) denotes the loss factor (-); furthermore, η 0 � 2cω 0 /(E 0 e), in which c denotes the viscosity (Pa · s) and e represents the base of the natural logarithm.

FRF Mobility Based on Curve-Fitted Complex
Wavenumber Data. We applied the curve-fitted data of the dimensionless complex wavenumber Z � X + iY to equations (8a)-(8e) to obtain the dimensionless mobility. e results were then compared with the original FRF data (see Figures 15 and 16). Furthermore, in these figures, the standard deviation StdReg of the difference between the thick and thin gray lines was calculated with respect to the absolute value and phase angle of W. ese are shown as (thick line) ± StdReg values via dashed lines. However, (thick line)− StdReg has a small lower limit to prevent the dashed line from being negative. In these figures, the thin lines are marked with the convergence points of the inverse calculations of equations (8a)-(8e), which are densely present in the frequency range between the antiresonance and resonance points. Here, the dimensionless mobilities based on the measured data (gray lines) and FRF based on the curve-fit of the complex wavenumber (thick lines) agree well. However, in the frequency range that crosses the resonance points, the convergence points of the inverse calculation of equations (8a)-(8e) are sparse, and the deviation between the thick and thin lines is conspicuous. In Figures 15 and 16, the thin gray lines fall within the range of (thick line) ± StdReg in most frequency ranges.

Viscoelasticity with Curve-Fitted Complex
Wavenumbers. We converted the curve-fitted dimensionless complex wavenumbers to dimensional complex wavenumbers β * � β R − iα, where β R � 2X/L and α � − 2Y/L. e results of determining the viscoelasticity by applying equations (6a) and (6b) to these complex wavenumbers are shown in Figures 17 and 18. In these figures, the thin gray lines are the results estimated directly from the FRF data, while the thick lines are based on the curve-fitted complex wavenumbers. Furthermore, the standard deviations of the thick and thin lines for the storage modulus and loss factor were calculated as StdReg. In Figures 17 and 18, the values of (thick line) ± StdReg are drawn as dashed lines. However, (thick line)− StdReg has a small lower limit to prevent negative values. In Figures 17(b) and 18(b), the thick lines have maximum values near 2 kHz and 2.5 kHz, respectively. 4.6.12. Abnormal Data. In Figure 17(b), an alternate long and short dashed line is shown near the resonance point at 5.6 kHz. Near this frequency, the absolute value of mobility in Figure 6(a) rose abnormally from the minimum value. e circles near 1 kHz in Figures 17(a) and 17(b) correspond to the spike-like noise of the phase angles in Figure 6(a).

Effect of Resonance/Antiresonance Points on Fluctuations of Storage Modulus and Loss Factor.
e resonant and antiresonant frequencies of the beam specimen are indicated by markers on the frequency axes of Figures 17  and 18. Here, the estimated storage modulus and loss factor measured via direct inverse calculations from the FRF data (thin lines) are high near the antiresonant frequency and low near the resonant frequency. e method of determining viscoelasticity based on the complex wavenumber is not originally affected by modes such as resonance and antiresonance. However, the inverse calculations for obtaining the complex wavenumbers from the FRF mobility converge well in the frequency range from the antiresonance to resonance point, whereas the number of convergence points becomes sparse in a frequency range exceeding the resonance point. In Figures 6 and 8 for the FRF data, the falling edge of the phase angle near the resonant frequency represents a sudden change close to a step. According to the qualitative study of equations in Appendix D, the inverse calculation converges poorly because |dW/dZ| ∝ 1/η 2 holds true near the resonance point. However, near the antiresonant frequency, noise may occur due to mixing of natural vibrations, which is concerning. e dimensionless mobility W fit obtained by applying the curve-fitted dimensionless complex wavenumber to equations (8a)-(8e) has a value different from the original measurement data W I as shown in Figures 15 and 16. Figures 19 and 20 show the real and imaginary parts of W fit � u fit + iv fit and W I � u I + iv I , respectively. Here, we consider the difference ΔW � W I − W fit as noise in the FRF measurement data.

Estimation of the Magnitude of Noise Contained in the
where u and v denote the real and imaginary parts of the dimensionless mobility W(− ), respectively. Subscripts I and fit denote direct inverse calculation value from the measurement and the curve-fitted complex wavenumber, respectively.

Noise in Dimensionless Complex Wavenumber.
Consider changes in the dimensionless real wavenumber X and dimensionless attenuation constant − Y when the FRF W changes by ∆W. Since W expressed by equation (8a) is differentiable at points other than at the resonance point, the following relationship holds: erefore, the following relationships among the change ∆W of W, that of ∆X of X, and that of ∆Y of Y hold: erefore, where ΔW � Δu + iΔv,

16
Shock and Vibration Equations (34a) and (34b) show that the fluctuation of the FRF affects both the real and imaginary parts of the complex wavenumber. On the contrary, the difference between the nondimensional complex wave number calculated directly from the FRF W measurement data and the curvefitted nondimensional complex wave number can be expressed by Let us consider δX or ∆X as the fluctuation of X and δY or ∆Y as the fluctuation of Y.
Slope dW/dZ. Figure 21 shows the gradient dW/dZ when Z � X fit + i Y fit is applied to equations (D.7a) and (D.7b). To make the figure clear, the vertical axis of Figure 21(a) is cut off by ±200.

Fluctuation of Dimensionless Complex Wavenumber.
ΔX/X fit and δX/X fit are considered the degree of noise in the data of real wavenumber X. From Figures 22(a) and 23(a), |ΔX|/X fit and |δX|/X fit are in the range of about 5-7% in most frequency regions. However, δX/X fit increases sharply above 5.5 kHz in Figures 22(a) and 23(a). is agrees with the fact that the relation of dimensionless real wavenumber X vs. the frequency relationship near 5.5 kHz in Figures 12(a) and 12(b) is not smooth and corresponds to the spike-like noise of the phase angle near 5.5 kHz in Figures 8 and 9. noise magnitude. In Figures 22(b) and 23(b), the ratio of attenuation constant to real wavenumber − Y fit /X fit is shown by dashed lines, which are the envelopes above ΔY/X fit and δY/X fit . For the measurement data of this study, the dimensionless attenuation constant − Y I obtained by direct inverse calculation from the FRF measurement data was small (− Y I ≈ 0), where the frequency was slightly higher than the resonance point. erefore, δY ≈ − Y fit ≥ 0 in this frequency range. From Figures 22(b) and 23(b), the noise magnitudes in the attenuation constant ΔY/X fit and δY/X fit are in the range of approximately 2-4% in most frequency regions below 5.5 kHz.

Effect of Phase Angle Correction on Loss Factor Estimates.
e measurement data (thin line) of the phase angle of W in Figure 16(b) were replaced with data based on the curvefitted complex wavenumber, represented by a dashed line in Figure 24(a). Nonetheless, the absolute value of W remained on the thick line in Figure 24(a). e loss factor estimated from the inverse calculation at this instant is depicted in Figure 24(b), wherein the fluctuation of the thin line is smaller than that in Figure 18

Conclusion
In this study, we presented a method for estimating the frequency dispersion of viscoelasticity from the measured FRF mobility data using the center impedance method for a beam specimen. To remove the noise in the measured data, the relationship between the complex wavenumber and frequency was curve fitted with smooth functions having mechanical meanings. Furthermore, the viscoelasticity was determined from the curve-fitted complex wavenumber. Using the proposed method, viscoelasticity could be determined over a wide frequency range using only the data obtained from FRF measured for several seconds. is proposed method can be utilized to materials comprising polymer alloys, fibers, etc. as it is difficult to expand the estimated frequency range by simply using TTS for these materials. is method can be applied to materials with large loss factors of 1 or more.
Only FRF mobility from existing measurement systems is utilized here. In the measured FRF data, the phase angle drops in steps on the high frequency side of the resonance point. On the contrary, for the FRF analytical solution, this drop has a smooth frequency dependence. However, the analytical solution only deals with bending vibrations, and vibration modes such as torsion can also cause noise if the beam test specimen is not installed properly. We confirmed that our data contains noise near the resonance points and hypothesized that one of the causes is random excitation by white noise, which is obtained by imposing impulses at random time intervals following a Poisson distribution. By adjusting the average and variance of the number of impulses per unit time, we anticipate that the noise in FRF data may be reduced. We would like to discuss and cooperate with measuring instrument manufacturers to test this hypothesis and optimize the excitation force containing white noise.

A. Derivation of the Inverse Formulas for the Complex Wavenumber Formula of Viscoelastic Beams
e inverse formula (equations (6a) and (6b)) is derived here. For the complex wavenumber equations (4) and (5) where P (m − 4 ) represents a variable related to the storage modulus E and Q (− ) a variable related to the loss factor η. e ratio R between the fourth power of the attenuation constant and that of the real wavenumber is expressed as Considering Q > 0, the variable Q can be expressed as Meanwhile, from the second part of equation (A.1), 1/(1 + η 2 ) can be expressed as follows with respect to Q: erefore, the loss factor η can be expressed in terms of the variable Q as  , the following relationship holds: (A.10) is implies that

B. Calculation Method for Determining Dimensionless Complex Wavenumber Z from Dimensionless Mobility W
In the inverse analysis of equations (8a)-(8e), it is necessary to consider the multivalence of W. erefore, the ramp response method developed in this study was used as a numerical calculation method [9].
Calculation procedure of the ramp response method: (A) e initial value of Z is set as Z (0) (B) W and the slope dW/dZ are calculated for the current value of Z using equations (8a)-(8e) and (D.7a) and (D.7b) (C) e temporary increment of W toward the target value is calculated using equations (B.2a)-(B.2c), which has the same form as the explicit difference approximation of the ramp response (D) Z value is updated from Z (n) to Z (n+1) using equation (B.1) (E) e updated W value W (n+1) for Z � Z (n+1) is calculated using equation (8a) By iterating the calculation procedure, it is expected that Z (n) will converge to the final value Z G as W (n) approaches the target value W G .
Calculations used in the ramp response method: .
(B.1) e superscript (n) represents the n th repetition. e temporary increment of the dimensionless mobility W is expressed as where C oefW denotes a small positive number corresponding to the ratio between the time step and time constant in the difference approximation of the ramp response. e real part of the initial value Z (0) must be assigned a value between the dimensionless real wavenumbers for the n R− 1 th and n R th resonance points.

C. Method for Determining Parameter ω I in Equations (23), (24a), and (24b)
e slope k f (ω) defined by equation (20) is curve fitted using equation (23). We show how to determine the parameters ω I in equation (23), which can be rewritten as Originally, x FI � l n ω I was defined as a constant in equation (23). However, if the slope k f (ω) is replaced with the value k fM (ω) from the measured data in equation (C.2), x FI will fluctuate with frequency. erefore, we estimated ω I as ( Figure 25)