JPH11118661A - Vibration characteristics analyzer - Google Patents

Abstract

PROBLEM TO BE SOLVED: To enhance the identification accuracy even if noise is mixed by an arrangement wherein by reducing the frequency response functions corresponding to intrinsic to the number of intrinsic modes, determining a weighting function while simultaneously taking account of a large number of frequency response functions and applying the weighting function to the sum of residual squares. SOLUTION: The vibration characteristics analyzer comprises a load cell 3 for detecting force, an acceleration sensor 4 and an operating unit 5. The load cell 3 is fixed to an impulse hammer 2 for shaking an object 1. The acceleration sensor 4 is fixed to an arbitrary part of the object 1 in order to measure the response to vibration caused by the impulse hammer 2. The operating unit 5 receives output signals from the load cell 3 and the sensor 4 andydetermines frequency response data and then identifies the mode characteristics by the least squares method for minimizing the sum of the residual squares of the frequency response data and the theoretical value of frequency response function. The operating unit 5 also determines a weighting function while taking account of a large number of frequency response functions simultaneously and applies the weighting function to the sum of residual squares.

Description

DETAILED DESCRIPTION OF THE INVENTION

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a vibration characteristic analyzing apparatus, and more particularly to a vibration characteristic analysis for obtaining a response by exciting a test object such as an engine block and identifying a mode characteristic of the test object from a response function. It concerns the device.

In order to analyze the dynamic characteristics of an actual mechanical structure, it is necessary to identify the dynamic characteristics of the structure from the vibration response characteristics obtained by performing a vibration test on the structure. In this case, it is a parameter that directly expresses the vibration phenomenon, and it is easy to understand the phenomenon directly, and the method of identifying modal parameters (modal parameters such as natural frequency, mode damping ratio, and natural mode shape) with good convergence is currently used. It has become mainstream.

[0002]

2. Description of the Related Art FIG. 1 schematically shows the structure of a conventionally known vibration characteristic analyzing apparatus.
And an impulse hammer 2 for exciting the DUT 1 are prepared.

A load cell 3 for force detection is attached to the hammer 2, and a three-dimensional (x, y, z) vibration response of the hammer 2 due to the vibration of the hammer 2 is provided at an arbitrary position on the DUT 1.
An acceleration sensor 4 for detecting directional acceleration is attached.

Then, an output signal of the acceleration sensor 4 and an output signal of the load cell 3 are supplied to the arithmetic unit 5 to perform FFT processing to obtain frequency response function (hereinafter, may be abbreviated as FRF (compliance)) data. It has a configuration for calculating and outputting characteristics.

[0005] In this case, the vibration experiment is usually performed with the vibration place of the DUT 1 fixed by the hammer 2, and the acceleration sensor 4 is sequentially moved to give a plurality of response functions to the arithmetic unit 5. . The DUT 1 is fixed by a soft spring (not shown) or the like.

[0006] Many algorithms have been proposed as algorithms for the arithmetic unit (experimental mode analyzer) 5 of such a vibration characteristic analyzer. Among them, the "variation iterative method" in the frequency domain has a high accuracy. Currently widely used for its goodness.

However, the above-described improved iterative method has the following two disadvantages. Especially when there are a plurality of response functions, the parameter {γ}
Also increases as the response function increases, so that the calculation time of the inverse matrix in the equation for obtaining the correction vector {Δγ} for correcting the parameter {γ} sharply increases. Since the number of variables increases, iterative calculations are difficult to converge and diverge easily, and in some cases calculation becomes impossible.

The inventor of the present invention disclosed in Japanese Patent Laid-Open Publication No. Hei 8-5450, in order to improve the problem of the essence of the algorithm, which is easy to diverge and takes a long time to calculate, the arithmetic unit 5 compares the measured data with the measured data. The deviation iterative method in which the parameter is repeatedly changed until the residual from the theoretical value by the mode characteristic parameter is minimized is calculated separately for the linear term and the nonlinear term, and even if the linear term and the nonlinear term change, It has been proposed to identify the mode characteristics by providing a condition under which the residual does not diverge.

