JP3335765B2 - Vibration characteristic analyzer - Google Patents

Description

DETAILED DESCRIPTION OF THE INVENTION

[0001]

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 modal characteristic of the test object from its 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.

[0003]

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.

[0004] A load cell 3 for force detection is attached to the hammer 2, and a three-dimensional (x, y, z) acceleration of a vibrating force by the hammer 2 is detected at an arbitrary position on the DUT 1. Pickup 4 is installed.

[0005] The output signal of the pickup 4 and the output signal of the load cell 3 are supplied to an arithmetic unit 5 and subjected to FFT processing to perform a frequency response function (transfer function: compliance).
And calculates and outputs the mode characteristics.

In this case, the vibration experiment is usually performed with the vibration place of the DUT 1 fixed by the hammer 2, and the pickup 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.

[0007] 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 the most widely used for its goodness.

[0008] In brief, this partial deviation iterative method is to perform an iterative calculation while testing an error between experimental data and a theoretical value, and identify a most probable mode characteristic by a least squares method. Furthermore, weighting can be freely performed in the frequency domain, and the influence of mode characteristics other than the mode characteristics of interest can be considered.

[0009] To explain <deviance iteration method> below the deviance iterative method, first, the response function of the modal properties imparted from the pickup 4 when the cuff DUT 1 by a hammer 2 to the calculating unit 5 ( Since the theoretical value has an amplitude and a phase, it is well known that it is represented by a complex number including a real part G _R (ω) and an imaginary part G _I (ω) as in the following equation.

[0010]

(Equation 1)

In this case, the angular frequency is represented by ω, and the number of eigenmodes (degrees of freedom) within the angular frequency range of interest is n (corresponding to the number of peaks indicating resonance points in FIG. 7 described later). And
Note that U _r and V _r indicate the real part and the imaginary part of the r-th residue with respect to the measurement point.

The sum of squares of the error between the complex number and the experimental data obtained from the pickup 4 (data with the frequency after the FFT calculation processing as the horizontal axis and the response function (compliance = displacement / force) as the vertical axis) ( The parameters (constants) in equation (1) are determined so that the residual sum of squares is minimized, and the mode characteristics are determined from the parameters.

That is, from the experimental data corresponding to the equation (1), the most likely parameters {γ} ^T = ｛‥, σ _r , ω _dr , U _r ,
To determine V _r , ‥, C, D, E, F, ‥｝ (｛‥｝ represents a vector), the residual square sum S is calculated assuming that the error of the experimental data with respect to the theoretical value is normally distributed. If the parameter {γ} is sequentially corrected from the initial value (this is a value obtained by reading the frequency of the resonance point of the experimental data E in FIG. 7) so as to be minimum, the experimental data E (solid line) as shown in FIG. ) Theoretical value T (likelihood) closest to (likelihood) frequency response characteristic
Is obtained. This residual sum of squares S is expressed by the following equation.

[0014]

(Equation 2)

When a plurality of response functions are processed simultaneously, their real and imaginary parts may be arranged in series, and the number of response functions (the number of measurement points) is m and one response function ( Assuming that the number of data during one excitation) is N, p = mN
Becomes _Wi is a weight and is proportional to the reciprocal of the variance of each measured data.

In order to obtain the parameter {γ} that minimizes the residual sum of squares S, the denominator of the response function (theoretical expression) f in the equation (2) in the parameter {γ} includes σ _r and ω _dr . Since a nonlinear term is included, a modified vector {γ} for correcting the parameter {γ} is obtained by applying a simplified improved iterative method (Gauss-Newton method) based on linear approximation as follows: The calculation may be repeatedly performed until the sum of squared differences S converges to the minimum value.

[0017]

｛△ γ｝ = ([A] ^T [W] [A]) ⁻¹ · [A] ^T [W] ｛y _i −f _i・ ・ ・ (3) where A _ij = ∂f _i / ∂γ _j , which is the theoretical formula
(1) shows how much the theoretical equation (1) changes when one of γ in (1) is changed, which forcibly converts the denominator including the nonlinear term into the numerator. "Approximation". Also it shows the [W] is a diagonal matrix with diagonal elements of W _i.

[0018]

However, the above-described partial deviation iteration 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 of equation (3) rapidly increases (see characteristic A in FIG. 3). Since the number of variables increases, iterative calculation is difficult to converge and diverges easily, and the calculation may not be possible on the way (see FIG. 8).

Accordingly, the present invention provides a load cell for force detection attached to an impulse hammer for exciting a test object, and a three-dimensional acceleration of the exciting force attached to an arbitrary position of the test object. And an arithmetic unit that receives the output signal of the pickup and the output signal of the load cell, obtains a frequency response function, calculates a mode characteristic parameter, and identifies a mode characteristic of the device under test. It is an object of the present invention to provide an improved vibration characteristic analyzing apparatus provided with the essential nature of an algorithm that is easy to diverge and takes a long time to calculate.

[0020]

In order to achieve the above object, a vibration characteristic analyzing apparatus according to the present invention is arranged such that an arithmetic unit operates until the residual between the measured data and the theoretical value based on the mode characteristic parameter is minimized. The mode characteristic is identified by calculating the deviation iterative method in which parameters are repeatedly changed by dividing into a linear term and a non-linear term, and providing a condition that the residual does not diverge even when the linear term and the non-linear term change. It is characterized by doing.

In the arithmetic unit, the condition that the divergence does not occur is that the partial derivative of the residual by the linear term and the nonlinear term becomes zero, and the second-order partial differential term at this time can be omitted. .

Further, it is preferable that the arithmetic unit adds a reduction condition within a threshold value to a correction vector of the nonlinear term when the deviation iterative method is executed, as the condition under which the divergence does not occur.

[0023]

In the vibration characteristic analyzing apparatus according to the present invention, a mode characteristic parameter is divided into a linear term and a non-linear term.

That is, the number of nonlinear term parameters is invariant even if the number of response functions increases, because the position (frequency) of the peak shown in FIG. 7 does not change, and is not affected by the number of response functions. This is a feature, and if this nonlinear term is found, a linear term can be easily found.

Therefore, only for the nonlinear term, first, an initial value (temporary value) of the nonlinear parameter is set, and the linear term parameter is calculated by this initial value.

The response functions G _R and G _I are calculated by the equation (1) using the linear term thus obtained and the nonlinear term of the initial value.

The response functions G _R , G thus calculated
The residual sum of squares S of _I and the experimental value is calculated according to equation (2),
When the residual sum of squares S is found to be the minimum value (minimum value), it is determined that the linear term and the nonlinear term obtained as described above are the final values of the mode characteristics to be identified.

When it is determined that the residual sum of squares S is not the minimum value, the residual sum of squares S finds a new nonlinear term in order to minimize the residual sum of squares S.

In order to obtain only the nonlinear term, a correction vector {β} is obtained in the same manner as in the above equation (3), and iterative calculation is performed until convergence. A constraint condition for preventing divergence is provided near the minimum point of the sum of squared differences S, and the correction vector {β} of the nonlinear term is obtained by the constraint condition.

Note that the correction vector {β} is obtained,
If used as it is, divergence may occur, so a side constraint is further set.

The correction vector {β} thus obtained and the previous nonlinear term parameter {β} (initial value are initially added) are added to calculate the new nonlinear term parameter, and this new nonlinear term By repeatedly executing using the parameters, it is possible to determine {α} and {β} showing the optimal mode characteristics.

[0032]

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS The configuration of the 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, but the difference from the conventional example is the arithmetic processing in the arithmetic unit 5. Hereinafter, this calculation process will be described with reference to the flowchart shown in FIG.

(1) Parameter division: When considering a plurality of response functions simultaneously, if the number of the response functions is m and the number of eigenmodes to be identified is n as above, the number of parameters {γ} Is 2n + (2n + 4) m. The parameter {γ}
In Equation (1), the nonlinear term in the denominator is 2n, and the linear term in the numerator is (2n + 4) m. This parameter {γ} is divided into a linear term {α} and a nonlinear term {β} as in the following equation.

[0034]

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

Here, σ _r in the nonlinear term {β} represents each peak (n possible peaks shown in FIG. 7), that is, the attenuation rate in each mode characteristic, and ω _dr is each peak, that is, each mode. The damped natural angular frequency in the characteristic is shown. This damping rate and natural frequency are the response points (the mounting position of the pickup 4).
Does not change even if changes. That is, the number 2n of the parameter {β} of the non-linear term is characterized in that the position (frequency) of the peak shown in FIG. It turns out that it is not affected by the number of functions. Then, if the nonlinear term is obtained, only the linear term can be easily obtained.

Therefore, it is devised to obtain only the nonlinear term first. For this reason, first, an initial value of the nonlinear parameter {β} is set (step S1 in FIG. 2). This is performed by reading the peak frequency of the frequency response as shown in FIG. 7 and setting the value to ω _dr, and setting σ _r to an appropriate value.

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

[0038]

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

Here, as shown in equation (4), {α} ^T
= ｛‥, U _r , V _r , ‥, C, D, E, F, ‥｝, so that D _ij
Is as follows.

First, when i = 1 to p, f _i = G
It is the _{R (ω i), D ij} = ∂G R (ω i) / ∂α j , and the following equation is obtained for each linear terms.

[0041]

(Equation 6)

Further, a _{r and} b _r in the above equation (6) are also _{expressed by the following} equations.
As shown in (1), it can be expressed by nonlinear terms σ _r, ω _dr ,
Since the nonlinear terms σ _{r and} ω _dr are initialized as described above, D _ij can be calculated by inserting their initial values and ω _i .

Similarly, when i = p + 1 to 2p, f
_{Since i} = G _I (ω _i ), D _ij = ∂G _I (ω _i ) / ∂α _j , and the following equation is obtained for each linear term.

[0044]

(Equation 7)

Also in this equation, D _ij can be calculated by inserting the initial value of the nonlinear term and ω _i .

Using the linear term {α} and the nonlinear term {β} obtained in this manner, the response functions G _R and G
_I is calculated (step S3).

The response functions G _R , G thus calculated
The residual sum of squares S of _I and the experimental value is calculated according to equation (2) (step S4).

Then, it is determined whether or not the residual sum of squares S is smaller than the value obtained last time (step S).
5). 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 not smaller than the previous value, it means that the residual sum of squares S has reached the minimum value. The nonlinear term {β} is output as the final value of the mode characteristic to be identified (step S6), and this routine ends.

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. Is obtained as follows.

