JP2009204517A - Oscillation analytic method and oscillation analytic equipment - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide an oscillation analytic method or oscillation analytic equipment capable of identifying features of a rigid body while hardly deteriorating identification accuracy even if measuring points to be actually measured are decreased. <P>SOLUTION: The degree of freedom of a characteristic matrix to be identified is made to be larger than the measuring degree of freedom to be actually measured at measuring points in an oscillating test, and thus, the number of values of mode characteristic computed based on the measuring data is made to be smaller than the number of values of mode characteristic computed based on the characteristic matrix. When identifying the characteristic matrix, the characteristic matrix is identified by altering values of the component of the characteristic matrix so as to coincide with the value of the mode characteristic computed based on the characteristic matrix to the value of the mode characteristic computed based on measured data only for a mode characteristic corresponding to the mode characteristic computed based on the measuring data among the mode characteristics computed based on the characteristic matrix. <P>COPYRIGHT: (C)2009,JPO&INPIT

Description

  The present invention relates to a vibration analysis method and a vibration analysis apparatus.

  In complex structures such as automobiles and engines, the rigid body characteristics (mass, center of mass position, three main moments of inertia, and directions of the three main shafts corresponding to these) are the motion performance and vibration noise performance of the structure. It is one of the largest factors that affect the performance of vibration and noise control. In fact, the value of the rigid body characteristics is the most necessary parameter when starting analysis, design, and optimization of various dynamic behaviors such as motion analysis, vibration analysis, vibration isolation mechanism analysis, control performance analysis, etc. for design support. It is. Therefore, it is very important to easily identify (estimate) the rigid body characteristics of a complex structure with practical accuracy.

  When identifying the rigid body characteristics of a structure such as an automobile, it is difficult to apply the classic rigid body characteristic measurement method by suspending the structure with a wire. Therefore, a pendulum bench type measuring device has been developed and used at present. However, in this measuring device, since it is necessary to handle the pendulum bench portion on which the automobile is mounted and the automobile to be measured as a rigid body, setting takes time. In addition, when measuring the moment of inertia, it is necessary to identify the position of the center of mass of the automobile in advance and install the automobile so that the center of mass is accurately located on the rotation center axis of the measuring device. As described above, the pendulum bench type measuring device has several problems in measuring the rigid body characteristics.

  The inventor of the present application has already proposed a new method named “Experimental characteristic matrix identification method” as a solution to the above problem (Patent Document 1). In this method, it is not necessary to treat the object to be identified as a rigid body. Rather, the rigid body characteristics can be obtained from an appropriate number of resonance vibration characteristics from the first to the second or third order as an elastic body. In addition, all rigid body characteristics can be identified by performing a single-point excitation multipoint response test once as an excitation test.

JP 2001-350741 A

  By the way, in the method proposed by the inventor, a single-point excitation multipoint response test is performed as an excitation test. Here, in order to identify the rigid body characteristic with high accuracy for a complex structure, it is necessary to set a certain number of measurement points depending on the complexity and size.

  However, when the number of measurement points for actual measurement increases, the time required for the measurement work of the vibration test including calibration of the sensor for measurement and setting of the measurement system increases. In addition, since the frequency response function is identified for each measurement point, the analysis load by the computer increases as the number of measurement points actually measured increases. On the other hand, when the vibration test is performed with the number of measurement points reduced, the identification accuracy of the rigid body characteristics is lowered.

  Therefore, in view of the above problems, an object of the present invention is to provide a vibration analysis method or vibration analysis apparatus capable of identifying rigid body characteristics without substantially reducing identification accuracy even when the number of measurement points for actual measurement is reduced. Is to provide.

  In order to solve the above-mentioned problem, in the first invention, a step of requesting coordinate data of a plurality of measurement points set on the analysis object, and an excitation test for exciting at least one point on the analysis object A step of calculating a mode characteristic based on measurement data at each measurement point obtained at the time of the measurement, and a constraint condition between components constituting a characteristic matrix representing the vibration of the analysis object based on the coordinate data In the vibration analysis method, comprising the steps of: calculating an expression; and identifying the characteristic matrix using the constraint condition expression based on the value of the mode characteristic calculated based on the measurement data. The degree of freedom of the characteristic matrix is greater than the degree of freedom of measurement actually measured at the measurement point in the test, and therefore the number of mode characteristic values calculated based on the measurement data is greater than the characteristic. Less than the number of mode characteristic values calculated based on the matrix, and in the step of identifying the characteristic matrix, the mode characteristics calculated based on the measurement data among the mode characteristics calculated based on the characteristic matrix For only the mode characteristics corresponding to, change the value of each component of the characteristic matrix so that the value of the mode characteristic calculated based on the characteristic matrix matches the value of the mode characteristic calculated based on the measurement data. By doing so, the characteristic matrix is identified.

  In order to solve the above problem, in the second invention, a step of requesting coordinate data of a plurality of actual measurement points and virtual measurement points set on the analysis object, and excitation at at least one point on the analysis object Performing a measurement only at actual measurement points when the vibration test is performed, calculating a mode characteristic based on measurement data at each actual measurement point obtained by vibration, and the analysis object based on the coordinate data A step of calculating a constraint equation between components constituting the characteristic matrix representing the vibration of the vibration, and identifying the characteristic matrix using the constraint equation based on the calculated mode characteristic value based on the measurement data The characteristic matrix is created with degrees of freedom equal to the degrees of freedom at both the actual measurement point and the virtual measurement point, and the characteristic matrix is identified. The mode characteristic value calculated based on the characteristic matrix is only measured for the mode characteristic corresponding to the mode characteristic calculated based on the measurement data among the mode characteristics calculated based on the characteristic matrix. The characteristic matrix is identified by changing the value of each component of the characteristic matrix so as to match the value of the mode characteristic calculated based on the data.

  In a third invention, in the second invention, the virtual measurement point is provided such that a measurement point closest to the virtual measurement point is an actual measurement point.

  According to a fourth invention, in the second or third invention, the method further comprises a step of requesting an input of a structural coupling state between measurement points. In the step of calculating the constraint condition expression, in addition to the coordinate data, A constraint equation is calculated based on the structural coupling state between the measurement points.

  In a fifth invention, in any one of the second to fourth inventions, the characteristic matrix identified in the step of identifying the characteristic matrix is a mass matrix, and the mass identified in the step of identifying the characteristic matrix The method further includes calculating a rigid body characteristic based on the matrix.

  In order to solve the above problems, in the sixth aspect of the invention, a measuring device that measures the displacement of the position, velocity, or acceleration at the measurement point on the analysis object, and an excitation that performs excitation at at least one point on the analysis object. In the vibration analysis apparatus comprising a vibrator and an analysis unit that receives measurement data and vibration data from the measurement device and the vibrator, respectively, and performs arithmetic processing based on the data, the analysis unit includes the above A component that constitutes a characteristic matrix representing the vibration of the analysis object based on the coordinate data of each measurement point based on the measurement data obtained at the measurement point obtained by the vibration test with the vibrator The constraint matrix is calculated, and the characteristic matrix is identified using the constraint formula based on the value of the mode characteristic calculated based on the measurement data. A mode in which the degree of freedom of the characteristic matrix is greater than the degree of freedom of measurement actually measured, and thus the number of mode characteristic values calculated based on the measurement data is calculated based on the characteristic matrix When identifying the characteristic matrix, which is smaller than the number of characteristic values, only the mode characteristic corresponding to the mode characteristic calculated based on the measurement data among the mode characteristics calculated based on the characteristic matrix. The value of each component of the characteristic matrix is changed so that the value of the mode characteristic calculated based on the characteristic matrix matches the value of the mode characteristic calculated based on the measurement data. Identify the matrix.

  In order to solve the above-described problem, in the seventh invention, a measuring instrument that measures a displacement of a position, velocity, or acceleration at an actual measurement point among actual measurement points and virtual measurement points set on the analysis object, and the analysis object A vibration exciter that excites at least one point on the object, and an analysis unit that receives measurement data and vibration data from the measuring instrument and the vibration exciter, and performs arithmetic processing based on the data, respectively. In the vibration analysis apparatus, the analysis unit calculates a mode characteristic based on measurement data at an actual measurement point obtained by an excitation test using the vibrator, and is based on the coordinate data of the actual measurement point and the virtual measurement point. Calculate the constraint equation between the components that make up the matrix, and identify the characteristic matrix using the constraint equation based on the value of the mode characteristic calculated based on the measurement data The characteristic matrix is created with a degree of freedom equal to the degrees of freedom at both the real measurement point and the virtual measurement point, and when identifying the characteristic matrix, among the mode characteristics calculated based on the characteristic matrix, Only for the mode characteristic corresponding to the mode characteristic calculated based on the measurement data, the value of the mode characteristic calculated based on the characteristic matrix matches the value of the mode characteristic calculated based on the measurement data. The characteristic matrix is identified by changing the value of each component of the characteristic matrix.

  In an eighth aspect based on the seventh aspect, the virtual measurement point is provided such that a measurement point closest to the virtual measurement point is an actual measurement point.

  In a ninth invention, in the seventh or eighth invention, when calculating the constraint condition expression, the constraint condition expression is calculated based on a structural coupling state between measurement points in addition to the coordinate data. .

  In a tenth invention, in any one of the seventh to ninth inventions, the characteristic matrix identified is a mass matrix, and the analysis unit applies the mass matrix identified in the step of identifying the characteristic matrix. Based on this, the rigid body characteristics are calculated.

  According to the present invention, even if the number of measurement points at which measurement is actually performed is reduced, the rigid body characteristics can be estimated without substantially reducing the estimation accuracy.

  Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings.

  FIG. 1 is a block diagram of the analysis apparatus of this embodiment. In FIG. 1, the analysis device 10 is a device that can analytically calculate the rigid body characteristics and the like of the analysis object 1, and includes a measurement unit 20 and an analysis unit 30. The measurement unit 20 includes a hammer 20 having a force transducer 21 at the tip, a plurality of measuring instruments (for example, accelerometers) 23a to 23c attached to the analysis target 1, a force transducer 21 and each measuring instrument 23a. Amplifiers 24 and 25a to 25c connected to .about.23c, and an external interface unit 26 connected to these amplifiers 24 and 25a to 25c.

  On the other hand, the analysis unit 30 accepts input from a user, a removable media device (RM device) 31 that writes to removable media, a removable media control unit (RM control unit) 32 that controls the removable media device 31, and the like. An input device 33, an input control unit 34 for controlling the input device 33, an output device 35 for outputting to a display, an output control unit 36 for controlling the output device 35, and a central processing unit 37 for performing various arithmetic processes And a storage device 38 and a bus 39 for interconnecting these external interface unit 26, removable media control unit 32, input control unit 34, output control unit 36, central processing unit 37, storage device 38 and the like.

  The striking hammer 22 is used to vibrate the analysis object 1. The magnitude of the force applied at the time of vibration is converted into an electric signal by a force converter 21 provided on the hammering hammer 22. This electric signal is amplified by the amplifier 24 and then output to the external interface 26. In this embodiment, an unsteady wave is given to the analysis object 1 using the hammering hammer 22 as a vibrator, but the vibrator and the vibration method are not limited to this. As long as the analysis object 1 can be vibrated and the magnitude of the force applied to the analysis object 1 is known, any vibration exciter may be used instead of the striking hammer 22. Accordingly, a hydraulic or piezoelectric exciter, a non-contact exciter using sound pressure or a magnetic field, or the like can be used as the exciter. Further, instead of the non-stationary wave, various excitation waveforms such as a stationary wave and a periodic wave may be excited. Furthermore, in this embodiment, a single point vibration test is performed, but a multipoint vibration test may be performed.

  Each of the accelerometers 23a to 23c is arranged at a measurement point on the analysis object 1, measures the magnitude of the acceleration generated in the analysis object 1 at the measurement point, converts this to an electrical signal, and outputs it. The electric signal is amplified by the amplifiers 25a to 25c and then output to the external interface 26. The measurement points are determined by the user so as to have an appropriate distribution and number to the extent that the eigenmodes up to the desired order can be sufficiently geometrically expressed. In the present embodiment, the triaxial accelerometers 23a to 23c are used as measuring instruments, but various parameters at measurement points such as a non-contact displacement meter, a laser Doppler velocimeter, a strain gauge, etc. using eddy currents or laser beams. Any measuring instrument may be used instead of the accelerometers 23a to 23c as long as the displacement of (position, velocity, acceleration, etc.) can be measured.

  The external interface unit 26 converts electrical signals input from the force transducer 21 and the accelerometer 23 via the amplifiers 24 and 25 into digital signals, and outputs the digital signals to the bus 39. The removable media device 31 is controlled by a removable media control unit 32. In the removable media device 31, a program (for example, a program for performing an analysis described later) recorded on the removable media 2 can be installed in the analysis unit 30, and an analysis result can be recorded on the removable media 2. As the removable medium 2, any storage medium such as FD, CD, DVD, and MO may be used.

  The input device 33 is an input device such as a keyboard, a mouse, and a digitizer, and is controlled by the input control unit 34. With this input device 33, the user can input various data such as the shape of the analysis object and the coordinate data of each measurement point. The output device 35 is an output device such as a display or a printer, and is controlled by the output control unit 36. The output device 35 outputs various data input to the input device 33, a menu for executing a program, an analysis result, and the like.

  The central processing unit 37 is connected to the external interface unit 26, the removable media control unit 32, the input control unit 34, the output control unit 36, and the storage device 38 via the bus 39, and gives instructions to these control units. Various processes such as Fourier transform and experimental characteristic matrix identification method are performed in accordance with each program stored in the storage device 129. The storage device 38 includes a RAM (Random Access Memory), a ROM (Read Only Memory), and the like, and stores programs, temporary data during execution of each program, analysis results, and the like.

  Next, an analysis process performed by the analysis apparatus 10 configured as described above will be described.