The improved variant iterative method will be briefly described below. The desired mode characteristic is represented by {γ}. The theoretical formula of the frequency response function FRF at the frequency point ω = ω _i is expressed by f
_i and the corresponding experimental data are denoted by y _i . When the number of FRFs to be referred to at the same time is represented by m and the number of frequency points is represented by N, the total number of data is 2Nm since the FRF is a complex number.

If the distribution of errors mixed in y _i is a normal distribution, the most likely mode characteristic is {γ} that minimizes the weighted residual sum of squares S of the following equation.

[0011]

(Equation 1)

Here, it is known that the weight coefficient W _i is proportional to the reciprocal (1 / σ ² ) of the population variance in the equation representing the probability density distribution of y _i .

All variables {γ} are divided into a linear term {α} (element u) and a nonlinear term {β} (element v). Here is an idea of eliminating {α}. A search is made for a region where the following equation is satisfied, using only {β} as a variable in the entire variable space.

[0014]

(Equation 2)

In the minimum region of {α} where the equation (2) is satisfied, the variable {α} of the theoretical formula f _i ({α}, {β}) of the FRF is represented by {β} as H _i. Then, the following equation is obtained.

[0016]

(Equation 3)

Here, the following equation is obtained from the equation (1).

[0018]

(Equation 4)

Once the differential coefficient of equation (3) is determined, the correction vector {Δβ} can be obtained by the least squares method using the following linear approximation.

[0020]

[C] ^T [W] [C] ｛Δβ｝ = [C] ^T [W] ｛y _i− f _{i i} (5) [C]: ∂H _i / ∂β _k the a ik element matrix [W]: W _i diagonal elements and diagonal matrix

The matrix [C] in equation (5) is called a Jacobian matrix. The size of this Jacobian matrix [C] is v × 2Nm, and increases with the number m of FRFs.

In recent years, with the development of a multi-channel data processing device and a multi-axis speed sensor, the arithmetic unit 4 has been able to obtain a large number of frequency response functions (FRF) in a short time by experiments.

For example, in FIG. 1, three FRFs are simultaneously obtained from information of three channels of X, Y, and Z from the three-axis acceleration sensor 4 and information of one channel of the load cell 3 for detecting a vibration input. Obtainable. When a plurality of acceleration sensors 4 are provided in parallel in order to shorten the measurement time, etc.,
The channel increases by the increase.

As shown in FIG. 16, this FRF has a frequency on the vertical axis and measurement points (point 1, point 2,...) On the horizontal axis.
Can be treated as a matrix of 2N rows and m columns in which the axis data X, Y, and Z are sequentially arranged. For example, a memory capacity of 4 bytes is allocated to one data.

However, the number m of the referred frequency response functions
With the increase in the number of pieces of information, it is becoming difficult to obtain quick results with a limited calculation speed and memory capacity of a computer or a conventional calculation method.

In order to solve this problem, since data to be processed simultaneously is given as an FRF matrix, this matrix may be reduced.

As a method of reducing the FRF matrix,
Conventionally, a singular value decomposition method has been proposed.

In this singular value decomposition method, the procedure for reducing the FRF matrix [P] (2N rows and m columns) obtained in the experiment using singular value decomposition will be described with reference to FIGS. 17 (1) to 17 (4). explain.

The matrix [B] (m rows and m columns) is calculated by the following equation (FIG. 1A).

[0030]

[B] = [P] ^T [P] Expression (6) [P] ^T : transposed matrix of matrix [P]

An eigenvalue λ and an eigenvector {x} of the matrix [B] are obtained (FIG. 2B). A set of the eigenvalue λ and the eigenvector {x} is arranged in descending order of the eigenvalue λ. Those whose eigenvalue λ is smaller than others are omitted. Although m eigenvalues of the matrix [B] are obtained, the number of modes n
Is reduced to a set of eigenvalues and eigenvectors {x _i } ((3) in the same figure).

From the n eigenvectors {x _i }, the FRF vector {P _i } reduced to n is calculated by the following equation ((4) in the figure).

[0033]

｛P _i ｝ = [P] ｛x _i ｉ i = 1１〜n Equation (7)