That is, in order to obtain only the nonlinear term, the correction vector {β} is calculated by the following equation using {γ} in the above equation (3).
And iterative calculation is performed until convergence.

[0052]

｛△ β｝ = ([C] ^T [W] [C]) ⁻¹ [C] ^T [W] · {y _i −f _i } (8) where C _ij = ∂f _i / ∂β _j (i = 1 to 2p, j = 1 to 2
n)

However, as it is, the linear term {α}
Since this is not taken into account, the following is devised in order to consider the linear term {α}.

The following equation is satisfied at the minimum point of the residual sum of squares S.

[0055]

９S / ∂α _j = 0 (j = 1 to (2n + 4) m) Expression (9)

If the following equation is satisfied in this vicinity, divergence can be prevented.

[0057]

１０ ² S / ∂α _j ∂β _k = 0 (j = 1 to (2n + 4) m, k = 1 to 2n) Expression (10)

That is, the meaning of the equation (10) is expressed as ∂S / ∂α _j
Shows how much the residual square sum S changes when the linear term is changed, and shows that ∂S / ∂α _j does not change even if the nonlinear term is changed by partially differentiating this with ∂β _k. (Because the right side is 0), this is a constraint for preventing divergence in consideration of the linear term {α}.

The number in the above equation (10) is 2n × (2n + 4) m, which is too large to eliminate the linear term {α}.
Since the objective is to reduce the variables rather than satisfying all the conditions of (0), if the derivative of the response function f by β _k under the equation (10) is _denoted by H _k , the derivative H _k becomes A calculation represented by an equation is performed (step S7).