First, the basic concept about the analysis process performed by the analysis apparatus 10 of this embodiment is demonstrated. For example, the analysis object as shown in FIG. 2A is modeled by a multi-degree-of-freedom system. In particular, in the example shown in FIG. 2A, there are n measurement points, and the degree of freedom at each measurement point is three degrees of freedom in the x, y, and z directions. Modeled with a system of degrees of freedom. In this case, if the mass matrix, the viscous damping matrix, and the stiffness matrix are set to [M], [C], and [K], respectively, the modeled equation of motion of the analysis object is expressed by the following equation (1). In the following formula (1), viscous damping is assumed. However, structural damping may be assumed as necessary, or an equation of motion incorporating both dampings may be used.

Here, {ψ} in equation (1) is a vector in which x-direction displacement x _i , y-direction displacement y _i , and z-direction displacement z _i at each point determined as a measurement point are vertically arranged, and {f} is a measurement exciting force fx _i in the x direction at each point which defines a point represents a vector that can be side by side y-direction exciting force fy _i, the z-direction exciting force fz _i vertically (the following formula (2) see ).

  The mass matrix [M], the viscous damping matrix [C], and the stiffness matrix [K] (hereinafter collectively referred to as “characteristic matrix”) are 3n × 3n symmetric matrices (that is, the number of rows and the number of columns are Symmetric matrix equal to 3n degree of freedom of measurement of the analysis object), i indicates the number of the measurement point. In the equation (1), {ψ ′} is a first derivative of a vector representing displacement, that is, a vector representing the velocity of a measurement point, and {ψ ″} is a second derivative of a vector representing displacement, that is, a measurement point. Is a vector representing the acceleration of. In this specification, a matrix is represented by [] and a vector is represented by {}.

  When the characteristic matrix such as the mass matrix [M] is obtained in the equation of motion (1), the state of vibration of the analysis object 1 can be grasped, and the rigid body characteristics of the analysis object 1 (for example, the analysis object) 1 mass, center of mass position, three main moments of inertia, and directions of three main axes corresponding to these, etc.). Therefore, in the present embodiment, the characteristic matrix of the analysis object 1 is obtained based on the measurement data at each measurement point when the analysis object 1 is vibrated by the striking hammer 22.

  The identification of the characteristic matrix of the analysis object 1 and the calculation of the rigid body characteristics based on the measurement data at the time of vibration are performed roughly according to four processing steps as shown in FIG.

  First, as shown in step S10, mode characteristics are identified based on measurement data during a single-point excitation test. That is, one point of the analysis object 1 is struck (vibrated) by the striking hammer 22, this vibration force is measured by the force transducer 21, and the vibration generated by this vibration is measured by each accelerometer 23. The Thereafter, a frequency response function and a coherence function are calculated for each measurement point from the measurement data measured by the force transducer 21 and the accelerometer 23 in this way. Next, the natural frequency, natural mode vector, and mode damping ratio (hereinafter collectively referred to as “mode characteristics”) are identified by a general method using the calculated frequency response function.

  Next, in step S20, a physical constraint condition expression between characteristic matrix components is calculated. That is, the value of a specific component of the characteristic matrix and the constraint condition between specific components are based on various mechanical principles, the coordinates of each measurement point, the structural coupling state between the measurement points, etc. Are irrelevant. For example, the arrangement of zero components in the characteristic matrix can be specified to some extent based on the structural coupling state between measurement points. In addition, the relationship (constraint condition) between non-zero components in the characteristic matrix can be specified to some extent based on various mechanical principles, the coordinates of measurement points, and the like.

  Thus, by calculating the constraint equation between the characteristic matrix components, the number of components in the characteristic matrix that must be identified based on the mode characteristics obtained in step S10 is reduced. For example, when the mass matrix [M] is a matrix of 3n rows and 3n columns, the mass matrix [M] includes 3n (3n + 1) / 2 components (the actual number of components is 3n × 3n). It can be said that the number is 3n (3n + 1) / 2 in consideration of the symmetry of the matrix). Here, when the constraint condition between the components of the mass matrix [M] is not obtained, the components in the mass matrix [M] that must be identified based on the mode characteristics obtained in step S10. Is 3n (3n + 1) / 2, and in reality, it is difficult to appropriately identify all the components of the mass matrix [M].

  However, in this step S20, a constraint condition expression between the components of the mass matrix [M] is obtained, so that the arrangement of the zero component in the mass matrix [M] and the specific component different from the specific component in the mass matrix [M] are determined. The relational expression with the component of is obtained. Therefore, by calculating the constraint equation between the mass matrix [M] components, the number of components in the mass matrix [M] that must be identified based on the mode characteristics obtained in step S10 is considerably reduced. . Specifically, the number obtained by subtracting the number of zero components in the mass matrix [M] and the number of created constraint conditions from the number 3n (3n + 1) / 2 of all the components in the mass matrix [M]. Are identified on the basis of the mode characteristics obtained in step S10.

