CN104536941B - A kind of frequency domain load recognition method based on Tikhonov regularizations - Google Patents

Abstract

The invention discloses a kind of frequency domain load recognition method based on Tikhonov regularizations.Frequency response function of this method for matrix inversion method in the presence of practical engineering application obtains the ill-conditioning problem in problem and calculating, frequency response function is obtained using finite element method, and frequency of use receptance function Matrix condition number evaluates the pathosis of plant equipment, load is identified based on different regularization methods when pathosis are different.And suitable ill-condition number threshold value is determined.The beneficial effects of the invention are as follows：Ill-conditioning problem in by the method matrix inversion method being overcome to be difficult to obtain plant equipment frequency response function and calculate in practical engineering application, improves the load identification precision in frequency domain, with good engineering application value.

Description

A kind of frequency domain load recognition method based on Tikhonov regularizations

Technical field

It is more particularly to a kind of to be based on Tikhonov regularizations the invention belongs to the identification of mechanical structure load and vibration analysis field Frequency domain load recognition method.

Background technology

With the development of modern project technology, it is more and more interesting issue that various vibration problems turn into engineering circles.And shake The identification in dynamic source, that is, our identification problems of load for saying, are again the key points of this problem.People are to acting on machine Load degree of concern in tool structure is increasing, it is studied gradually deeply.In real work, act in mechanical structure Influence of the dynamic loading to structure is very big, and with destructive and Unpredictability, aerodynamic loading that aircraft such as in-flight is subject to, Wave stroke that steamer is subject to etc..Therefore, the determination of dynamic loading is significant with research in the analysis of mechanical structure, exactly Identification dynamic loading is to ensure that the important guarantee of engineering structure reliability and security.

For most of engineering structures, the dynamic loading suffered by it is often to be difficult direct measurement, even impossible Measurement.Identification technology based on structure actual measurement response inverting dynamic loading is to determine the degradation pathways of dynamic loading, i.e. dynamic loading is known Other technology.So-called Dynamic Load Identification technology is vibrated according to the dynamic characteristic parameter of system and the dynamic system response reverse of actual measurement Dynamic load suffered by system.The theoretical method comparative maturity of current load identification, but can be used to instruct the knowledge of engineering practice Other method is actually rare, and its main cause is the presence of following application difficult point：On the one hand, due to research object, experimental field environment With the difference of point layout, can all there be certain change, the meter of load identification in the optimal measuring point number of participation load identification and position Easily there is ill phenomenon in calculation process, so as to cause recognition result insincere；On the other hand, in face of the knot of some complexity in practice During structure, the cost of actual test be it is very high, sometimes even not possible with.

The content of the invention

According to the technical problem present on, the invention discloses a kind of frequency domain load based on Tikhonov regularization methods Lotus recognition methods, the described method comprises the following steps：

S100, mould measurement is carried out to plant equipment, calculate damping matrix, so as to obtain the FEM model containing damping；

S200, use acceleration transducer collection machinery vibration equipment response signal, and using containing damping finite element mould Type calculates required frequency response function matrix；

The frequency response function matrix that S300, basis are obtained, calculates the Matrix condition number under different frequency；

S400, the Tikhonov regularization sides that different regularization parameters are chosen according to the Matrix condition number under different frequency Method carries out load identification.

The method disclosed in the present has the characteristics that：

(1) damping matrix is calculated using mould measurement, frequency response function is calculated by the FEM model containing damping, so that gram The drawbacks of taking complex mechanical equipment inconvenience or even test cannot be dismantled.

(2) pathosis of plant equipment are evaluated using frequency response function Matrix condition number, in plant equipment pathosis not Load is identified based on different regularization methods simultaneously, and suitable ill-condition number threshold value is determined.

Brief description of the drawings

Fig. 1 is the experimental bench of the analog mechanical equipment real work situation built, wherein A₁-A₇Represent acceleration transducer, F₁, F₂It is the unit sinusoidal force for applying；

Fig. 2 is that experimental bench external applied load recognition result figure abscissa is frequency, and unit is Hz, and ordinate is amplitude, and unit is G, danger signal is impact point vibration signals measured, and blue signal is load recognition result.

Specific embodiment

It is an object of the invention to provide a kind of frequency domain load recognition method based on Tikhonov regularization methods, for External applied load suffered by identification plant equipment is calculated in frequency domain.The present invention obtains structure frequency receptance function using finite element method, and Frequency of use receptance function Matrix condition number evaluates the pathosis of plant equipment, when pathosis are different based on different regularizations Method is identified to load.And suitable plant equipment ill-condition number threshold value is determined.The method can overcome Matrix Calculating Inverse method is difficult to the ill-conditioning problem in obtaining plant equipment frequency response function and calculating in practical engineering application, improves in frequency domain Interior structural loads accuracy of identification.