Therefore, in this method, it takes time to calculate the eigenvalues and eigenvectors of the equations (6) and (7). In order to calculate equations (6) and (7), matrix [P] (2N rows m
Column) and matrix [B] (m rows and m columns) need to be stored.

Also, since the eigenvalues are calculated, the eigenvalues and the eigenvectors are arranged in descending order, and those having small coefficients are truncated and compressed into n-mode data.
The physical meaning of the compression information was unknown.

Therefore, the present inventor has disclosed in Japanese Patent Application No. Hei 8-68251.
In the publication, a load cell for force detection attached to an impulse hammer for exciting a test object, and a vibration caused by the impulse hammer of the test object attached to an arbitrary place of the test object Receiving an acceleration sensor for measuring a response, an output signal of the load cell, and an output signal of the sensor to obtain frequency response function data; and a weighted residual square sum of the frequency response function data and a theoretical value of the frequency response function. A computing device that identifies the modal characteristics of the device under test by calculating the modal characteristic parameters by the least squares method that minimizes, in a vibration characteristic analysis device comprising: For the purpose of compression, the arithmetic unit may add +1 or-to the polarity of the imaginary part corresponding to the selected eigenmode of the frequency response function.
It has been proposed that the frequency response functions corresponding to the eigenmodes are respectively obtained by superimposing the frequency response functions which have been multiplied by 1 to make them positive, and the number of the frequency response functions is reduced to the number of the eigenmodes.

The physical meaning of this method will be described below. First, the FR in which the FRF of each point is arranged as a column vector
The F matrix [P] is reduced by the following equation.

[0038]

[Q] = [P] [K] Expression (8) [P]; FN matrix obtained by experiment 2N rows m
Column [Q]; Reduced FRF matrix 2N rows and n columns [K]; Matrix for reduction m rows and n columns n; Number of adopted eigenmodes (5 to 10) (m> n)

The method of creating the matrix [K] is as follows.
From (2N rows and m columns), only the imaginary part corresponding to the resonance peak of each eigenmode is taken out, and a matrix arranged in the row direction is created. Next, the sign of each component of the matrix is the same, and only the size is changed to 1 (see FIG. 2A).

A matrix [K] is created by transposing the matrix [K] ^T created in this manner (FIG. 2B).

Then, a reduced FRF matrix [Q] is created based on the FRF matrix [P] and the matrix [K] by equation (8) (FIG. 3 (3)).

In the calculation of equation (8), all FRFs are superimposed so that the sign of the imaginary part near the resonance peak of a certain eigenmode is the same.

This means that the direction of the force at each point is adjusted so that a certain eigenmode is most excited by setting the excitation point at the time of measuring the FRF matrix as a response point (measurement point). This corresponds to measuring the FRF when shaking. Therefore, the number of FRFs equal to the number of eigenmodes is replaced by the original FRF
Matrix can be reduced.

In practice, it is very difficult to correctly excite each point so that a certain mode is most excited. But,
If a large number of FRFs are obtained by a hit test or the like for each point, a measurement result corresponding to simultaneous impact excitation can be easily obtained using equation (8).

The processing of equation (8) applies -1 or + to each FRF.
Multiply by 1 and add. The absolute value of the addition coefficient is 1
Therefore, in the FRF after the addition, the random noise is canceled and reduced.

Even if the number of FRFs is large, the number of FRFs can be reduced to the same number as the number of eigenmodes to be extracted by the equation (8). The computation time and storage capacity for this are far less than methods using singular value decomposition or principal component analysis.
In addition, since all the FRFs are added, the random noise is canceled, and the variation in data is reduced.

[0047]

As described above, the inventor of the present invention has made it possible to efficiently perform multipoint reference for identifying a mode characteristic by using a large number of frequency response functions (FRFs) at the same time. Has been proposed to reduce the FRF while maintaining the following equation. In order to execute the maximum likelihood estimation method for identifying the most probable mode characteristic statistically, as shown in the above equation (1), , It is necessary to use the reciprocal of the variance as a weight function in curve fitting.

If the FRF is not reduced, the original FRF
Can be used as it is to obtain the variance for each FRF. However, when the input data obtained by reducing the FRF is used, the variance cannot be obtained by this method, so that the weight function cannot be obtained.