Next, in step S30, identification of independent unknown components of the mass matrix [M] and the stiffness matrix [K] is performed by a mathematical optimization method. That is, the relationship between the mass matrix [M] and the stiffness matrix [K], the natural angular frequency Ω _t and the natural mode vector {ψ} _t can be expressed as the following equation (3) (where t represents the order). .

As can be seen from Equation (3), when each component of the mass matrix [M] and the stiffness matrix [K] is regarded as a variable and its value is changed, the natural angular frequency Ω _t and the natural mode vector { The value of ψ} _t will also change. Accordingly, in this step S30, the natural angular frequency Ω _t and the natural mode vector {ψ} _t obtained by the above equation (3) are made to coincide with the natural angular frequency and the natural mode vector obtained in step S10. The values of the independent unknown components of the mass matrix [M] and the stiffness matrix [K] gradually increase while satisfying the constraint equation between the components of the mass matrix [M] and the stiffness matrix [K] obtained in step S20. It can be changed. Analysis target in which the mass matrix [M] and the stiffness matrix [K] when the natural angular frequencies and the natural mode vectors substantially match each other are modeled by a multi-degree-of-freedom system as shown in FIG. Identified as an object mass matrix [M] and a stiffness matrix [K].

  In the following, matching the natural angular frequencies and the natural mode vectors in this way is described as “matching” of the natural angular frequencies and the natural mode vectors, and the natural angular frequencies and the natural mode vectors are matched. The calculation and identification of the optimum values of each component of the mass matrix [M] and the stiffness matrix [K] will be described as “optimization” of the mass matrix [M] and the stiffness matrix [K].

  As a mathematical method for optimizing the mass matrix [M] and the stiffness matrix [K] in this way, many methods such as “the steepest descent method belonging to Newton's method” can be applied. Any method can be used as long as optimization can be performed efficiently and stably.

  Next, in step S40, based on the mass matrix [M] identified in step S30, the rigid body characteristics of the object 1 to be analyzed, that is, the mass of the object 1 to be analyzed, the center of mass position, the three main moments of inertia, and corresponding to these. The orientation of the three main axes is required.

  By the way, in the above-described characteristic matrix and rigid body characteristic identification method, an appropriate number n of measurement points are set on the analysis object 1, and vibrations in the x, y, and z directions at these measurement points during an excitation test are set. The rigid body characteristics are estimated by measuring the frequency response function and identifying the characteristic matrix of the corresponding degrees of freedom. When the analysis object 1 is a complex structure, it is necessary to set a certain number of measurement points depending on the complexity and size of the shape.

  However, in the above-described characteristic matrix and rigid body characteristic identification method, measurement point setting, accelerometer calibration, and measurement system are performed between the measurement point setting on the analysis object 1 and the calculation of the frequency response function. Therefore, the process from setting the measurement point to calculating the frequency response function is the most time-consuming and time-consuming work. Therefore, when the number of measurement points increases, the time and effort required to calculate the frequency response function increase accordingly, and as a result, the time required to identify the characteristic matrix and the rigid body characteristic increases. Further, since the frequency response function is obtained by the identification process performed in the central processing unit 37, when the number of measurement points increases, the work load in the central processing unit 37 increases, and the convergence speed of the frequency response function greatly decreases. . Therefore, from the viewpoint of quickly identifying the characteristic matrix and the rigid body characteristics, it is preferable that the number of measurement points to be actually measured is small.

  Therefore, in the embodiment of the present invention, the number of measurement points to be actually measured is reduced without reducing the number of measurement points, that is, without reducing the degree of freedom of the characteristic matrix. That is, a sufficient number of measurement points are set on the analysis object 1 as in the conventional case, and some of the measurement points are actually measured to calculate a frequency response function (hereinafter referred to as “this”). Such measurement points are referred to as “real measurement points”, and the remaining measurement points are not actually measured (hereinafter, such measurement points are referred to as “virtual measurement points”). In other words, in the present embodiment, a characteristic matrix having a degree of freedom larger than the degree of freedom of the measurement point at which measurement is actually performed is identified.

  In this embodiment, after setting actual measurement points and virtual measurement points in this way, in step S30, identification of independent unknown components of the mass matrix [M] and the stiffness matrix [K] is performed by a mathematical optimization method. In doing so, only the components of the eigenmode vector corresponding to the actual measurement point in the eigenmode vector are matched.

  Here, even if the inventor of the present application introduces the concept of virtual measurement points in this way, when a constraint equation is appropriately obtained based on the structural coupling state between the measurement points, etc. If only the component corresponding to the actual measurement point in the eigenmode vector can be matched, the component corresponding to the virtual measurement point can also be subordinately (or automatically) matched. It has been discovered that mass matrix [M] and stiffness matrix [K] can be identified. In other words, the inventor of the present application has found that the characteristic matrix can be optimized with sufficient accuracy without necessarily matching all eigenmode vector components when matching the eigenmode vectors. It was.

  Below, the procedure of the analysis process performed by the analyzer 10 which concerns on embodiment of this invention is demonstrated in detail. Note that the analysis processing in the present embodiment is also basically performed according to the processing procedure shown in FIG. Therefore, hereinafter, the processing procedure of the analysis processing according to the present embodiment will be described with reference to FIGS. 4 to 7 in a form corresponding to the processing procedure of the analysis processing.

[Identification of mode characteristics]
First, as shown in step S10, the mode characteristics are identified based on the measurement result of the single-point excitation test. FIG. 4 is a flowchart showing a processing procedure of mode characteristic identification processing according to the present invention. The identification of the mode characteristic is roughly classified into a vibration measurement process that performs a vibration test to calculate a frequency response function and the like, and a characteristic identification process that identifies the mode characteristic based on the calculated frequency response function and the like.

  Prior to the execution of the processing procedure shown in FIG. 4, the user sets the measurement system. In setting the measurement system, the analysis object 1 is arranged to be supported by an elastic body such as a rubber string or a tire tube by the user. Thereby, the surrounding free state of the analysis object 1 is approximately realized. Note that the analysis object 1 may be supported in any manner as long as the surrounding free state is realized in a pseudo manner.

  In the mode characteristic identification process as shown in FIG. 4, first, in step S <b> 11, the central processing unit 37 displays the output screen 35 so as to input initial data such as the shape and dimensions of the analysis object 1 to the user. To request through. Such initial data includes coordinate data of all measurement points including actual measurement points and virtual measurement points.

  The user sets the actual measurement point and the virtual measurement point at arbitrary positions on the analysis object 1. At this time, the virtual measurement points are set so that the virtual measurement points are not adjacent to each other as much as possible, that is, the measurement point closest to the virtual measurement point is the actual measurement point. In other words, the virtual measurement points are set to be arranged between the actual measurement points in at least one certain direction. The virtual measurement points are set so that the number of virtual measurement points does not become 1/2 or more of the total number of measurement points (preferably, it does not become 1/3 or more).

  FIG. 2B shows an example of setting actual measurement points and virtual measurement points. In FIG. 2B, circled measurement points indicate actual measurement points at which measurement is actually performed, and measurement points not circled indicate virtual measurement points at which measurement is not actually performed. ing. That is, in the illustrated example, the measurement is not actually performed at the measurement points Nos. 3, 6, 7, and 11.

  The user determines the setting position of each measurement point under such conditions, and inputs the coordinate data of all the measurement points in response to the request from the central processing unit 37. In addition, an accelerometer 23 is attached to the actual measurement point setting position.

Next, the user strikes (vibrates) the analysis object 1 at any one position (single point) of the positions where the accelerometer 23 is attached by the striking hammer 22.
The hit applied to the analysis object 1 by the hitting hammer 22 by the user is measured by the force transducer 21, and the measurement data is input to the external interface unit 26 via the amplifier 24. Further, the vibration generated in the object 1 to be analyzed by this impact is also measured by the accelerometer 23, and the measurement data is input to the external interface unit 26 via the amplifier 25. These input measurement data are taken into the central processing unit 37. In this way, a single point vibration test is performed.

  Next, in step S13, a frequency response function or the like is calculated based on the measurement data. That is, the central processing unit 37 performs a Fourier transform (for example, a fast Fourier transform) on the measurement data measured by the force transducer 21 and the accelerometer 23 to restore the frequency spectrum. Thereafter, the central processing unit 37 has a frequency spectrum related to the input to the analysis object 1 (frequency spectrum based on the output of the force transducer 21) and a frequency spectrum related to the output from the analysis object 1 (frequency based on the output of the accelerometer 23). The power spectrum and the cross spectrum are calculated from the spectrum), and the frequency response function and the coherence function are calculated from the power spectrum and the cross spectrum. These frequency response function and coherence function are calculated for each actual measurement point.

  Here, the coherence function (relevance function) is an inner product of an input, an output, and indicates a correlation between input and output. When the value is 0, the input and the output are uncorrelated, and the correlation increases as the value approaches 1. This coherence function is also used for identifying mode characteristics based on a frequency response function described later (used in step S15).

  From the viewpoint of obtaining more reliable data, it is preferable to repeat this single-point excitation test a plurality of times, average the individual power spectra and cross spectra, and use the average value.

  The processes in steps S11 to S13 are vibration measurement processes, and a frequency response function and a coherent function at each actual measurement point are calculated based on measurement data obtained by a single-point excitation test.

  As described above, after the frequency response function is calculated, the characteristic identification process for identifying the mode characteristic based on the calculated frequency response function or the like is performed.

  In the characteristic identification process, first, as shown in step S14, the output screen 35 is set so that the central processing unit 37 inputs a frequency band (hereinafter referred to as “identification frequency band”) for identifying the mode characteristic to the user. To request through. The user inputs the identification frequency band from the input device 33. As the identification frequency band is wider, the number of resonance peaks appearing in the band increases. For this reason, when the identification frequency band is widened, many order mode characteristics can be identified.

  Next, in step S15, mode characteristics up to a predetermined order are identified based on the frequency response function and coherent function calculated in step S13. This predetermined order changes in accordance with the identification frequency band input in step S14. Note that various known mode identification methods can be used as a method for identifying mode characteristics. Examples of such a mode identification method include a multipoint reference method (poly-reference method). , Partial iteration method, Proly's method, circle fit method and the like.

