CN104536941A - Frequency domain load identification method based on Tikhonov regularization - Google Patents

Abstract

The invention discloses a frequency domain load identification method based on Tikhonov regularization. The method aims at solving the problems existing in frequency response function acquisition in practical engineering application of a matrix inversion method and the morbidity problems existing in calculation. A frequency response function is acquired through a finite element method, the morbidity of mechanical equipment is evaluated through the number of matrix conditions of the frequency response function, and loads are identified based on different regularization methods when the morbidity is different. A proper mobility condition number threshold value is determined. The method has the advantages that the problem that the frequency response function of the mechanical equipment is difficult to acquire in practical engineering application of the matrix inversion method and the morbidity problems existing in calculation can be solved, the load identification accuracy in a frequency domain is improved, and high engineering application value is achieved.

Description

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

Technical field

The invention belongs to the identification of physical construction load and vibration analysis field, particularly a kind of frequency domain load recognition method based on Tikhonov regularization.

Background technology

Along with the development of modern project technology, various vibration problem becomes the more and more concerned problem of engineering circles.And the identification of vibration source, namely our identification problem of load of saying, be again the key point of this problem.People are increasing to the load degree of concern acted in physical construction, to its research gradually deeply.In real work, act on the dynamic loading in physical construction very large on the impact of structure, and there is destructiveness and Unpredictability, the wave stroke etc. that the aerodynamic loading be subject to as aloft aircraft, steamer are subject to.Therefore, significant in the analysis that dynamic loading fixes on physical construction really and research, identify that dynamic loading guarantees the important guarantee of engineering structure reliability and security exactly.

For most of engineering structure, the dynamic loading suffered by it is difficult to directly measurement often, or even can not measure.The recognition technology of structure based actual measurement response inverting dynamic loading determines the degradation pathways of dynamic loading, i.e. Dynamic Load Identification technology.So-called Dynamic Load Identification technology is the dynamic load suffered by the dynamic characteristic parameter of system and the dynamic system response reverse vibrational system of actual measurement.The theoretical method comparative maturity of current load identification, but the recognition methods that can be used to guide engineering practice is actually rare, its main cause applies difficult point below existing: on the one hand, due to the difference of research object, experimental field environment and point layout, all can there be certain change the best measuring point number and the position that participate in load identification, very easily there is ill phenomenon in the computation process of load identification, thus causes recognition result insincere; On the other hand, in the reality during structure of some complexity, the cost of actual test is very high, sometimes or even impossible.

Summary of the invention

According to the above technical matters existed, the invention discloses a kind of frequency domain load recognition method based on Tikhonov regularization method, said method comprising the steps of:

S100, mould measurement is carried out to plant equipment, calculate damping matrix, thus obtain the finite element model containing damping;

S200, use acceleration transducer collection machinery vibration equipment response signal, and utilize containing the frequency response function matrix needed for the finite element model calculating of damping;

The frequency response function matrix that S300, basis obtain, calculates the Matrix condition number under different frequency;

S400, the Tikhonov regularization method choosing different regularization parameters according to the Matrix condition number under different frequency carry out load identification.

The method disclosed in the present has following characteristics:

(1) utilize mould measurement to calculate damping matrix, calculate frequency response function by the finite element model containing damping, thus overcome the drawback that complex mechanical equipment inconvenience even cannot dismantle test.

(2) utilize frequency response function Matrix condition number to evaluate the pathosis of plant equipment, based on different regularization methods, load is identified when plant equipment pathosis is different, and determine suitable ill-condition number threshold value.

Accompanying drawing explanation

Fig. 1 is the experiment table of the analog mechanical equipment real work situation of building, wherein A ₁-A ₇represent acceleration transducer, F ₁, F ₂for the unit sinusoidal force applied;

Fig. 2 is experiment table external applied load recognition result figure horizontal ordinate is frequency, and unit is Hz, and ordinate is amplitude, and unit is g, and danger signal is impact point vibration signals measured, and blue signal is load recognition result.