Therefore, in the above method, a value obtained empirically (for example, a value using the square of the reciprocal of the amplitude of the FRF) is used as the weighting function. Could not be guaranteed.

Accordingly, the present invention provides a load cell for force detection attached to an impulse hammer for exciting a DUT, and a load cell for force detection attached to an arbitrary position of the DUT. An acceleration sensor that measures the excitation response caused by the load cell, receives the output signal of the load cell and the output signal of the sensor, obtains frequency response function data, and calculates the frequency response function data and the theoretical value of the frequency response function. An arithmetic unit for identifying a mode characteristic of the device under test by calculating a mode characteristic parameter by a least square method that minimizes a weighted residual square sum, wherein the arithmetic unit selects the frequency response function. The polarity of the imaginary part corresponding to the eigenmode is multiplied by +1 or −1, respectively, and the frequency response function is set to a positive polarity. The purpose of the present invention is to theoretically derive a weighting function that can cope with the reduction of the frequency response function in a vibration characteristic analysis device that obtains corresponding frequency response functions and reduces the number of the frequency response functions to the number of the eigenmodes. I do.

[0051]

According to the present invention, a weighting function matrix is derived in accordance with statistical theory. In the maximum likelihood estimation method, it is necessary to use a reciprocal of a population variance of a response as a weight function. Since the population variance cannot be determined accurately in a finite number of experiments, the sample variance is used approximately. The real and imaginary parts of the averaged FRF obtained by averaging the FRFs obtained by q repeated experiments are
Assuming that the measurement point is j and the angular frequency is ω _i , it is expressed by the following equation.

[0052]

(Equation 9)

(Equation 10)

If no error is mixed in the input, the coherence function is given by the following equation.

[0054]

[Equation 11]

The real and imaginary parts of the sample variance are represented by the following equations.

[0056]

(Equation 12)

(Equation 13)

Ideally, the coherence function and the averaging F
It suffices if the real part and the imaginary part of the variance can be separately derived from the RF, but it is difficult with only the above information amount. Therefore, assuming that the real part and the imaginary part have the same weight for convenience, the following equations are derived from the above equations (9) to (13).

[0058]

[Equation 14]

Therefore, the following weight function is obtained as the reciprocal of equation (14).

[0060]

(Equation 15)

Next, the case of multi-point reference will be described. It is thought that there are many error factors in the FRF data obtained by the experiment. It is known that the errors follow a normal distribution, given that they are of similar magnitude and contribute independently of each other. According to statistical theory, the sample sum of m FRFs follows a normal distribution.

In the calculation of the above equation (8), since the imaginary parts of the points corresponding to the resonance peaks of the respective FRFs are added so as to have the same sign, the variance of the reduced FRF data is expressed by the equation (1)
The following equation is obtained from 4).

[0063]

(Equation 16)

Therefore, the weight function can be derived as in the following equation.

[0065]

[Equation 17]

Therefore, the arithmetic unit in the vibration characteristic analyzing apparatus according to the present invention further obtains a weighting function considering simultaneously a large number of the frequency response functions and applies the weighting function to the residual sum of squares. If used, the least squares method is identical to the maximum likelihood estimation method for estimating the most likely mode characteristics.

[0067]

DESCRIPTION OF THE PREFERRED EMBODIMENTS The configuration of a vibration characteristic analyzing apparatus according to the present invention can be the same as that of the conventional example shown in FIG. 1 described above. It is. Hereinafter, this calculation process will be described with reference to the flowchart shown in FIG.

<Partitioning of Parameters> When considering a plurality of response functions at the same time, if the number of response functions is m and the number of eigenmodes to be identified is n in the same manner as above, the parameter {γ}
Is 2n + (2n + 4) m. Among the parameters {γ}, the nonlinear term is 2n and the linear term is (2n + 4) m
It is. This parameter {γ} is divided into a linear term {α} and a nonlinear term {β} as in the following equation.

[0069]

１８γ｝ ^T = ｛｛α｝ ^T , ｛β｝ ^T ｛｛α｝ ^T = ｛‥, U _r , V _r , ‥, C, D, E, F, ‥｝ ｛β｝ ^T = ｛‥, σ _r , ω _dr , ‥｝ ・ ・ ・ Equation (18)