ここで、式（４）中のＹは慣性項定数パラメータ、Ｚは剰余項定数パラメータであり、また同定周波数帯域内に存在する共振峰の数をａとする。さらに、そのａ個の共振峰に対応する固有角振動数をΩ_t、モード減衰比をζ_t、加振点ｑ、応答点（測定点）ｐの固有モード成分をそれぞれψ_qt、ψ_pt（ｔは次数）とする。 Below, the basic principle at the time of identifying a mode characteristic based on a frequency response function is demonstrated easily. A frequency response function relating to a limited frequency band (that is, the identification frequency band) of the vibration system in which the damping is proportional viscosity damping can be approximately expressed by the following equation (4).

Here, Y in equation (4) is an inertia term constant parameter, Z is a remainder term constant parameter, and the number of resonance peaks existing in the identification frequency band is a. Furthermore, the natural angular frequency corresponding to the a resonance peaks is Ω _t , the mode damping ratio is ζ _t , the excitation point q, and the natural mode components of the response point (measurement point) p are ψ _qt and ψ _pt ( t is the order).

As described above, the frequency response function H _pq expressed by the equation (4) can be obtained from the vibration test in step S13 if the natural angular frequency Ω _t , the mode damping ratio ζ _t , and the natural mode components ψ _pt and ψ _qt are appropriately selected. The frequency response function based on Therefore, in this mode characteristic identification, the natural angular frequency Ω _t , the mode damping ratio ζ _t , such that the frequency response function H _pq calculated by the equation (4) substantially matches the frequency response function calculated in step S13, The eigenmode components ψ _pt and ψ _qt are calculated by an optimization method or the like, and these mode characteristics are identified. The coherence function is used in the process of identifying the mode characteristics.

  The process of steps S14 and S15 is a characteristic identification process, and the mode characteristic is identified based on the frequency response function.

[Identification of physical constraints between characteristic matrix components]
Next, as shown in step S20 in FIG. 3, a physical constraint condition expression between characteristic matrix components is calculated. FIG. 5 is a flowchart showing a processing procedure for calculating a physical constraint condition expression between characteristic matrix components of the present invention.

  In calculating the physical constraint condition expression between the characteristic matrix components, the central processing unit 37 requests via the output screen 35 to input the structural coupling state (connectivity) between the measurement points in step S21. Here, the structural coupling state between measurement points means whether a certain measurement point and another measurement point are structurally connected, that is, a certain measurement point and another measurement point are structurally adjacent to each other. It means the state of whether or not there is a matching relationship.

  For example, considering the structure shown in FIG. 2A, the structural coupling state between the measurement points is expressed as shown in FIG. In the figure, measurement points that are structurally connected to each other are connected to each other by straight lines. In the example shown in FIG. 2A, for example, the measurement point 1 is structurally coupled to the measurement points 2, 3, 4 and the measurement point 3 is structurally coupled to the measurement points 1, 2, 4 and 5. Is bound to. In the present embodiment, the structural coupling state between the measurement points is determined based on the subjectivity of the user based on the arrangement of the measurement points, the distance between the measurement points, and the like, and is input by the user using the input device 33. However, for example, the user inputs the shape of the analysis object 1 and the arrangement of the measurement points, and the central processing unit 37 calculates the structural coupling state between the measurement points based on the input data. Also good.

  Next, in step S22, the arrangement of zero components of the mass matrix [M] and the stiffness matrix [K] is specified based on the structural coupling state between the measurement points input in step 21. Hereinafter, the arrangement of the zero component will be described.

Here, a model in which a simple uniform rigid body element and a rotary spring as shown in FIG. Assuming that each rigid element has a uniform shape and has a length l _i , a center of mass, a mass m _i , and a rotational spring constant k _i , the degree of freedom of bending vibration expression (the number of measurement points in the above embodiment) Is equivalent to the translational displacement of points 1-5. In this case, the mass matrix [M] and the stiffness matrix [K] are represented by the following equations (5) and (6), respectively.

  Here, * in the equations (5) and (6) means that it is not a zero component, that is, a non-zero component. Thus, in the model as shown in FIG. 6, the mass matrix [M] is represented as a tridiagonal matrix, and the stiffness matrix [K] is represented as a five diagonal matrix. The reason for this will be described with reference to FIGS. In FIG. 6, the rotational elastic pair is drawn by increasing the dimension between the rigid elements, but it is assumed that there is no dimension between the rigid elements. That is, the right end of the rigid element I and the left end of the rigid element II are considered to be at the same position.

Referring to FIG. 6 (b), when a motion equation is created by a free body diagram for the center of mass of rigid body element I (centroid of rigid body element I), the rotational motion is translated as in equation (7). Is expressed as in equation (8). When these formulas (7) and (8) are collectively displayed as a determinant, the formula (9) is obtained.

ここで、剛体要素Ｉが均一形状であって均一材料密度分布であるとすると、下記式（１２）のように表せる。なお、式（１２）においてｍ₁は剛体要素１の質量である。

The conversion rule between the translational displacement degree of freedom at both ends of the rigid element I and the degree of freedom of displacement at the center of mass can be expressed as the following equation (10). Therefore, substituting Equation (10) into Equation (9) and multiplying by [T ₁ ] ^T from the left side of both sides yields Equation (11) below. [T ₁ ] ^T represents a transposed matrix of T ₁ .

Here, assuming that the rigid element I has a uniform shape and a uniform material density distribution, it can be expressed as the following formula (12). In Equation (12), m ₁ is the mass of the rigid element 1.

Similarly, with reference to FIG. 6C, the equation of motion is created and arranged in the same way as in the case of the rigid element I by the free body diagram with respect to the center of mass of the rigid element II (element centroid). (13) is derived.

In this way, the equations of motion of rigid body element III, rigid body element IV,... Are obtained in order, and all the motion equations are combined into a matrix expression to create the motion equation of the entire system, and the following equation (14) is obtained. Thus, the mass matrix [M] is a tridiagonal matrix, and the stiffness matrix [K] is a pentadiagonal matrix. Thus, the arrangement of the non-zero components of the mass matrix [M] and the stiffness matrix {K] can be determined by the structural coupling state between the measurement points.

Considering the model of the object to be analyzed that is in a structurally coupled state between the measurement points as shown in FIG. 2C based on the dynamic principle described above, for example, the arrangement of zero components of the stiffness matrix [K] Becomes the following equation (15).

  [N] in the stiffness matrix [k] of the above formula (15) is a small matrix portion that may be a non-zero component of 3 rows and 3 columns, and [0] is a small matrix portion that is a zero component of 3 rows and 3 columns. Represents the matrix part. Here, [N] and [0] are small matrices of 3 rows and 3 columns because the measurement degrees of freedom at each measurement point are 3 degrees of freedom in the x, y, and z directions. If the measurement freedom at the measurement point is 6 degrees of freedom in addition to the x, y and z directions plus the rotation directions around the x, y and z axes, [N] and [0] are 6 rows 6 It becomes a small matrix of columns.

  In the example shown in FIG. 2C, the measurement point 1 is physically coupled to the measurement points 2 and 3, and these measurement points 2 or 3 are physically coupled to the measurement points 4 and 5. Accordingly, as shown in the equation (17), in the stiffness matrix [K], the non-zero small matrix components are arranged in the first column corresponding to the measurement point 1 in the first to fifth rows. In this way, the zero component arrangement of the characteristic matrix is automatically determined by the computer based on the dynamic principle according to the structural coupling model designed by the user.

Next, in step S23, for the mass matrix [M] of the three-dimensional structure, the mass matrix [M] for any n-degree-of-freedom model can be arbitrarily set in the structure by the calculation of the left side of the following equation (16). Based on the dynamic principle that it can be converted to a 6-by-6 rigid mass matrix M _rigid ] regarding 6-degree-of-freedom displacement of one point as a set point representing the rigid body motion, a set of non-zero components of the mass matrix [M] A linear equation is created as a physical constraint equation. Here, the matrix [Φ] in the equation (16) is a 3n × 6 matrix (hereinafter referred to as “rigid body displacement mode matrix”) represented by the following equation (17), and the rigid mass matrix [M _rigid ]] is a 6 × 6 matrix represented by the following formula (18).