Embodiment

The object of this invention is to provide a kind of frequency domain load recognition method based on Tikhonov regularization method, identifying external applied load suffered by plant equipment for calculating in frequency domain.The present invention utilizes Finite Element Method to obtain structure frequency response function, and frequency of utilization response function Matrix condition number evaluates the pathosis of plant equipment, identifies based on different regularization methods when pathosis is different to load.And determine suitable plant equipment ill-condition number threshold value.The method can overcome matrix inversion method and be difficult to obtain plant equipment frequency response function and the ill-conditioning problem in calculating in practical engineering application, improves the structural loads accuracy of identification in frequency domain.

Below in conjunction with accompanying drawing, content of the present invention is described in further detail:

With reference to shown in Fig. 1, build the experiment table of analog mechanical equipment real work situation.In figure shown in 1, the shell for cantilever installs two vibrators additional and swashs, and makes it with the exciting force work of 160Hz, two vibration sources when working as equipment, both two external applied load (F ₁, F ₂).And use force snesor collection signal.Shell being arranged, 7 acceleration transducers carry out the collection (A of vibration response signal ₁, A ₂, A ₃, A ₄, A ₅, A ₆, A ₇).

With reference to shown in Fig. 2, horizontal ordinate is frequency, and unit is Hz, and ordinate is amplitude, and unit is g, and danger signal is impact point vibration signals measured, and blue signal is load recognition result.

The present invention implements according to the following steps:

S100, mould measurement is carried out to plant equipment, calculate damping matrix, thus obtain the finite element model containing damping;

S200, use acceleration transducer collection machinery vibration equipment response signal, and utilize containing the frequency response function matrix needed for the finite element model calculating of damping;

The frequency response function matrix that S300, basis obtain, calculates the Matrix condition number under different frequency;

S400, the Tikhonov regularization method choosing different regularization parameters according to the Matrix condition number under different frequency carry out load identification.

Calculate damping matrix in described step S100 to be specially:

S101, according to mode test result, by following formulae discovery machinery equipment finite element 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)/({\mathrm{\ω}}_2^2-{\mathrm{\ω}}_1^2)$

In above formula, α is the Tuned mass damper coefficient of described plant equipment; β is the rigidity quality coefficient of described plant equipment; ξ ₁and ξ ₂the ratio of the actual damping obtained by described mould measurement and critical damping; ω ₁and ω ₂for the natural angular frequency under different modalities in described mould measurement;

S102, according to damping ratio parameter calculate damping matrix;

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

In above formula: [C] is damping matrix; The gross mass matrix that [M] is plant equipment; The global stiffness matrix that [K] is plant equipment.

Described [M] and [K] are generated by ANSYS software automatically by the material properties of plant equipment.

Described step S200 is specially:

S201, arrange in plant equipment acceleration transducer and select energized position;

S202, utilize the vibration response signal of each measuring point of acceleration transducer collection machinery equipment;

S203, utilization calculate the frequency response function matrix of each measuring point to energized position containing the finite element model of damping.

Described step S203 is specially:

Applying unit sinusoidal force containing energized position corresponding on the finite element model of damping, according to the following equation harmonic responding analysis carried out to the finite element model containing damping:

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