Present disclosure is described in further detail below in conjunction with the accompanying drawings：

Shown in reference picture 1, the experimental bench of analog mechanical equipment real work situation is built.It is the shell of cantilever in figure shown in 1 Install two vibrators additional to swash, it is worked with the exciting force of 160Hz, two vibration sources when being worked as equipment, both two carried outward Lotus (F₁, F₂).And gather signal with force snesor.Arrange that 7 acceleration transducers carry out the collection of vibration response signal on shell (A₁, A₂, A₃, A₄, A₅, A₆, A₇)。

Shown in reference picture 2, abscissa is frequency, and unit is Hz, and ordinate is amplitude, and unit is g, and danger signal is target Point vibration signals measured, blue signal is load recognition result.

The present invention is implemented according to the following steps：

S100, mould measurement is carried out to plant equipment, calculate damping matrix, so as to obtain the FEM model containing damping；

S200, use acceleration transducer collection machinery vibration equipment response signal, and using containing damping finite element mould Type calculates required frequency response function matrix；

The frequency response function matrix that S300, basis are obtained, calculates the Matrix condition number under different frequency；

S400, the Tikhonov regularization sides that different regularization parameters are chosen according to the Matrix condition number under different frequency Method carries out load identification.

Damping matrix is calculated in the step S100 to be specially：

S101, according to mode test result, by following equation calculating machine equipment FEM model proportional damping parameter；

In above formula, α is the Tuned mass damper coefficient of the plant equipment；β is the rigidity quality coefficient of the plant equipment；ξ₁ And ξ₂The ratio between the actual damping obtained by the mould measurement and critical damping；ω₁And ω₂It is different moulds in the mould measurement Natural angular frequency under state；

S102, according to damping ratio parameter calculate damping matrix；

[C]=α [M]+β [K]

In above formula：[C] is damping matrix；[M] is the gross mass matrix of plant equipment；[K] is the global stiffness of plant equipment Matrix.

[M] is automatically generated by the material properties of plant equipment with [K] by ANSYS softwares.

The step S200 is specially：

S201, arrangement acceleration transducer and selection energized position in plant equipment；

S202, the vibration response signal using each measuring point of acceleration transducer collection machinery equipment；

S203, the frequency response function matrix using each measuring point of FEM model calculating containing damping to energized position.

The step S203 is specially：

Corresponding energized position applies unit sinusoidal force on the FEM model containing damping, to the finite element containing damping Model carries out harmonic responding analysis according to the following equation：