In the rigid mass matrix, I _xx , I _yy , and I _zz are moments of inertia around axes that pass through the rigid body motion expression points and are parallel to the coordinate axes, and I _yx , I _zx , and I _yz are inertias between these three axes. It is a product. A, B, and C are values resulting from a deviation between the rigid body motion expression point and the center of gravity. The mass of the analysis object 1 is m, the center of gravity coordinates are (x _G , y _G , z _G ), and the rigid body motion expression point. When coordinates are (x _R , y _R , z _R ), A = m (x _G −x _R ), B = m (y _G −y _R ), and C = m (z _G −z _R ). [M _rigid ] is a symmetric matrix as shown in the above equation (18).

Hereinafter, the derivation method of the rigid body displacement mode matrix [Φ], the rigid body mass matrix [M _rigid ] and the above equation (18) will be briefly described.
The rigid body displacement mode matrix [Φ] is composed of six independent rigid body motion vectors. Specifically, using the measurement point coordinate value that is a part of the input data, the measurement point of the measurement point on the structure when the rigid body motion expression point undergoes unit translational displacement in the x-axis direction, the y-axis direction, and the z-axis direction. The vectors representing the displacement are applied to the first, second, and third columns, respectively, and the rigid body motion expression points pass through the points around the axis parallel to the x axis, around the axis parallel to the y axis, and parallel to the z axis. By applying a vector that linearly expresses the translational displacement of the measurement point on the structure relating to the rigid body angular displacement motion that is slightly angularly displaced about an arbitrary axis to the fourth, fifth, and sixth columns, the rigid body displacement mode matrix [Φ ] Can be easily created.

For example, referring to FIG. 7, the rigid body motion expression point coordinates set in the rigid body are R (x _R , y _R , z _R ), and the x, y, z coordinates of the measurement point i are (x _i , y _i , z _i ), when the rigid body is simply displaced in translation in the x-axis direction by δ _Rx , the displacement at the measurement point i is accompanied by a transformation matrix (this is a constant composed of δ _Rx as shown in the following equation (19)). In this case, it can be expressed by multiplication by a column vector).

The displacement of the measurement point i by rotating the rigid body to the right hand by the angle θ _x around the axis parallel to the x axis at the rigid motion expression point is expressed algebraically from the geometry as shown in the following formula (21). Can do.

Here, considering that the translational amplitude and the angular displacement vibration of the target mechanical vibration are small values, approximations of sin θx≈θx and cos θx≈1 can be sufficiently applied. Can be transformed into the following equation (21). Similar to Equation (19), the displacement of the measurement point i can be expressed by multiplying the angular displacement of the motion expression point by a transformation matrix composed of constants based on the coordinate data.

Considering all rigid body motions with 6 degrees of freedom in this way, the component structure of the rigid body displacement mode matrix [Φ] is expressed in the following formula (22) as a form in which only the x, y, z direction displacement components of the measurement point i are specified. ) Can be created.

  Therefore, the rigid body displacement mode matrix [Φ] can be obtained by inputting the coordinate data of each measurement point. The experimental identification methods such as Fourier transform and mode characteristic identification method are all data processing methods (signal processing methods), and essentially require measurement point coordinates as geometric data (or structural data). It can be considered that there has never been an "experimental identification method", and this method has a great feature in this respect.

Next, the rigid mass matrix [M _rigid ] will be described. It is known that the rigid mass matrix [M _rigid ] on the right side of the above equation (16) has a unique component arrangement relationship as in the above equation (18) for any target structure. For example, in the rigid body mass matrix [M _rigid ], the 1 row 1 column component, the 2 row 2 column component, and the 3 row 3 column component are always the same value, and the value is the value of the mass of the analysis object. The 2 row 1 column component and the 1 row 2 column component are always zero, and the 5 row 1 column component and the 4 row 2 column component have the same absolute value and opposite signs. There is a significant component arrangement relationship.

The rigid body displacement mode matrix [Φ], the rigid body mass matrix [M _rigid ], and the mass matrix [M] expressed as described above have a relationship as expressed by the above equation (16). Here, the values of the non-zero components (m, A, Ixx, etc.) themselves of the equation (18) are unknown before identification, but using the relationship between these non-zero components, It is possible to automatically generate eleven equality constraint expressions related to the components of the left-side mass matrix [M].

  If a rigid body characteristic value such as the mass (weight) of the analysis object 1 is given as a known value, this is substituted into equation (18) to generate a constraint equation between mass matrix components. Also good. If a known rigid body characteristic value is obtained in this way, there is an advantage that the known rigid body characteristic value need not be expressed by other components.

  In the simultaneous equations generated in this way, the number of non-zero components that are unknowns is always greater than the equations independent of each other. Therefore, it is impossible to specify all the non-zero components of the mass matrix [M] only from the simultaneous equations generated in this way. However, as many unknowns (non-zero components) as the number of equations can be dependent variables of other unknowns (non-zero components). In the following description, a component expressed in a form dependent on other components of the characteristic matrix by the constraint expression obtained in step S23 and the like is a dependent component, and a form dependent on another component also by the constraint expression The component which is not expressed is demonstrated as an independent uncomponent.

ここで、ｍ_i,jは質量行列のｉ行ｊ列の成分を表し、左辺の成分が従属成分、右辺の各項の成分が独立未知成分である。なお、どの成分が独立未知成分に割り当てられ、どの成分が従属成分に割り当てられるかはプログラミングに依存する。例えば、式（２３）の最上段の式の場合、ｍ_1,1を独立未知成分とすることもでき、この場合ｍ_2,1が従属成分となる。 When such an equation constraint conditional expression is automatically generated by the above method, an equation such as the following expression (23) is created.

Here, m _{i, j} represents the component of i-th row and j-th column of the mass matrix, the component on the left side is the dependent component, and the component of each term on the right side is the independent unknown component. Note that which component is assigned to an independent unknown component and which component is assigned to a dependent component depends on programming. For example, in the case of the uppermost expression of Expression (23), m _1,1 can be an independent unknown component, and in this case, m _2,1 is a dependent component.

Next, in step S24, with respect to the stiffness matrix [K] and the damping matrix [C] of the three-dimensional structure, no strain occurs in any position of the structure in the rigid body motion state, and all the stresses based on the strain are applied. Based on the basic dynamic principle of being zero in the part, simultaneous linear equations are created as physical constraint equations between the non-zero components of the stiffness matrix [K] and the damping matrix [C]. From the dynamic principle, the following formula (24) is derived. Note that [0] on the right side of the equation (24) is a zero matrix of n rows and 6 columns.

  By this equation (24), simultaneous linear equations relating to non-zero components in each characteristic matrix are created. As a result, as in step S23, a plurality of equality constraint expressions relating to the components are generated for the stiffness matrix [K] and the viscosity matrix [C], and the number of independent unknown components is reduced by the number of generated equations. be able to.

  As described above, in steps S21 to S24, a constraint condition expression indicating the relationship between each component is obtained for each characteristic matrix based on the physical principle or the like. As a result, the number of independent unknown components included in each characteristic matrix can be reduced, and dynamic (physical) component related characteristics can be added, and the optimization accuracy in the optimization process in step S30 can be increased. Can do.

[Identification of characteristic matrix]
Next, as shown in step S30 in FIG. 3, the identification of independent unknown components in the mass matrix [M], the stiffness matrix [K], and the viscosity matrix [C], and thus the identification of these characteristic matrices are performed. As described above, many methods can be adopted as a mathematical method for optimizing the mass matrix [M] and the stiffness matrix [K], and one of them is described below with reference to FIGS. One example will be described. 8 and 9 are flowcharts showing the processing procedure of the characteristic matrix identification processing of the present invention.

The characteristic matrix identification process of this embodiment is roughly divided into three processes. The first processing is such that the mass matrix [M] into which an appropriate value is assigned as an independent unknown component becomes a positive definite matrix, and the stiffness matrix [K] into which an appropriate value is assigned as an independent unknown component. Is a positive definite value / non-negative definite value processing for correcting the value of the independent unknown component so that becomes a non-negative definite matrix. In the second process, the natural angular frequency Ω _t and the natural mode vector {ψ} _t obtained by the above equation (3) are converted into the natural angular frequency and the natural mode vector (hereinafter referred to as the natural frequency vector obtained in step S10). each "target natural angular frequency Omega _tt", to match the "called target eigenmode vector {ψ} _tt), gradually changes the value of the independent unknown component of the mass matrix [M] and stiffness matrix [K] The third process is to multiply the normalized mass matrix [M] and stiffness matrix [K] by the gain translation coefficient α, and finally the mass matrix [M]. 8 and 9, each of these three large processes is surrounded by a broken line.

  In the identification process of the characteristic matrix, first, in step S31, the physical constraint equation is created in step S20 (specifically, step S24) among the components of the mass matrix [M] and the stiffness matrix [K]. A random number is assigned as an initial value to the component that becomes the component. Further, the dependent component is calculated using the constraint condition equation created in step S20 based on the assigned random number, and the mass matrix [M] and the stiffness matrix [K] as initial values are created as numerical matrices.

  However, since the mass matrix [M] and the stiffness matrix [K] as initial values are created from random numbers, there are cases where the conditions that must be physically satisfied are not satisfied. That is, since the mass and the rotational inertia are always positive values, the mass matrix is a positive definite matrix. Similarly, the stiffness matrix is basically a positive value, but since the strain is zero in rigid body motion, there are up to six eigenvalues (zero eigenvalues exist for the rigid body displacement degree of freedom). Therefore, the stiffness matrix is a non-negative definite matrix. However, the mass matrix [M] and the stiffness matrix [K] as the initial values may not satisfy such a condition.

  Therefore, in this embodiment, in steps S32 to S34, the mass matrix is converted to a positive definite matrix, and the stiffness matrix is converted to a non-negative definite matrix. In other words, for the mass matrix [M], the independent unknown component of the mass matrix [M] is corrected so that all 3n eigenvalues become positive values, and all 3n eigenvalues of the stiffness matrix [K] are negative. The independent unknown component of the stiffness matrix [K] is corrected so as not to become a value. Here, as a representative, a method for forming a positive definite matrix of a mass matrix will be described.