In above formula: ω is angular frequency, sets as required; { u ₁it is the real displacement at point position place; { u ₂it is the virtual displacement at point position place; { F ₁for energized position apply exciting force real part; { F ₂for energized position apply exciting force imaginary part;

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

ω ²({ u ₁+ i{u ₂) be the acceleration frequence responses function of measuring point to energized position, be the frequency response function of needs.

The frequency response function component frequency response function matrix of described different measuring points position; In described frequency response function matrix, the element of capable n-th row of m, is the frequency response function of m point position to the n-th energized position.

Described step S300 is specially:

According to the Matrix condition number under following formulae discovery different frequency;

$n_i=H_i^H{H}_i$

In above formula: n _ifor the Matrix condition number under frequency i; H _ifor the frequency response function matrix under frequency i; for the associate matrix of the frequency response function matrix under frequency i.

Described use sentences method for distinguishing based on conditional number, selects suitable regularization method to the identification of load according to different frequency response function matrix conditional numbers; And the ill-condition number threshold value provided is as follows:

Direct matrix is utilized to solve physical construction load in Practical Project:

${\left\{F\right\}}_{n\×1}={\left[H\right]}_{m\×n}^{-1}{\left\{X\right\}}_{m\×1}$

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 ill posed often, the plant equipment frequency response function matrix [H] of morbid state _{m × n}in the process of directly inverting, measuring error less in response signal seriously can be amplified, the load substantial deviation that reverse is obtained is actual, makes load identification become meaningless.In order to solve this problem caused by plant equipment morbid state, before beginning LOAD FOR, first Matrix condition number being calculated, obtaining the conditional number of plant equipment on each Frequency point, i.e. pathological situation; Secondly, when assumed (specified) load, each Frequency point enters the computation cycles of this Frequency point:

Step S400 is specially:

S401, to choose Matrix condition number threshold value be 1000, during conditional number > 1000, and plant equipment Very Ill-conditioned, during conditional number≤1000, plant equipment is slightly ill;

S402, when conditional number is greater than 1000, select the optimum regularization parameter that normal crossing check addition is determined in Tikhonov regularization method, otherwise, select L curve method to carry out choosing of optimum regularization parameter;

S403, based on regularization parameter, use Tikhonov regularization method to calculate the load of this Frequency point, enter the load identification of next Frequency point afterwards.

Described step S402 is specially:

If regularization parameter is λ

Normal crossing check addition:

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

In above formula: || || be Euclidean norm; N is response point number; X is the vibratory response of measuring; H is that 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(λ) be the diagonal angle item of Matrix C (λ);

Work as V ₀(λ) λ when obtaining minimum value is the optimum regularization parameter that normal crossing check addition is determined;

L curve method:

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

In above formula: X is the vibratory response of measuring; H is that frequency response function is called; F is load to be identified;

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

When K (λ) gets maximal value, determine optimum regularization parameter λ by this maximum point.

Embodiment 1:

According to the experiment table of the equipment of analog mechanical shown in Fig. 1 real work situation.

1. pair this structure carries out hammering mould measurement, and sets up the finite element model of shell structure.By contrasting with finite element modal analysis result, correction model boundary condition, and calculate structure proportion damping parameter α and β by mode test result and formula (2), (3), obtain the finite element model containing damping.

2., with reference to shown in Fig. 1, determined response point arranges acceleration transducer, gathers the vibration response signal of each measuring point under vibrator 160Hz sine excitation power.And utilize the finite element model containing damping to calculate the frequency response function matrix of each measuring point to vibrator position.

3., according to the plant equipment frequency response function matrix obtained, calculate the Matrix condition number under different frequency.

4., according to the Matrix condition number obtained, on each Frequency point, when conditional number is greater than 1000, select normal crossing method of calibration determination regularization parameter; Otherwise, select L curve method to carry out choosing of best regularization parameter.After determining regularization parameter, Tikhonov regularization method is used to calculate the load of this Frequency point.Obtain recognition result as shown in Figure 2, recognition effect is all better over the entire frequency band.

Above to method provided by the present invention, be described in detail, apply specific case herein and set forth principle of the present invention and embodiment, the explanation of above embodiment just understands method of the present invention and core concept thereof for helping; Meanwhile, for one of ordinary skill in the art, according to thought of the present invention, all will change in specific embodiments and applications, in sum, this description should not be construed as limitation of the present invention.

Claims (9)

1., based on a frequency domain load recognition method for Tikhonov regularization, it is characterized in that: said method comprising the steps of:

S100, mould measurement is carried out to plant equipment, calculate damping matrix, thus obtain the finite element model containing damping;

S200, use acceleration transducer collection machinery vibration equipment response signal, and utilize containing the frequency response function matrix needed for the finite element model calculating of damping;

The frequency response function matrix that S300, basis obtain, calculates the Matrix condition number under different frequency;

S400, the Tikhonov regularization method choosing different regularization parameters according to the Matrix condition number under different frequency carry out load identification.

2. method according to claim 1, is characterized in that, preferably, calculates damping matrix and be specially in described step S100:

S101, according to mode test result, by following formulae discovery machinery equipment finite element 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)/({\mathrm{\ω}}_2^2-{\mathrm{\ω}}_1^2)$

In above formula, α is the Tuned mass damper coefficient of described plant equipment; β is the rigidity quality coefficient of described plant equipment; ξ ₁and ξ ₂the ratio of the actual damping obtained by described mould measurement and critical damping; ω ₁and ω ₂for the natural angular frequency under different modalities in described mould measurement;

S102, according to damping ratio parameter calculate damping matrix;

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

In above formula: [C] is damping matrix; The gross mass matrix that [M] is plant equipment; The global stiffness matrix that [K] is plant equipment.

3. method according to claim 1, is characterized in that, described step S200 is specially:

S201, arrange in plant equipment acceleration transducer determination measuring point and select energized position;

S202, utilize the vibration response signal of each measuring point of acceleration transducer collection machinery equipment;

S203, utilization calculate the frequency response function matrix of each measuring point to energized position containing the finite element model of damping.

4. method according to claim 3, is characterized in that, step S203 is specially:

Applying unit sinusoidal force containing energized position corresponding on the finite element model of damping, according to the following equation harmonic responding analysis carried out to the finite element model containing damping:

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

In above formula: ω is angular frequency, sets as required; { u ₁it is the real displacement at point position place; { u ₂it is the virtual displacement at point position place; { F ₁for energized position apply exciting force real part; { F ₂for energized position apply exciting force imaginary part;

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

ω ²({ u ₁+ i{u ₂) be the acceleration frequence responses function of measuring point to energized position, be the frequency response function of needs.

5. method according to claim 4, is characterized in that, the frequency response function component frequency response function matrix of described different measuring points position; In described frequency response function matrix, the element of capable n-th row of m, is the frequency response function of m point position to the n-th energized position.

6. method according to claim 5, is characterized in that, described step S300 is specially:

According to the Matrix condition number under following formulae discovery different frequency;

n _i＝H _i ^HH _i

In above formula: n _ifor the Matrix condition number under frequency i; H _ifor the frequency response function matrix under frequency i; for the associate matrix of the frequency response function matrix under frequency i.

7. method according to claim 6, is characterized in that, described step S400 is specially:

S401, to choose Matrix condition number threshold value be 1000, during conditional number > 1000, and plant equipment Very Ill-conditioned, during conditional number≤1000, plant equipment is slightly ill;

S402, when conditional number is greater than 1000, select the optimum regularization parameter that normal crossing check addition is determined in Tikhonov regularization method, otherwise, select L curve method to carry out choosing of optimum regularization parameter;

S403, based on regularization parameter, use Tikhonov regularization method to calculate the load of this Frequency point, enter the load identification of next Frequency point afterwards.

8. method according to claim 2, is characterized in that: described [M] and [K] are generated by ANSYS software automatically by the material properties of plant equipment.

9. method according to claim 7, is characterized in that: described step S402 is specially:

If regularization parameter is λ

Normal crossing check addition:

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

In above formula: || || be Euclidean norm; N is response point number; X is the vibratory response of measuring; H is that frequency response function is called; F is load to be identified; I is unit matrix; for diagonal matrix, diagonal matrix is by 1/ (1-C _kk(λ)) try to achieve, C _kk(λ) be the diagonal angle item of Matrix C (λ);

Work as V ₀(λ) λ when obtaining minimum value is the optimum regularization parameter that normal crossing check addition is determined;

L curve method:

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

In above formula: X is the vibratory response of measuring; H is that frequency response function is called; F is load to be identified;

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

When K (λ) gets maximal value, determine optimum regularization parameter λ by this maximum point.