[0060]

[Equation 11]

[0061] In order to determine the differential coefficient H _k is, ∂α
_It is necessary to calculate _j / ∂β _k (j = 1 to (2n + 4) m, k = 1 to 2n).

Therefore, the following equation is obtained by substituting the residual sum of squares S of equation (2) into equation (10).

[0063]

(Equation 12)

Substituting H _k of equation (11) for ∂f _i / ∂β _{k in} the first term on the left side of equation (12), and further rearranging and displaying the matrix, the following equation is obtained.

[0065]

(Equation 13)

Therefore, by transforming equation (13), {a} =
By setting ｛b｝ [B] ⁻¹ , ∂α _j / ∂β _k can be obtained.

[B] is a square matrix of (2n + 4) m × (2n + 4) m. When s and t belong to different response functions in _Bst , [B] becomes As shown, since B _st = 0, it can be divided into m (2n + 4) × (2n + 4) matrices.

[0068]

[Equation 14]

Therefore, the magnitude of the inverse matrix [B] ^-1 for obtaining {a} is independent of the number m of the response functions.

Then, if the differential coefficient H _k of the equation (11) is obtained, the correction vector 非 線形
β｝ is obtained (step S8).

[0071]

｛△ β｝ = ([C] ^T [W] [C]) ⁻¹ [C] ^T [W] · {y _i −f _i } (15) where C _ij = H _j (f = f _i ) (i = 1 to 2p, j = 1 to 2)
n)

<Setting of Side Constraint> Although the correction vector {β} of the nonlinear term is obtained by the equation (15), it may diverge if used as it is despite the constraint of the above equation (10). Further, side constraints are set.

As a side constraint of the divergence prevention, there is a method of reducing the size by applying a reduction factor to the correction vector and keeping the direction the same, but in this case, if the mode is different, the influence of the non-linear terms is reduced. Therefore, the following constraint is set for each nonlinear term (step S9).

[0074]

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

In this case, the coefficient K ₁ is 0.05 and the coefficient K ₂ is 0.8. If K ₂ <1 and the initial value of σ _r is positive, σ
_r remains positive after repeated calculations.

The correction vector {β} thus obtained is added to the previous nonlinear term parameter {β} (the initial value is an initial value) to calculate a new nonlinear term parameter {β ^* } this time ( Step S10).