このようにして算出された固有値λの中に負値のものがなければ、質量行列［Ｍ］はこの時点の各成分の値で正定値行列である。一方、固有値λに負値のものがある場合には、これら負値の固有値λを正値へ変化させる必要がある。このために、負値の固有値λを正値へ変化させるための１次感度解析を行う。この感度解析としては、例えば、最急降下法、共役勾配法といったニュートン法に代表されるニュートン法類の手法等、一般に知られている各種の感度解析の手法を利用することができる。以下では、その一例について示す。質量行列［Ｍ］内の独立未知成分ｐ_iに関する固有値λについての１次感度ｓ_t,iは上記式（２５）を独立未知成分ｐ_iで微分した下記式（２６）で算出される。

First, after normalizing the mass matrix [M] with respect to the maximum absolute value of the component, that is, with the maximum absolute value of the component of the mass matrix [M] being 1, The following equation (25), which is an eigenvalue problem, is analyzed to obtain all 3n eigenvalues λ and eigenvectors {x}, and negative eigenvalues are selected.

If there is no negative value among the eigenvalues λ calculated in this way, the mass matrix [M] is a positive definite matrix with the values of the respective components at this time. On the other hand, when there are negative eigenvalues λ, it is necessary to change these eigenvalues λ to positive values. For this purpose, a primary sensitivity analysis is performed to change the negative eigenvalue λ to a positive value. As this sensitivity analysis, for example, various well-known sensitivity analysis methods such as Newton's methods such as the steepest descent method and the conjugate gradient method can be used. Below, the example is shown. The primary sensitivity s _{t, i} for the eigenvalue λ relating to the independent unknown component p _i in the mass matrix [M] is calculated by the following equation (26) obtained by differentiating the above equation (25) with the independent unknown component p _i .

Here, λ _t and {x} _t in the equation (26) indicate a t-order eigenvalue and an eigenvector, respectively. Further, p _i is one of q independent unknown components in the mass matrix [M] (i = 1 to q). Increment Δλ of eigenvalue λ _t when q independent unknown components are increased by a small amount Δp _i (i = 1 to q) using the primary sensitivity s _{t, i} calculated by the above equation (26). _t can be estimated approximately by the following equation (27).

Therefore, this equation is generated for all negative eigenvalues λ. At this time, for example, a positive value larger than the absolute value of the negative eigenvalue is substituted into the eigenvalue increment Δλ _t . Then, to calculate the fine increment Delta] p _i independent unknown component that satisfies this way is generated wherein the updating of the independent unknown component by using a very small increment Delta] p _i of the thus independent unknown component calculated by Do. By repeating such an operation, the mass matrix [M] can be made positive definite. Note that the non-negative definite value of the stiffness matrix [K] is also performed in the same manner as the positive definite value of the mass matrix [M].

Based on this concept of positive definite value, in step S32, a negative value is selected from all the eigenvalues λ _t calculated based on the above equation (25). Next, in step S33, eigenvalue primary sensitivity st _{, i} for each independent unknown component is calculated for each negative eigenvalue λ selected in step S32 using the above equation (26). Then, in step S34, by using the first-order sensitivity calculated in step S33 is calculated minute increment Delta] p _i independent unknown component by the above formula (27), the following equation using the small increment Delta] p _i calculated (28 ) Updates the independent unknown component.

  Next, in step S35, it is determined whether or not the positive definite value of the mass matrix [M] and the non-negative definite value of the stiffness matrix [K] are realized as a result of updating the independent unknown component in step S34. If it is determined that the matrix [M] is not positively definite and the stiffness matrix [K] is not negatively definite, steps S32 to S34 are repeated. On the other hand, when it is determined in step S35 that the positive definite value of the mass matrix [M] and the nonnegative definite value of the stiffness matrix [K] are realized, the process proceeds to step S36, that is, the optimization process. The processing from step S32 to step S35 corresponds to the positive definite value / non-negative definite value processing described above.

In the characteristic matrix identification process, an optimization process is performed next. In the optimization process, the natural angular frequency (undamped natural angular frequency) Ω _t and the natural mode vector {ψ} _t calculated by solving the general eigenvalue problem of the following equation (29) were obtained through experiments. A process of matching the natural angular frequency (target natural angular frequency Ω _tt ) and the natural mode vector (target natural mode vector {ψ} _tt ) is performed. For this reason, in this embodiment, a primary sensitivity analysis is performed in order to match these natural angular frequencies and natural mode vectors. Also in this sensitivity analysis, various known sensitivity analysis methods can be used.

The first-order sensitivity for the natural angular frequency Ω is obtained by differentiating the above equation (29) with respect to the independent unknown component p _i and modifying the _t-th natural angular frequency Ω _t and the natural mode vector {ψ} _t. It is calculated based on the following formula (30) obtained.

Here, p _i in formula (30) is assumed to be one of the q Separate unknown component into the mass matrix [M] and stiffness matrix [K] (i = 1~q) . The eigenvector sensitivity is obtained by performing a differentiation operation using the above equation (29), the normalization equation for the mass matrix of the eigenvector (the following equation (31)), and the equation for the orthogonality condition of the eigenmode (the following equation (32)). Is obtained by the following equation (33).

Using the first-order sensitivity obtained in this way, the increment of the square value of the natural angular frequency when q independent unknown components are increased by a small amount Δp _i (i = 1 to q), Steps based on a value obtained by squaring the natural angular frequency Ω _t calculated by solving the general eigenvalue problem of the above equation (29) based on the current mass matrix [M] and stiffness matrix [K] and an excitation test When q independent unknown components are increased by a small amount Δp _i so as to be a difference (Ω _tt ² −Ω _t ² ) from the square value of the target natural angular frequency Ω _tt calculated in S15 Eigenmode vector {ψ} _t and excitation obtained by solving the general eigenvalue problem of Equation (29) based on the current mass matrix [M] and stiffness matrix [K] Based on the test, the target eigenmode vector calculated in step S15 { [psi} _tt and the difference - so that _{({ψ} tt {ψ}} t), increment Delta] p _i of the small amount is calculated (the following equation (34)).

Then, to update the independent unknown component by the same formula as the formula (28) using a very small increment Delta] p _i of the thus independent unknown component calculated by. By repeating such operations, the mass matrix [M] and the stiffness matrix [K] are optimized, and the components of the mass matrix [M] and the stiffness matrix [K] are identified.

In the lower expression of the above expression (34), in the present embodiment, only the eigenmode component corresponding to the actual measurement point degree of freedom is arranged in the column vector (target eigenmode vector) {ψ _tt } on the left side, and the virtual The eigenmode component corresponding to the measurement point degree of freedom is not arranged. Therefore, in the present embodiment, only the eigenmode component corresponding to the actual measurement point degree of freedom is matched with the target value by the above equation (34), and the eigenmode component corresponding to the virtual measurement point degree of freedom is set to the target value. Do not match to. That is, in this embodiment, the degree of freedom corresponding to the degree of freedom of the actual measurement point from among the eigenmode vectors composed of 3n components obtained by eigenvalue analysis of the mass matrix [M] and stiffness matrix [K] to be identified. Only the components are extracted to create a column vector {ψ _t } to optimize the above equation (34).

  As described above, the mass matrix [M] and the stiffness matrix [K] can be appropriately optimized even when the eigenmode components corresponding to the actual measurement point degrees of freedom are matched with the target values. That is, in the present embodiment, a constraint condition expression based on the structural coupling state between the measurement points obtained in the above steps S23, S24, etc. is used. For this reason, if the actual measurement points are sufficiently adequately arranged on the analysis object 1 so that the shape of the eigenmode can be appropriately expressed, the displacement of the virtual measurement points arranged between the actual measurement points and the eigenmode components Can be expressed by interpolation based on the displacement of the actual measurement point and the eigenmode component. Therefore, if only the eigenmode components corresponding to the actual measurement point degrees of freedom are matched to the target value, the eigenmode components corresponding to the virtual measurement point degrees of freedom arranged between the actual measurement points are subordinately appropriate values. Will be adjusted to be.

  Note that matching is performed for eigenmode components corresponding to the degrees of freedom of all measurement points set on the analysis object 1 by performing matching for only the eigenmode components corresponding to the actual measurement point degrees of freedom. Compared to the case, the calculation load of the optimization analysis of the mass matrix [M] and the stiffness matrix [K] is greatly reduced, and thus the calculation speed of the identification analysis can be improved. In particular, since a part of the measurement points is set as the virtual measurement point as described above, the time required for setting the accelerometer 23 and the like can be significantly reduced. Therefore, in this embodiment, the analysis time can be greatly shortened. .

  In addition, since a sufficient degree of freedom can be ensured by arranging the virtual measurement points, the shape of the eigenmode can be expressed smoothly, so the degrees of freedom of all the measurement points set on the analysis object 1 The identification accuracy of the mass matrix [M] and the stiffness matrix [K] hardly deteriorates when matching is performed for the eigenmode components corresponding to. Therefore, according to the present embodiment, the analysis time can be significantly reduced without substantially reducing the identification accuracy of the mass matrix [M] and the stiffness matrix [K].