Here, U _r and V _r in the linear term {α} indicate the real part and the imaginary part of the r-th residue with respect to the measurement point.
C and 1 / D are the parameters introduced as the residual mass, and the effects of the eigenmodes, which are not noted, are used as the residual stiffness.

Σ _r in the nonlinear term {β} indicates each peak (n possible peaks), that is, the attenuation rate in each mode characteristic, and ω _dr is each peak, that is, the attenuation characteristic in each mode characteristic. Indicates angular frequency.

Since the response function (theoretical value) of the mode characteristic given to the arithmetic unit 5 when the DUT 1 is vibrated by the hammer 2 has an amplitude and a phase, the real part G _R is given by the following equation. A theoretical model represented by a complex number including (ω) and an imaginary part G ₁ (ω) is well known.

[0073]

[Equation 19]

The damping rate σ _r and the natural frequency ω _dr do not change even if the response point (the mounting position of the acceleration sensor 4) changes. That is, the number 2n of the parameter {β} of the nonlinear term is characterized in that the position (frequency) of the peak does not change even if the number m of the response functions increases, so that it is invariant. It turns out that it is not affected.

When the nonlinear term is obtained, only the linear term can be easily obtained.

Therefore, it is necessary to first obtain only the nonlinear term.
Husband. Therefore, first, the weighting coefficient W in the equation (1)_i
Is calculated by equation (15), and the nonlinear parameter
Damped natural angular frequency ω as {β}_drAnd mode decay rate σ
_rIs set (step S1 in FIG. 3).

Since the initial value of the damping converges at any value, only the initial value of the damping natural angular frequency may be obtained from the peak of the function obtained by summing the square of the FRF amplitude.

Then, a frequency response function y _i , which is experimental data, is input (step S 2).

Next, in order to reduce the following calculation load, the input frequency response function FRF is reduced using equation (8) (step S3). The reduction ratio is n (number of adopted modes) / m (number of FRFs).

When the nonlinear term parameter {β} is initialized, the linear term parameter {α} can be obtained by one calculation using the following equation (step S4).

[0081]

｛Α｝ = ([D] ^T [W] [D]) ⁻¹ [D] ^T [W] ｛y _i式 (20) where D _ij = ∂f _i / ∂ α _j (i = 1 to 2p, j = 1 to (2
n + 4) m)

Using the linear term {α} and the nonlinear term {β} obtained in this manner, the frequency response function F
The theoretical RF value f _i is calculated (step S5). The value of the first non-linear term is the value initially set in step S1.

The residual sum of squares S of the theoretical value and the experimental value of the response function thus calculated is calculated according to the equation (1) (step S6).

Then, it is determined whether or not this residual sum of squares S is smaller than the value obtained last time (step S
7). In this case, the initial residual sum of squares S may be set to an appropriate value.

As a result, when the residual sum of squares S is smaller than the previous value, it means that the residual sum of squares S has not yet reached the minimum value. The following non-linear term parameter {β} is obtained as described below.

Here, the above equation (5) is calculated by
Since it increases with the number m of F, it places a large load on both the storage capacity of the computer and the calculation time. To reduce such a load, the calculation method of the formula should be improved as follows. Can be.

As a method of solving the equation (5), there are (1) a method of solving as a simultaneous equation of v elements and v rows called a normal equation, and (2) an orthogonal transformation of the Jacobian matrix [C] by the Householder method or the like. The method of solving is known so far, and it is generally said that the latter (2) has better calculation accuracy.

When the calculations were actually performed, the results of the calculations in the case of double precision in this case were the same. In addition, the former (1), which solves a normal equation, is more advantageous in terms of calculation time and storage capacity.

Therefore, the present embodiment focuses on improving the calculation method using a normal equation. First, the normal equation of Expression (5) can be rewritten as the following expression.

[0090]

[G] ｛{β} = {d} (21)

[0091] expressions in (21) [G], when calculating the {d} from Jacobian matrix [C], a derivative of the theoretical equation H _i of the frequency response function, the residual sum of squares S of using Newton's method To the derivative of

By substituting equation (3) into equation (5), equation (2)
The right side of 1) is as follows.