(-ω[M]+iω[C]+[K])({u₁}+i{u₂)={ F₁}+i{F₂}

In above formula：ω is angular frequency, is set as needed；{u₁It is the real displacement at point position；{u₂It is point position The virtual displacement at place；{F₁It is the exciting force real part of energized position applying；{F₂It is the exciting force imaginary part of energized position applying；

When applying unit sinusoidal force, { F in above formula₁It is 1, { F₂For 0 when, { the u being calculated by above formula₁}+i{u₂Be It is the displacement frequency response function of measuring point to energized position；

ω²({u₁}+i{u₂) measuring point to the acceleration frequence responses function of energized position is, the frequency for as needing is rung Answer function.

The frequency response function component frequency receptance function matrix of the different measuring points position；The frequency response function square The element of the row of m rows n-th, as m-th point position to n-th frequency response function of energized position in battle array.

The step S300 is specially：

Matrix condition number under different frequency is calculated according to below equation；

In above formula：n_iIt is the Matrix condition number under frequency i；H_iIt is the frequency response function matrix under frequency i；For under frequency i Frequency response function matrix associate matrix.

Described use sentences method for distinguishing based on conditional number, suitable according to different frequency response function matrix conditional number selections Identification of the regularization method to load；And the ill-condition number threshold value that is given is as follows：

Mechanical structure load is solved using direct matrix in Practical Project：

In formula：

{F}_n×1--- external applied load column vector；

--- the inverse matrix of frequency response function matrix；

{X}_m×1--- work response column vector.

M --- work response signal quantity；

N --- external applied load quantity.

Load identification in Practical Project is often ill posed, the plant equipment frequency response function matrix [H] of morbid state_m×n During directly inverting, less measurement error in response signal can seriously be amplified so that the load that reverse is obtained is serious Deviate actual, recognizing load becomes meaningless.In order to solve the problems, such as it is this caused by plant equipment morbid state, starting carry Before lotus calculates, Matrix condition number is calculated first, obtain conditional number of the plant equipment on each Frequency point, be i.e. disease State situation；Secondly, in assumed (specified) load, the calculating that the Frequency point is entered on each Frequency point is circulated：

Step S400 is specially：

S401, selection Matrix condition number threshold value are 1000, during conditional number ＞ 1000, plant equipment Very Ill-conditioned, conditional number When≤1000, plant equipment is slightly ill；

S402, when conditional number is more than 1000, during selection normal crossing check addition determines Tikhonov regularization methods Optimal regularization parameter, conversely, selection L-curve method carries out the selection of optimal regularization parameter；

S403, based on regularization parameter, the load of the Frequency point is calculated using Tikhonov regularization methods, enter afterwards The load identification of next Frequency point.

The step S402 is specially：

If regularization parameter is λ

Normal crossing check addition：

In above formula：| | | | it is Euclidean norms；N counts out for response；X is the vibratory response of measurement；H is frequency response Function is called；F is load to be identified；I is unit matrix；C (λ)=H (H^HH+λI)^-1H^H；B (λ) is diagonal matrix, and diagonal matrix is by 1/ (1-C_kk(λ)) try to achieve, C_kk(λ) is the diagonal item of Matrix C (λ)；

Work as V₀λ when (λ) obtains minimum value is the optimal regularization parameter of normal crossing check addition determination；

L-curve method：

If ρ (λ)=| | HF-X | |, η (λ)=| | F | |

In above formula：X is the vibratory response of measurement；H is called for frequency response function；F is load to be identified；

When K (λ) takes maximum, optimal regularization parameter λ is determined by the maximum point.

Embodiment 1：

According to the experimental bench of the real work situation of analog mechanical equipment shown in Fig. 1.

1. pair structure carries out hammering mould measurement, and sets up the FEM model of shell structure.By with finite element modal Analysis result is contrasted, correction model boundary condition, and is calculated structure ratio by mode test result and formula (2), (3) Example damping parameter α and β, obtain the FEM model containing damping.

2. shown in reference picture 1, acceleration transducer is arranged in identified response point, gather each measuring point in vibrator Vibration response signal under 160Hz sine excitation power.And calculate each measuring point to vibrator position using the FEM model containing damping The frequency response function matrix put.

3., according to the plant equipment frequency response function matrix for obtaining, the Matrix condition number under different frequency is calculated.

4., according to the Matrix condition number for obtaining, on each Frequency point, when conditional number is more than 1000, selection is common to hand over Fork method of calibration determines regularization parameter；Conversely, selection L-curve method carries out the selection of optimal regularization parameter.It is determined that just Then change after parameter, the load of the Frequency point is calculated using Tikhonov regularization methods.Be identified result as shown in Fig. 2 Recognition effect is preferable over the entire frequency band.

Above to method provided by the present invention, it is described in detail, specific case used herein is to the present invention Principle and implementation method be set forth, the explanation of above example is only intended to help and understands the method for the present invention and its core Thought is thought；Simultaneously for those of ordinary skill in the art, according to thought of the invention, in specific embodiment and model is applied Place and will change, in sum, this specification content should not be construed as limiting the invention.

Claims (8)

1. a kind of frequency domain load recognition method based on Tikhonov regularizations, it is characterised in that：Methods described includes following step Suddenly：

S100, mould measurement is carried out to plant equipment, calculate damping matrix, so as to obtain the FEM model containing damping, wherein, Damping matrix is calculated to be specially：

S101, according to mode test result, by following equation calculating machine equipment FEM model proportional damping parameter；

$\mathrm{\α}=2{\mathrm{\ω}}_1{\mathrm{\ω}}_2({\mathrm{\ξ}}_1{\mathrm{\ω}}_2-{\mathrm{\ξ}}_2{\mathrm{\ω}}_1)/({\mathrm{\ω}}_2^2-{\mathrm{\ω}}_1^2)$

$\mathrm{\β}=2({\mathrm{\ξ}}_2{\mathrm{\ω}}_2-{\mathrm{\ξ}}_1{\mathrm{\ω}}_1)/\left({\mathrm{\ω}}_2^2-{\mathrm{\ω}}_1^2\right)$

In above formula, α is the Tuned mass damper coefficient of the plant equipment；β is the rigidity quality coefficient of the plant equipment；ξ₁And ξ₂ The ratio between the actual damping obtained by the mould measurement and critical damping；ω₁And ω₂For under different modalities in the mould measurement Natural angular frequency；

S102, according to damping ratio parameter calculate damping matrix；

[C]=α [M]+β [K]

In above formula：[C] is damping matrix；[M] is the gross mass matrix of plant equipment；Area] it is the global stiffness matrix of plant equipment；

S200, use acceleration transducer collection machinery vibration equipment response signal, and using containing damping FEM model meter Frequency response function matrix needed for calculating；

The frequency response function matrix that S300, basis are obtained, calculates the Matrix condition number under different frequency；

S400, the Tikhonov regularization methods for choosing different regularization parameters according to the Matrix condition number under different frequency enter Row load is recognized.

2. method according to claim 1, it is characterised in that the step S200 is specially：

S201, in plant equipment arrange acceleration transducer determine measuring point and selection energized position；

S202, the vibration response signal using each measuring point of acceleration transducer collection machinery equipment；

S203, the frequency response function matrix using each measuring point of FEM model calculating containing damping to energized position.

3. method according to claim 2, it is characterised in that step S203 is specially：

Corresponding energized position applies unit sinusoidal force on the FEM model containing damping, to the FEM model containing damping Harmonic responding analysis are carried out according to the following equation：

(-ω[M]+iω[C]+[K])({u₁}+i{u₂)={ F₁}+i{F₂}

In above formula：ω is angular frequency, is set as needed；{u₁It is the real displacement at point position；{u₂For at point position Virtual displacement；{F₁It is the exciting force real part of energized position applying；{F₂It is the exciting force imaginary part of energized position applying；

When applying unit sinusoidal force, { F in above formula₁Be₁, { F₂For 0 when, { the u being calculated by above formula₁}+i{u₂It is measuring point To the displacement frequency response function of energized position；

ω²({u₁}+i{u₂) measuring point to the acceleration frequence responses function of energized position is, the frequency response letter for as needing Number.

4. method according to claim 3, it is characterised in that the frequency response function composition frequency of the different measuring points position Rate receptance function matrix；The element of the row of m rows n-th, as m-th point position to n-th in the frequency response function matrix The frequency response function of individual energized position.

5. method according to claim 4, it is characterised in that the step S300 is specially：

Matrix condition number under different frequency is calculated according to below equation；

n_i=H_i ^HH_i

In above formula：n_iIt is the Matrix condition number under frequency i；H_iIt is the frequency response function matrix under frequency i；It is the frequency under frequency i Ring the associate matrix of Jacobian matrix.

6. method according to claim 5, it is characterised in that the step S400 is specially：

S401, selection Matrix condition number threshold value are 1000, during conditional number ＞ 1000, plant equipment Very Ill-conditioned, conditional number≤ When 1000, plant equipment is slightly ill；

S402, when conditional number is more than 1000, selection normal crossing check addition determines optimal in Tikhonov regularization methods Regularization parameter, conversely, selection L-curve method carries out the selection of optimal regularization parameter；

S403, based on regularization parameter, the load of the Frequency point is calculated using Tikhonov regularization methods, afterwards into next The load identification of Frequency point.

7. method according to claim 1, it is characterised in that：[M] and [K] are passed through by the material properties of plant equipment ANSYS softwares are automatically generated.

8. method according to claim 6, it is characterised in that：The step S402 is specially：

If regularization parameter is λ

Normal crossing check addition：

$V_0\left(\mathrm{\λ}\right)=\frac{1}{n}\left|\right|X-HF|{|}^2=\frac{1}{n}|\left|B\right(\mathrm{\λ}\left)\right(I-C\left(\mathrm{\λ}\right)\left)X\right|{|}^2$

In above formula：| | | | it is Euclidean norms；N counts out for response；X is the vibratory response of measurement；H is frequency response function It is called；F is load to be identified；I is unit matrix；C (λ)=H (H^HH+λI)^-1H^H；B (λ) is diagonal matrix, and diagonal matrix is by 1/ (1- C_kk(λ)) try to achieve, C_kk(λ) is the diagonal item of Matrix C (λ)；

Work as V₀λ when (λ) obtains minimum value is the optimal regularization parameter of normal crossing check addition determination；

L-curve method：

If ρ (λ)=| | HF-X | |, η (λ)=| | F | |

In above formula：X is the vibratory response of measurement；H is called for frequency response function；F is load to be identified；

$K\left(\mathrm{\λ}\right)=\frac{|{\mathrm{\ρ}}^{\′}{\mathrm{\η}}^{\′\′}-{\mathrm{\ρ}}^{\′\′}{\mathrm{\η}}^{\′}|}{{({\mathrm{\ρ}}^{\′2}+{\mathrm{\η}}^{\′2})}^{3/2}}$

When K (λ) takes maximum, optimal regularization parameter λ is determined by the maximum point.