Based on such an optimization concept, in step S36, the primary sensitivity of each natural frequency for each independent unknown component is calculated for each natural frequency Ω by the above equation (33). Next, in step S37, by using the first-order sensitivity calculated in step S36 is calculated minute increment Delta] p _i independent unknown component by the above formula (34), said equation using the small increment Delta] p _i calculated (28 The independent unknown component is updated by the same formula as ().

In step S38, the natural angular frequency Ω _t and the natural mode vector {ψ} _t calculated based on the mass matrix [M] and the stiffness matrix [K] including the independent unknown component updated in step S37 are used as the target natural angle. It is determined whether the frequency Ω _tt and the eigenmode vector {ψ} _tt match within an allowable range. The determination in step S38 is performed by calculating a value indicating the correlation between the mode characteristic obtained from the experimental result and the mode characteristic obtained by the calculation, and determining whether or not this value satisfies the allowable range.

  If it is determined in step S38 that the natural angular frequencies and the natural mode vectors do not match within the allowable range, steps S36 and S37 are repeated. On the other hand, when it is determined that the natural angular frequencies and the natural modes match within the allowable range, the process proceeds to step S39, that is, the gain adjustment process. The processing from step S36 to step S38 corresponds to the optimization processing described above.

  In the characteristic matrix identification process, a gain adjustment process is performed next. In the identification method of the mass matrix [M] and the stiffness matrix [K] by the above-described iterative calculation, the normalization condition for setting the maximum absolute value of the components of the mass matrix [M] to 1 is set for convenience of mathematical processing of optimization. Optimization processing is performed while holding. For this reason, the frequency response function corresponding to the measurement result of the vibration test calculated using the characteristic matrix obtained in the optimization process is the frequency response function obtained by the vibration test and the frequency of the resonance peak ( That is, the natural angular frequency) and the natural mode match, but the gain (vibration amplitude level) does not match. However, as shown in FIG. 10, the frequency response function (solid line X in FIG. 10) calculated using the characteristic matrix obtained by the optimization process and the frequency response function obtained by the vibration test (FIG. 10) The broken line Y) in FIG. 10 is in such a relationship that when the mass matrix [M] obtained by the optimization process is appropriately multiplied by a constant, they coincide with each other. Therefore, an optimization operation is performed so that the gains are also optimally matched.

  Here, in order to make the frequency response function X based on the characteristic matrix coincide with the frequency response function Y based on the vibration test, the values of both frequency response functions at a specific frequency are calculated, and based on the ratio of the calculated values. It is conceivable to calculate the gain. However, since the change in the value of the frequency response function is large in the frequency band near the resonance frequency (Z in FIG. 10), the gain is calculated based on the ratio of the values of both frequency response functions in the frequency band other than the vicinity of the resonance frequency. Is preferred.

Therefore, in this embodiment, the value of the frequency response function obtained by the following equation (35) in the frequency band excluding the vicinity of the resonance frequency optimally matches the value of the frequency response function actually obtained in the vibration test. In this way, the value of the gain adjustment parameter α is determined. In Expression (35), {f (ω)} is an external force vector for calculating the frequency response function at the frequency ω, and {h (ω)} is a vector representing the frequency response function at the frequency ω.

  Here, since not all measurement points are actually measured in the vibration test, the frequency response function based on the vibration test is calculated while the frequency response function based on the characteristic matrix is obtained for all measurement points. It is only required for actual measurement points. Therefore, in the present embodiment, the value of the gain adjustment parameter α is determined so that only the frequency response function based on the characteristic matrix corresponding to the actual measurement point matches the frequency response function based on the vibration test.

  Finally, the gain adjustment parameter α calculated in this way is multiplied by the mass matrix [M] and the stiffness matrix [K] calculated by the optimization process, and the matrix (α [M] and Let α [K]) be the mass matrix and stiffness matrix as the final identification results, respectively.

  Based on such a concept of gain adjustment, in step S39, a frequency response function is calculated using equation (35) based on the mass matrix [M] and stiffness matrix [K] identified in steps 36 to 38. The Next, in step S40, the gain response parameter α is calculated by comparing the frequency response function calculated in step 39 with the frequency response function obtained in the vibration test. Next, in step S41, the mass matrix [M] and the stiffness matrix [K] identified in steps 36 to 38 are multiplied by α to calculate a mass matrix and a stiffness matrix as a final identification result ([M ] = Α [M], [K] = α [K]).

  As described later, since only the mass matrix [M] is necessary for identifying the rigid body characteristics, the attenuation matrix [C] is not identified in the above embodiment. However, since the attenuation matrix [C] can be identified, the identification method of the attenuation matrix [C] will be briefly described below.

In identifying the damping matrix [C], first, the initial value of the damping matrix [C] is created by copying the stiffness matrix [K] and [C] = β [K]. Then, the optimum β is determined by the least square method with respect to the mode damping ratio ζ _t from the first order to the predetermined order to be optimally matched. That is, the correction coefficient is adjusted so that the mode damping ratio ζ _t calculated by the following equation (36) _matches the mode damping ratio ζ _tt identified in step S15 based on the vibration test for all orders to be matched. β is calculated.

Using the correction coefficient β calculated in this way, a temporary attenuation matrix is created by [C] = β [K]. Next, for further optimization, the mode attenuation ratio ζ _t calculated based on the attenuation matrix is determined based on the vibration test using the attenuation matrix thus created as an initial value. The primary sensitivity analysis is performed so as to match the ratio ζ _tt . Also in this sensitivity analysis, various known sensitivity analysis methods can be used.

The primary sensitivity for the mode damping ratio is calculated by the following equation (37) with respect to the independent unknown components p _i (i = 1 to q) constituting the damping matrix for the _t-th mode damping ratio ζ _t . Using the first-order sensitivity calculated in this manner, the increment of the mode attenuation ratio when q independent unknown components are increased by a small amount Δp _i (i = 1 to q) is the current attenuation matrix [C ] be the difference between the calculated modal damping ratio ζ _tt (ζ _tt -ζ _t) based on the mode pressure and damping ratio zeta _t vibration test is calculated based on the incremental Delta] p _i of the small amount of Is calculated (the following equation (38)).

Then, to update the independent unknown component by the same formula as the formula (28) using a very small increment Delta] p _i of the thus independent unknown component calculated by. Then, by repeating such operations, the attenuation matrix [C] is optimized, and the components of the attenuation matrix [C] are identified.

[Identification of rigid body properties]
Next, as shown in step 50 in FIG. 3, the rigid body characteristics are calculated based on the mass matrix [M] identified based on the identification process. FIG. 11 is a flowchart showing a processing procedure of rigid body characteristic calculation processing based on the mass matrix [M] of the present invention.

  By the way, when the mass matrix [M] identified based on the identification process is substituted into the left side of the equation (16) and calculated, a 6 × 6 matrix (rigid mass matrix) having a component configuration shown on the right side is calculated. Is obtained as a numerical matrix. Each component of the rigid mass matrix has the component configuration shown in the above equation (18). Therefore, the mass m can be directly obtained based on the rigid mass matrix calculated in this way. Further, when the rigid body displacement mode matrix [Φ] is created using the rigid body motion expression point as the coordinate origin and the left side of the above equation (16) is calculated, the coordinates of the center of mass (A / m, B / m, C / m) Is also calculated directly.

On the other hand, the calculation of the main moment of inertia and the direction of the main spindle of inertia is performed as follows. The rigid body rotational motion vectors in the 4th to 6th columns of the rigid body displacement mode matrix [Φ] are expressed not by the rigid body motion expression point but by the coordinates of the center of mass calculated as described above, and the coordinate origin of the above equation (16). Perform the left-hand side operation again. Then, the rigid mass matrix [M _rigid ] of the calculation result is as shown in the following formula (39).

For the so-called inertia tensor (3 rows and 3 columns) from the 4th row and the 4th column to the 6th row and the 6th column of the equation (39), the standard eigenvalue problem of the following equation (40) is solved. The three eigenvalues calculated in this way are the values of the three principal moments of inertia of the object 1 to be analyzed, and the eigenvectors normalized to the length 1 obtained corresponding to them are the inertial principal axes corresponding to the respective principal moments of inertia. It is a direction cosine vector (representing the direction of the principal axis of inertia). That is, the following equation (41) is obtained as a solution to the eigenvalue problem of equation (40).