[0093]

(Equation 22)

The following equation is obtained by developing and organizing the equation (22).

[0095]

(Equation 23)

Since {α} satisfies the minimum condition at a certain stage in the iterative calculation, {S / ∂α _s = 0
(S = 1 to u) holds. Therefore, the following equation is obtained.

[0097]

(Equation 24)

Next, the quantity obtained by omitting the second derivative term in the equation (4) is represented by the following equation.

[0099]

(Equation 25)

The component G _jk of the coefficient matrix [G] on the left side of Expression (21) is obtained from Expression (5) as follows.

[0101]

(Equation 26)

By substituting equation (3) into equation (26), the following equation is obtained.

[0103]

[Equation 27]

The right side of equation (27) = A + B + C, and by substituting equations (1) and (21), A, B, C
The following equations are obtained respectively.

[0105]

[Equation 28]

[0106]

(Equation 29)

[0107]

[Equation 30]

Here, {x} ^T , {a}, {b} are defined as follows.

[0109]

(Equation 31)

By substituting equation (31) into equation (30), the following equation is established.

[0111]

(Equation 32)

Therefore, equations (28), (29) and (32)
Is substituted into G _ik = A + B + C to obtain the following equation.

[0113]

[Equation 33]

Until now, the number of multiplications of the matrix [F] ⁻¹ of the second term on the right side of the equation (3) has been large, and this has been a burden on calculation time. Moreover, the calculation of the second term on the right side is v × 2Nm
Had been going around.

However, in the above equation (33), G
Since な has been removed to calculate _jk , this calculation is required only three times, and the whole [G] is 3v ² times.

The nonlinear term is a common term of all FRFs, and its number v is smaller than the normal frequency point N and the reference FRF number m, so that the calculation time can be reduced. Until now, the matrix [C] of 2Nm rows and v columns has been stored as it is, but [G] of v rows and v columns and {d} of v rows need only be stored, and the required storage capacity is reduced. I do.

That is, to calculate equation (33), S
It is only necessary to calculate the derivative of α and S with respect to α _s and β _k .

As described above, d _k is calculated using equation (24), and the differential coefficients of α and β of S and S ′ are calculated using equations (4) and (25). Then, _Gjk is calculated (step S8).

The correction vector {Δβ} is calculated by substituting the above result into the equation (21) (step S9).

<Setting of Side Constraint> The correction vector {Δβ} of the nonlinear term has been obtained as described above. However, if used as it is, divergence may occur. Therefore, a side constraint is further set.

As a side constraint of the divergence prevention, there is a method of reducing the same magnitude in the direction by applying a reduction factor to the correction vector. In this case, if the mode is different, the mutual effect of the nonlinear terms is reduced. The following constraint is set for each term.

[0122]

| △ ω _dr / ω _dr | ≦ K ₁ | △ σ _r / σ _r | ≦ K ₂ , σ _r > 0 Equation (34)

If {Δβ} does not meet this condition, it is corrected by equation (25) (step S11).

The correction vector {Δβ} obtained in this way and the previous nonlinear term parameter {β} (initial value are initially added) are added to calculate the new nonlinear term parameter {β ^* } this time. Step S12).

[0125]

３５β ^* ｝ = ｛β｝ + ｛Δβ｝ (35)

Using the new non-linear term parameter {β ^* } obtained in this way, the above-described step S4 and subsequent steps are repeatedly executed. If the residual sum of squares S calculated in step S6 of this iterative loop is not smaller than the previous value, it means that S has reached the minimum value,
The processing exits from the calculation loop of {β} (“N” in step S7).
o ").

Then, an FRF that does not shrink is input again (step S13), and the FRF and the already determined {β} are substituted into equation (19) to calculate a linear term {α} (step S14). ).

In this process, the linear least squares method can be applied, so that iterative calculation is unnecessary.

Then, {α} and {β} indicating the optimal mode characteristics are output (step S15).

[Verification by Model Data] Here, the vibration characteristic analyzing apparatus according to the present invention will be verified with model data as follows.