[0077]

１７β ^* ｝ = ｛β｝ + ｛Δβ｝ (17)

Using the new non-linear term parameter {β ^* } obtained in this way, the above-mentioned step S3 and subsequent steps are repeatedly executed, and {α} and {β} indicating the optimal mode characteristics are output (step S3). S6).

In the present invention, the second expansion of the response function f, which is the second item of the equation (12), can be similarly ignored and expanded.

[0080]

As described above, according to the vibration characteristic analyzing apparatus of the present invention, the deviation is repeatedly changed until the residual between the measured data and the theoretical value based on the mode characteristic parameter is minimized. The mode characteristic is identified by calculating the iterative method separately for a linear term and a non-linear term and providing a condition that the residual does not diverge even if the linear term and the non-linear term change. As shown in FIG. (1) Calculation time: A comparison of the calculation time is shown in FIG. 3 and Table 1 below. The number n of modes to be identified is 5, and the number of data in one compliance is 573. In single point reference, the computer use time is shorter in the conventional example,
In the conventional example, when m of the frequency response function increases, the computer use time rapidly increases. This is because the size of the inverse matrix to be solved is (2n + 4) m. In the improved type, the size of the inverse matrix is constant at 2n and increases almost linearly because the number of calculations using the inverse matrix only increases.

[0081]

[Table 1]

(2) Convergence: The compliance created by calculation is used to compare the convergence. In the model, the end of the first spring is fixed with five mass points and five springs. Table 2 shows the specifications and the undamped natural frequency. The number of mode characteristics to be identified is 5, and the data is 901 at every 1 Hz from 50 to 950 Hz. The modal attenuation ratio is 1 for all modal characteristics.
%.

[0083]

[Table 2]

In order to compare the convergence between the present invention and the conventional example, an iterative calculation was performed by changing the initial value of the nonlinear term parameter ω _dr . The initial value of the modal damping ratio is the correct value of 1.
/ 10 is 0.1%.

The results are shown in FIGS. 4 to 6.
The vertical axis is the normalized residual Z, which is represented by the following equation.

[0086]

(Equation 18)

FIG. 4 corresponds to the case where the second derivative term of the second term on the left side of the equation (12) is omitted, and FIG. 5 corresponds to the case where the second derivative term is not omitted. FIG. 6 shows the characteristics in the case where the second-order differential term is not omitted and the side constraint by the equation (16) is taken into consideration. The convergence is improved in this order, and the case of FIG. It can be seen from FIG.

[Brief description of the drawings]

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

FIG. 2 is a flowchart showing a calculation algorithm stored and executed in a calculation device in the vibration characteristic analysis device according to the present invention.

FIG. 3 is a graph for comparing the operation time of the present invention with that of a conventional example.

FIG. 4 is a graph showing convergence (1) of the vibration characteristic analyzing apparatus according to the present invention.

FIG. 5 is a graph showing convergence (No. 2) of the vibration characteristic analyzing apparatus according to the present invention.

FIG. 6 is a graph showing convergence (3) of the vibration characteristic analyzing apparatus according to the present invention.

FIG. 7 is a graph showing a measurement result and a theoretical value of compliance (response function) when a test object is vibrated with respect to a frequency.

FIG. 8 is a graph showing convergence of a conventional example.

[Explanation of symbols]

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

──────────────────────────────────────────────────続 き Continued on the front page (58) Field surveyed (Int.Cl. ⁷ , DB name) G01H 13/00 G01M 7/08 JICST file (JOIS)

Claims (3)

(57) [Claims] 1. A load cell for force detection attached to an impulse hammer for exciting an object to be tested, and a three-dimensional acceleration of the exciting force to be installed at an arbitrary position on the object to be tested. And a calculation device that receives the measured data from the pickup and the output signal of the load cell, determines a frequency response function, calculates a mode characteristic parameter, and identifies the mode characteristic of the device under test. In the vibration characteristic analyzing apparatus, the arithmetic unit divides a deviation iterative method in which the parameter is repeatedly changed until a residual between the measured data and a theoretical value based on the parameter is minimized into a linear term and a nonlinear term. A vibration characteristic characterized by identifying the mode characteristic by providing a condition in which the residual does not diverge even if the linear and nonlinear terms change. Gender analysis device. 2. The arithmetic unit according to claim 1, wherein the non-divergence condition is that partial differentiation of the residual by the linear and nonlinear terms becomes zero, and the second-order partial differential term at this time is omitted. The vibration characteristic analysis device according to claim 1. 3. The method according to claim 1, wherein the arithmetic unit adds a reduction condition within a threshold value to a correction vector of the nonlinear term when executing the partial repetition method as the non-divergence condition. Or the vibration characteristic analyzer according to 2.