Based on such a concept of mode characteristic identification, in step S51, a rigid body displacement mode matrix [Φ] created using the mass matrix [M] identified based on the identification process and the rigid body motion expression point as a coordinate origin. Based on the above, the mass m and the coordinates of the mass center (A / m, B / m, C / m) are calculated using the above equation (16). Next, in step S52, based on the mass matrix [M] identified based on the identification process and the rigid body displacement mode matrix [Φ] created using the coordinates of the center of mass calculated in step S51 as the origin, A rigid mass matrix [M _rigid ] is calculated using Equation (16). In step S53, based on the rigid mass matrix [M _rigid ] calculated in step S52, the values of the three main inertia moments of the object 1 to be analyzed, and the direction cosine vector of the inertia main axis corresponding to each main inertia moment are calculated. The

  According to the embodiment of the present invention, as described above, the identification of the characteristic matrix relating to the vibration of the analysis object 1 and the rigid body characteristics from the multipoint frequency response function based on the measurement data of the single-point excitation test of the analysis object. Can be identified. Therefore, in the embodiment of the present invention, since the characteristic matrix can be identified experimentally, a physical model of an analysis object to be directly subjected to vibration analysis can be expressed. Therefore, it can be directly applied to the partial structure synthesis method, optimum design, integration with the finite element method, various simulations, and the like.

  In particular, in the embodiment of the present invention, by setting some of the measurement points as virtual measurement points, a characteristic matrix having a degree of freedom larger than the measurement degree of freedom of the actual measurement points to be actually measured is identified. . As a result, as described above, the analysis time can be greatly reduced without substantially reducing the identification accuracy of the characteristic matrix and the rigid body characteristics. In other words, the identification accuracy of the characteristic matrix and the rigid body characteristic can be improved as compared with the case where the characteristic matrix having the same degree of freedom as the actual measurement point to be actually measured is identified.

  In the above-described embodiment, a characteristic matrix having a degree of freedom greater than the actual degree of freedom of measurement is identified by using some of the measurement points as virtual measurement points. However, the present invention can be applied if the degree of freedom of the identified characteristic matrix is larger than the degree of freedom of measurement actually measured in the vibration test.

  Therefore, for example, a characteristic matrix of 3n rows and 3n columns is identified based on n measurement points having two degrees of freedom (for example, only in the x-axis direction and the y-axis direction) at each measurement point, It is also possible to identify a characteristic matrix of 6n rows and 6n columns based on n measurement points whose measurement degrees of freedom of points are 3 degrees of freedom (for example, x-axis direction, y-axis direction, and z-axis direction). is there. The characteristic matrix of 6n rows and 6n columns is a characteristic matrix created on the assumption of rotational motions about the x, y, and z axes in addition to translational motions in the x, y, and z axis directions.

  In other words, according to the present invention, a characteristic matrix of degrees of freedom that combines the actual measurement degrees of freedom in which measurement is actually performed and the virtual measurement degrees of freedom in which measurement is not actually performed is used as the actual measurement freedom in which measurement is actually performed. The present invention can also be applied to the case where identification is performed based on measurement data of degrees.

It is a figure which shows the block diagram of the analyzer of this embodiment. It is a figure which shows the example of modeling in the multi-degree-of-freedom system of an analysis target object. It is a flowchart of an analysis process. It is a flowchart of the identification process of a mode characteristic. It is a flowchart which shows the process sequence of the calculation process of a physical constraint conditional expression. It is a figure for demonstrating the relationship between a structural coupling | bonding state and arrangement | positioning of a zero component. It is a figure for demonstrating a rigid body displacement mode matrix and a rigid body mass matrix. It is a part of flowchart which shows the process sequence of the identification process of a characteristic matrix. It is a part of flowchart which shows the process sequence of the identification process of a characteristic matrix. It is a figure for demonstrating a gain adjustment process. It is a flowchart which shows the process sequence of the calculation process of a rigid body characteristic.

Explanation of symbols

DESCRIPTION OF SYMBOLS 1 Analysis target object 10 Analyzing device 20 Measuring part 21 Force transducer 22 Impact hammer 23 Accelerometer 24, 25 Amplifier 26 External interface part 30 Analyzing part 31 Removable media device 33 Input device 35 Output device 37 Central processing unit 38 Storage device

Claims (10)

Requesting coordinate data of a plurality of measurement points set on the analysis object;
Calculating a mode characteristic based on measurement data at each measurement point obtained when an excitation test is performed to perform excitation at at least one point on the analysis target;
Based on the coordinate data, calculating a constraint equation between components constituting a characteristic matrix representing the vibration of the analysis object;
A step of identifying the characteristic matrix using the constraint condition expression based on the value of the mode characteristic calculated based on the measurement data,
The degree of freedom of the characteristic matrix is larger than the degree of freedom of measurement actually measured at the measurement point in the vibration test, and thus the number of mode characteristic values calculated based on the measurement data is greater than the above. Less than the number of modal characteristic values calculated based on the characteristic matrix,
In the step of identifying the characteristic matrix, the mode calculated based on the characteristic matrix only for the mode characteristic corresponding to the mode characteristic calculated based on the measurement data among the mode characteristics calculated based on the characteristic matrix A vibration analysis method for identifying the characteristic matrix by changing the value of each component of the characteristic matrix so that the value of the characteristic matches the value of the mode characteristic calculated based on the measurement data. Requesting coordinate data of a plurality of actual measurement points and virtual measurement points set on the analysis object;
A step of performing measurement only at an actual measurement point when an excitation test for exciting at least one point on the analysis object is performed, and calculating a mode characteristic based on measurement data at each actual measurement point obtained by excitation; ,
Calculating a constraint equation between components constituting a characteristic matrix representing the vibration of the analysis object based on the coordinate data;
A step of identifying the characteristic matrix using the constraint equation based on the calculated mode characteristic value based on the measurement data,
The above characteristic matrix is created with degrees of freedom equal to the degrees of freedom at both the real and virtual measurement points,
In the step of identifying the characteristic matrix, the mode calculated based on the characteristic matrix only for the mode characteristic corresponding to the mode characteristic calculated based on the measurement data among the mode characteristics calculated based on the characteristic matrix A vibration analysis method for identifying the characteristic matrix by changing the value of each component of the characteristic matrix so that the value of the characteristic matches the value of the mode characteristic calculated based on the measurement data.   The vibration analysis method according to claim 2, wherein the virtual measurement point is provided such that a measurement point closest to the virtual measurement point is an actual measurement point. Further comprising inputting a structural coupling state between the measurement points;
The vibration analysis method according to claim 2 or 3, wherein, in the step of calculating the constraint conditional expression, the constraint conditional expression is calculated based on a structural coupling state between measurement points in addition to the coordinate data. The characteristic matrix identified in the step of identifying the characteristic matrix is a mass matrix,
The vibration analysis method according to claim 2, further comprising a step of calculating a rigid body characteristic based on the mass matrix identified in the step of identifying the characteristic matrix. Measured from a measuring device that measures the displacement of the position, velocity, or acceleration at a measurement point on the analysis object, an excitation device that excites at least one point on the analysis object, and the measurement device and the excitation device. In a vibration analysis apparatus comprising an analysis unit that receives data and excitation data and performs arithmetic processing based on these data,
The analysis unit calculates a mode characteristic based on the measurement data at each measurement point obtained by the vibration test with the vibrator, and represents the vibration of the analysis object based on the coordinate data of each measurement point. Calculating a constraint expression between components constituting the characteristic matrix, identifying the characteristic matrix using the constraint expression based on the value of the mode characteristic calculated based on the measurement data,
The degree of freedom of the characteristic matrix is larger than the degree of freedom of measurement actually measured at the measurement point in the vibration test, and thus the number of mode characteristic values calculated based on the measurement data is greater than the above. Less than the number of modal characteristic values calculated based on the characteristic matrix,
When identifying the characteristic matrix, only the mode characteristic corresponding to the mode characteristic calculated based on the measurement data among the mode characteristics calculated based on the characteristic matrix is calculated based on the characteristic matrix. A vibration analysis device that identifies the characteristic matrix by changing the value of each component of the characteristic matrix so that the value of the mode characteristic matches the value of the mode characteristic calculated based on the measurement data. A measuring instrument that measures the displacement of the position, velocity, or acceleration at the actual measurement point among the actual measurement point and the virtual measurement point set on the analysis object, and a vibration exciter that excites at least one point on the analysis object. In the vibration analysis apparatus comprising an analysis unit that inputs measurement data and excitation data from the measuring instrument and the vibrator, respectively, and performs arithmetic processing based on these data,
The analysis unit calculates a mode characteristic based on measurement data at an actual measurement point obtained by an excitation test using the vibrator, and constitutes a characteristic matrix based on the coordinate data of the actual measurement point and the virtual measurement point And calculating the constraint equation between them, identifying the characteristic matrix using the constraint equation based on the value of the mode characteristic calculated based on the measurement data,
The above characteristic matrix is created with degrees of freedom equal to the degrees of freedom at both the real and virtual measurement points,
When identifying the characteristic matrix, only the mode characteristic corresponding to the mode characteristic calculated based on the measurement data among the mode characteristics calculated based on the characteristic matrix is calculated based on the characteristic matrix. A vibration analysis device that identifies the characteristic matrix by changing the value of each component of the characteristic matrix so that the value of the mode characteristic matches the value of the mode characteristic calculated based on the measurement data.   The vibration analysis apparatus according to claim 7, wherein the virtual measurement point is provided such that a measurement point closest to the virtual measurement point is an actual measurement point.   The vibration analysis apparatus according to claim 7 or 8, wherein when calculating the constraint condition formula, the constraint condition formula is calculated based on a structural coupling state between measurement points in addition to the coordinate data. The characteristic matrix identified above is a mass matrix,
The vibration analysis apparatus according to claim 7, wherein the analysis unit calculates a rigid body property based on the mass matrix identified in the step of identifying the property matrix.

Images