Table 1 shows the mass of each mass point and the value of the spring stiffness between each mass point in the 10-degree-of-freedom model having one end fixed in FIG. Assuming proportional viscous damping, damping coefficients c _i ,
It is assumed that the following relationship exists between the mass m _i and the spring stiffness k _i .

[0132]

[Table 1]

[0133]

C _i = αm _i + βk _i (α = 0.15, β = 3.0 × 10 ⁻⁶ ) Equation (36) The target frequency range and frequency resolution Δf are 0.1 to 6
40.1 (Hz) and Δf = 0.1 (Hz).

In this frequency domain, there are first to sixth order eigenmodes. The excitation point is mass point 1 and the response point is mass point 1.
From the simulation, 10 FRFs when mass point 10 is obtained are obtained by simulation and used as input data for multipoint reference.

FIGS. 5 to 10 show the results of identification of mass point 1 in the self-FRF. First, FIG. 5 (1) shows the phase, amplitude and coherence function of input data when identification is performed when noise is not included in the input (excitation force) and 4% normal distribution noise is mixed in the response. ) To (3).
6 (1) and 6 (2) show, using the input data of FIG. 5, the result of identification and the true value identified by, for example, the square of the reciprocal of the amplitude of FRF, which has been generally used as a weight function. Shown with phase. In FIG. 7 (1) and (2),
The result and the true value identified by the weight function of the equation (15) according to the present invention using the input data of FIG. 5 are shown together with the phase.

A comparison between FIG. 6 and FIG. 7 shows that when the input data contains no noise, the true value and the value of the identification result are almost the same in both methods, indicating that the true value can be identified.

Next, the excitation force and the response were respectively 2% and 4%.
Identification is performed by mixing noise according to the normal distribution of%. FIGS. 8A to 8C show the phase, amplitude, and coherence function of input data used for identification. FIGS. 9A and 9B show the result and the true value together with the phase when the identification is performed by the conventional weight function using the input data of FIG. FIGS. 10A and 10B show the result and the true value together with the phase when identification is performed by the weight function according to the present invention using the input data of FIG.

Comparing FIG. 9 and FIG. 10, when the identification is performed by the latter weighting function according to the present invention, the identification result and the true value are almost the same, whereas in the former case, the difference is large. Is recognized. Therefore, it can be seen that the identification accuracy of the present invention is clearly superior to that using the conventional weight function.

As described above, in the conventional identification using the weight function, the identification accuracy is good when the input does not include noise.
The identification accuracy decreases as the noise mixed into the input increases. On the other hand, in the identification by the weight function according to the present invention, since the weight function is obtained from the entire frequency response function, errors are accidentally canceled out, and the identification accuracy is good even if noise is mixed in the input.

[Identification Result of Simplified Bearing Beam] A vibration test is performed on a mild steel square bar shown in FIG. Identification is performed using the conventional weight function using the square of the reciprocal of the amplitude of the frequency response function and the weight function according to the present invention shown in Expression (15), and the two results are compared.

However, here, the third point is excited in the y-direction, and the first to fourteenth responses are adopted.

FIG. 12 shows the phase of the input data of the third self-response (FIG. 12 (1)) and the amplitude of the conventional method and the present invention (FIG. 12 (2)). As shown in the figure, it can be seen that the accuracy of the identification result according to the present invention is higher than that of the conventional method.

FIGS. 13 to 15 respectively show 1
The identification results of the mode shapes of the next eigenmode (first bending in the Z direction), the third eigenmode (third bending in the Z direction), and the fifth eigenmode (third bending in the Z direction) are shown. .

In FIGS. 13 and 14, the results of identification using both weighting functions can identify the true mode shape, but in FIG. 15, the results of the identification using the weighting function according to the present invention are better than the conventional weighting functions. Is close to a tertiary symmetric mode (center position is 391 mm), so it is considered to be close to the true value.

This is because the method of the present invention obtains a weighting function from all the FRFs and the coherence function, so that the influence of the noise included in the input data is canceled out, and a high-precision identification result is obtained in all the eigenmodes. Can be

[0146]

As described above, according to the vibration characteristic analyzing apparatus of the present invention, the least squares method for minimizing the weighted residual square sum of the frequency response function data and the theoretical value of the frequency response function is used. An arithmetic unit that identifies a mode characteristic of the DUT by calculating a mode characteristic parameter,
A frequency response function corresponding to the eigenmode is obtained by superimposing a positive frequency response function by multiplying the polarity of the imaginary part corresponding to the selected eigenmode of the frequency response function by +1 or −1, respectively. Reducing the number of the frequency response functions to the number of the eigenmodes,
Since the weighting function is obtained by simultaneously considering a large number of the frequency response functions and applied to the residual sum of squares, even if noise is mixed in the input, the identification accuracy is improved compared to the conventional weighting function. .

Further, since the weighting function according to the present invention is common to all the frequency response functions, the load on the storage capacity of the computer can be reduced when the input data increases in multipoint reference.

[Brief description of the drawings]

FIG. 1 is a block diagram showing a vibration characteristic analysis device common to the present invention and a conventional example.

FIG. 2 shows an F used in the vibration characteristic analyzer according to the present invention.
FIG. 3 is a diagram illustrating data reduction of an RF (frequency response function) matrix.

FIG. 3 is a flowchart showing an arithmetic algorithm stored and executed in an arithmetic unit in the vibration characteristic analyzing apparatus according to the present invention.

FIG. 4 is a diagram showing a 10-degree-of-freedom model having one end fixed.

FIG. 5 is a graph showing a phase, an amplitude, and a coherence function of input FRF data used for identification (No. 1: a signal obtained by mixing 4% normal distribution noise into a response without input noise). .

6 is a graph showing a phase and an amplitude when the input data shown in FIG. 5 (part 1) is identified using a conventional weight function.

FIG. 7 is a graph showing a phase and an amplitude when the input data (part 1) shown in FIG. 5 is identified by using the weight function of the present invention.

FIG. 8 is a graph showing the phase, amplitude, and coherence function of input FRF data used for identification (part 2: input and response mixed with 2% and 4% of normally distributed noise, respectively); is there.

9 is a graph showing a phase and an amplitude when the input data (No. 2) shown in FIG. 8 is identified using a conventional weight function.

FIG. 10 is a graph showing a phase and an amplitude when the input data (part 2) shown in FIG. 8 is identified using the weight function of the present invention.

FIG. 11 is a view showing a mild steel square bar to be subjected to a vibration test simulating a bearing beam.

12 is a graph showing the phase and amplitude of the input data of the self-response at the third point on the square bar in FIG. 11;

13 is a graph showing a result of identifying a mode shape of a first-order eigenmode (bending in the Y direction) when the square bar of FIG. 11 is excited.

14 is a graph showing a result of identifying a mode shape of a second-order eigenmode (bending in the Y direction) when the square bar of FIG. 11 is vibrated.

FIG. 15 is a graph showing a result of identifying a mode shape of a third-order eigenmode (third bending in the Y direction) when the square bar of FIG. 11 is vibrated.

FIG. 16 is a diagram showing an example of an FRF matrix.

FIG. 17 is a diagram showing an example of data reduction by a conventional singular decomposition method.

[Explanation of symbols]

 DESCRIPTION OF SYMBOLS 1 DUT 2 Impulse hammer 3 Load cell 4 Acceleration sensor 5 Arithmetic unit In a figure, the same numerals show the same or corresponding parts.

Claims (1)

[Claims] 1. A load cell for force detection attached to an impulse hammer for vibrating a test object, and a load cell attached to an arbitrary position of the test object and caused by the impulse hammer of the test object. An acceleration sensor for measuring vibration response, an output signal of the load cell, and an output signal of the sensor to obtain frequency response function data, and a weighted residual between the frequency response function data and a theoretical value of the frequency response function. An arithmetic unit for identifying a mode characteristic of the device under test by calculating a mode characteristic parameter by a least square method of minimizing a sum of squares, wherein the arithmetic unit operates on a selected eigenmode of the frequency response function. A frequency response function corresponding to the eigenmode is obtained by superimposing a positive frequency response function by multiplying the polarity of the corresponding imaginary part by +1 or -1, respectively. A vibration characteristic analysis device for calculating each of the frequency response functions and reducing the number of the frequency response functions to the number of the eigenmodes; A vibration characteristic analysis device characterized by being applied to summation.