A kind of wind induced structural vibration based on auto-correlation function responds efficient frequency domain estimation method
Technical field
The present invention relates to the efficient frequency domain estimation method of wind vibration response based on auto-correlation function.
Background technology
Wind-induced response is large and complex structure (such as Long Span Roof Structures, high-rise building etc.) wind force proofing design
One of important step.The vibration characteristics of large and complex structure is difficult to represent by the simple vibration shape, the contribution of high order mode to be considered
And the coupling between the vibration shape, therefore, the structural dynamic response under wind loads effect is extremely complex.Traditional time domain wind shakes
Response analysis, needs to input a large amount of wind load time-history information, uses Newmark-β integration method, and amount of calculation is relatively big, required during calculating
Memory space is relatively big, causes computational efficiency low.And traditional frequency domain wind-induced response, need also exist for inputting bulk information, bag
Include the covariance matrix of wind load, use numerical integrating, calculate the longest, cause computational efficiency low.Additionally, traditional frequency domain
Computational methods are difficult to consider modes coupling effect, calculate error big.
Summary of the invention
The invention aims to solve that existing wind-induced response method computational efficiency is low and to calculate error big
Shortcoming, and propose a kind of wind induced structural vibration based on auto-correlation function and respond efficient frequency domain estimation method.
A kind of wind induced structural vibration based on auto-correlation function responds efficient frequency domain estimation method detailed process:
Step one: the body structure surface pulsating wind pressure time course data p recorded based on wind tunnel testiT (), calculates each wind load and adds
The meansigma methods of the blast time course data of loading pointThe standard deviation sigma of blast time course datapi, the auto-correlation function of blast time course data
Rpi(τ) and the coherent function Coh of blast time course datapij(ω);
Wherein, i, j are any two point in wind load load(ing) point sum Q, and 1≤i≤Q, 1≤j≤Q, Q are that wind load loads
Point sum, for positive integer;
Step 2: use exponential function matching auto-correlation function, obtains the dimensionless blast spectrum of the wind load load(ing) point of i point
The crest frequency ω of curvemi;
From the definition of auto-power spectrum
In formula, Spi(ω) it is the blast Power spectral density of each load(ing) point;ω is frequency;τ is the time difference;
Then dimensionless blast stave is shown asCan draw,
In formula, Si(ω) it is dimensionless blast spectrum;
Step 3: use exponential function matching coherent function, obtain the relevant index k of this structurec, calculate i, j's 2
Wind load coherence factor;
Step 4: extract the mass matrix [M] of structure, stiffness matrix [K], damping matrix [C], wind load load(ing) point and knot
The transition matrix [R] of structure degree of freedom, the face that affects matrix [I] of response;
Structure is carried out model analysis, solves characteristic equationObtain each rank natural frequency of vibration of structure
ωnk, k=1,2 ..., N, k are the rank number of mode corresponding to the natural frequency of vibration, and N is the Degree of Structure Freedom number, for positive integer;
Take structural eigenvector matrixWhereinFor kth, r first order mode, k,
R=1,2 ... W;W is to calculate selected rank number of mode, and W is positive integer, 10≤W≤N;K, r are natural frequency of vibration ωnk、ωnrInstitute
Corresponding rank number of mode;MakeThe broad sense damping ratio of computation structure k first order mode
In formula, T is transposed matrix;
Step 5: calculate Analytical Integration according to step one, step 2, step 3 and step 4;Obtain [Σkr];
Step 6: calculate modal response covarianceForm modal response covariance square
Battle array;
Step 7: calculate dynamic respond covariance matrix [Σx]=[Ψ]T[Σy] [Ψ], and then it is the loudest to try to achieve it
The covariance matrix answered.
The invention have the benefit that
This method designs the large and complex structure efficient frequency domain estimation method of wind vibration response based on auto-correlation function, solves to pass
The system wind-induced response method shortcoming that computational efficiency is low and computational accuracy is poor in large and complex structure is analyzed.
First the wind loads that wind tunnel test obtains is processed as the statistic for calculating and frequency spectrum parameter by the method
Information, enormously simplify the input of wind load.Then, in conjunction with modal analysis method, to each vibration shape, to wind on each degree of freedom
Load frequency spectrum carries out the Analytical Integration simplified, and obtains modal response covariance.Finally, by the combination of coupled modes is obtained wind
The covariance of vibration response.The method substantially increases computational efficiency while ensureing modal coupling computational accuracy.
The present invention is ensureing on the premise of computational accuracy, is greatly improved large and complex structure wind and shakes the efficiency calculated, fall
Low calculate shared by memory space.Take full advantage of the Computing Principle of Analytical Integration, numerical integration is converted into parsing long-pending
Point algebraic operation, it is to avoid the error that numerical integration is brought, also substantially increase computational efficiency simultaneously.In implementation process,
Step is clear and definite, and operability is relatively strong, is found by instance analysis, and the present invention relatively traditional algorithm efficiency significantly increases, more practical.
Use this invention that certain fan-shaped plan Long-span Cantilever grid structure is carried out wind-induced response, compared with traditional analysis,
Maximum error is less than 1%, is much better than traditional simplification frequency domain algorithm error (16.7%), and efficiency relatively traditional algorithm improves 200
Times.
Accompanying drawing explanation
Fig. 1 is the flow chart of the present invention;
Fig. 2 a isTake-6.42, σpiTake the body structure surface point wind load time-history of 2.48 by the knot of step one statistical modeling
Really schematic diagram, vertical coordinate isFor normalized wind load time-history, abscissa is t, for time (second), pi(t)
For wind load time-history, piT () is Wind Loads Acting, σpiFor wind load standard deviation;
Fig. 2 b is Rpi(τ)=exp (-ωmi| τ |), ωmiThe body structure surface point wind load time-history of=2.40rad/s is by step
The result schematic diagram of a rapid statistical modeling, vertical coordinate is Rpi(τ), Rpi(τ) it is the auto-correlation function of wind load, ωmiFor wind load
Frequency, abscissa be τ, τ be the time difference (second);
Fig. 2 c isBody structure surface point wind load time-history by the knot of step one statistical modeling
Really schematic diagram, vertical coordinate isFor dimensionless wind load power spectrum, abscissa be ω, ω be frequency (arc
Degrees second), Spi(ω) it is wind load power spectral density function;
Fig. 3 a is that fan-shaped plan is encorbelmented the FEM (finite element) model schematic diagram of grid structure;
Fig. 3 b is that fan-shaped plan grid structure of encorbelmenting carries out, by step 2, the first order mode (ω that model analysis providesn1=
7.4rad/s) result schematic diagram;
Fig. 3 c is that fan-shaped plan grid structure of encorbelmenting carries out, by step 2, the second_mode (ω that model analysis providesn2=
7.5rad/s) result schematic diagram;
Fig. 3 d is that fan-shaped plan grid structure of encorbelmenting carries out, by step 2, the second_mode (ω that model analysis providesn3=
9.9rad/s) result schematic diagram;
Fig. 4 a is the modal response root-mean-square result schematic diagram that step 3 calculates, and RMS modal disp is that modal response is equal
Root value, modal order is rank number of mode, and frequency is natural frequency of structures, and bar diagram is modal response standard deviation,
Circle is the frequency of structure each order mode state;
Fig. 4 b is the modal response covariance matrix result schematic diagram that step 4 calculates;
Fig. 5 is algorithm and the traditional algorithm comparison in precision and efficiency of the present invention.
Detailed description of the invention
Detailed description of the invention one: combine Fig. 1 and present embodiment is described, present embodiment a kind of based on auto-correlation function
Wind induced structural vibration responds efficient frequency domain estimation method detailed process:
Step one: (such as Long Span Roof Structures is (such as two supports steel structural roof for the structure recorded based on wind tunnel test
Span is more than 36m, the cantilever span structure more than 10m, such as most common stadium, departure hall, conference and exhibition center's exhibition
The roof structure in shop etc.), high-rise building (such as highly high-level structure more than 100m)) surface pulsating wind pressure time course data pi
T (), calculates the meansigma methods of the blast time course data of each wind load load(ing) pointThe standard deviation sigma of blast time course datapi, blast time
The auto-correlation function R of number of passes evidencepi(τ) and the coherent function Coh of blast time course datapij(ω);
Wherein, i, j are any two point in wind load load(ing) point sum Q, and 1≤i≤Q, 1≤j≤Q, Q are that wind load loads
Point sum, for positive integer;
Step 2: use exponential function matching auto-correlation function, obtains the dimensionless blast spectrum of the wind load load(ing) point of i point
The crest frequency ω of curvemi;
From the definition of auto-power spectrum
In formula, Spi(ω) it is the blast Power spectral density of each load(ing) point;ω is frequency;τ is the time difference;
Then dimensionless blast stave is shown asCan draw,
In formula, Si(ω) it is dimensionless blast spectrum;
Step 3: use exponential function matching coherent function, obtain the relevant index k of this structurec, calculate i, j's 2
Wind load coherence factor;
Step 4: extract the mass matrix [M] of structure, stiffness matrix [K], damping matrix [C], wind load load(ing) point and knot
The transition matrix [R] of structure degree of freedom, the face that affects matrix [I] of response;
Structure is carried out model analysis, solves characteristic equationObtain each rank natural frequency of vibration of structure
ωnk, k=1,2 ..., N, k are the rank number of mode corresponding to the natural frequency of vibration, and N is the Degree of Structure Freedom number, for positive integer;
Take structural eigenvector matrixWhereinFor kth, r first order mode, k,
R=1,2 ... W;W is to calculate selected rank number of mode, and W is positive integer, determines according to engineering experience, generally, and 10≤W≤N;
K, r are natural frequency of vibration ωnk、ωnrCorresponding rank number of mode;MakeComputation structure k rank shake
The broad sense damping ratio of type
In formula, T is transposed matrix;
Step 5: calculate Analytical Integration according to step one, step 2, step 3 and step 4;Obtain [Σkr];
Step 6: calculate modal response covarianceForm modal response covariance square
Battle array;
Step 7: calculate dynamic respond covariance matrix [Σx]=[Ψ]T[Σy] [Ψ], and then it is the loudest to try to achieve it
The covariance matrix answered.
Detailed description of the invention two: present embodiment is unlike detailed description of the invention one: based on wind in described step one
The body structure surface pulsating wind pressure time course data p recorded is tested in holeiT (), calculates the blast time course data of each wind load load(ing) point
Meansigma methodsThe standard deviation sigma of blast time course datapi, the auto-correlation function R of blast time course datapi(τ) with blast time course data
Coherent function Cohpij(ω);Detailed process is:
Solved by the mean function in MATLAB;σpiSolved by the std function in MATLAB;Rpi(τ) by MATLAB
Xcorr function solve;Cohpij(ω) solved by the mscohere function in MATLAB.
Other step and parameter are identical with detailed description of the invention one.
Detailed description of the invention three: present embodiment is unlike detailed description of the invention one or two: adopt in described step 2
With exponential function matching auto-correlation function, obtain the crest frequency of the dimensionless blast spectral curve of the wind load load(ing) point of i point
ωmi;Detailed process is:
Use exponential function matching auto-correlation function Rpi(τ)=exp (-ωmi| τ |), obtain the wind load load(ing) point of i point
The crest frequency ω of dimensionless blast spectral curvemi。
Other step and parameter are identical with detailed description of the invention one or two.
Detailed description of the invention four: present embodiment is unlike one of detailed description of the invention one to three: described step 3
Middle employing exponential function matching coherent function, obtains the relevant index k of this structurec, calculate the wind load phase responsibility of i, j 2
Number;Detailed process is:
Use exponential function matching coherent functionObtain the relevant index k of this structurec,
Calculate the wind load coherence factor of i, j 2
In formula, DijRepresenting two wind load load(ing) point i, j spacings, U represents with reference to wind speed, cijWind load for i, j 2
Coherence factor, for nonnegative real number.
Other step and parameter are identical with one of detailed description of the invention one to three.
Detailed description of the invention five: present embodiment is unlike one of detailed description of the invention one to four: described step 5
Middle according to step one, step 2, step 3 and step 4 calculating Analytical Integration;Obtain [Σkr];Detailed process is:
Analytical Integration is calculated according to step one, step 2, step 3 and step 4;
In formula,
σijkrFor k, r rank coupled mode at the coupling wind vibration response root-mean-square of i, j point-to-point transmission;
σpiFor the wind load standard deviation of i point, σpjWind load standard deviation for j point;
ωmiThe crest frequency of dimensionless blast spectral curve for the wind load load(ing) point of i point;ωmjWind load for j point
The crest frequency of the dimensionless blast spectral curve of load(ing) point;
ωnkFor the structure k rank natural frequency of vibration;ωnrFor the structure r rank natural frequency of vibration;
ξnkFor structure k rank damping ratio;ξnrFor for structure r rank damping ratio;
cijWind load coherence factor for i, j 2;
Θ (ω) is the molecule multinomial of integrand;θqFor the polynomial coefficient of molecule;ω2qFor the 2q power of frequency, q
=0,1,2,3;
Denominator complex polynomails for integrand;For complex frequency;λsCoefficient for denominator polynomials;For the s power of complex frequency, s=1~7;;
I1Molecule determinant for integration;I0Denominator determinant for integration;
Form W2Individual Q rank square formation, is expressed as
In formula, i, j=1,2 ... Q;K, r=1,2 ... W;W is to calculate selected rank number of mode, and W is positive integer, according to
Engineering experience determines, generally, and 10≤W≤N;σ1QkrFor σijkrMiddle i=1, j=Q;σ2QkrFor σijkrMiddle i=2, j=Q;With this type of
Push away, σQQkrFor σijkrMiddle a=Q, b=Q;[Σkr] it is the intermediate variable matrix of k, r rank coupled mode;
First calculate the situation of k=r, obtain [Σkr], calculateWherein, k=1,
2,…,W;WillIt is ordered as from big to smallWherein, km=1,2 ..., W, m=1,2 ..., W is the sequence number of sequence;
In formula,Variance for kth rank modal response;[Σkk] be the intermediate variable matrix of kth order mode state, i.e. [Σkr]
The situation of middle k=r;
For front Z item so thatCalculate [the Σ of k < rkr], according to symmetry principle, [Σrk]=
[Σkr], obtain k > [the Σ of rkr];
Described, 1≤Z≤W;1≤n≤W;
For kthmThe variance of rank modal response, i.e.Take k=kmResult.
Other step and parameter are identical with one of detailed description of the invention one to four.
Detailed description of the invention six: present embodiment is unlike one of detailed description of the invention one to five: described molecule is many
The coefficient θ of item formulaqCoefficient lambda with denominator polynomialssParticularly as follows:
Molecule multinomial Θ (ω) about frequencies omega is merged multinomial for frequencies omega, obtains ω2qCoefficient θq, q
=0,1,2,3;
To about complex frequencyDenominator complex polynomailsFor complex frequencyMerge multinomial, obtainCoefficient
λs, s=1~7.
Other step and parameter are identical with one of detailed description of the invention one to five.
Detailed description of the invention seven: present embodiment is unlike one of detailed description of the invention one to six: described I1And I0Tool
Body formula is:
Work as cijWhen ≠ 0,
Work as cijWhen=0,
Other step and parameter are identical with one of detailed description of the invention one to six.
Detailed description of the invention eight: present embodiment is unlike one of detailed description of the invention one to seven: described step 6
Middle calculating modal response covarianceForm modal response covariance matrix;Detailed process
For:
Calculate modal response covarianceForm modal response covariance matrix
Other step and parameter are identical with one of detailed description of the invention one to seven.
Detailed description of the invention nine: present embodiment is unlike one of detailed description of the invention one to eight: described step 7
Middle calculating dynamic respond covariance matrix [Σx]=[Ψ]T[Σy] [Ψ], and then its covariance square arbitrarily responded can be tried to achieve
Battle array;Detailed process is:
Calculate dynamic respond covariance matrix [Σx]=[Ψ]T[Σy] [Ψ], and then its association side arbitrarily responded can be tried to achieve
Difference matrix [Σs]=[I]T[Σx][I]。
Other step and parameter are identical with one of detailed description of the invention one to eight.
Employing following example checking beneficial effects of the present invention:
Embodiment one:
The present embodiment efficient frequency domain estimation method of a kind of wind vibration response based on auto-correlation function is specifically according to following step
Rapid preparation:
Take Practical Project fan-shaped plan rack example of encorbelmenting to illustrate and verify.
Step one~three: blast time-histories wind tunnel test recorded carries out statistical modeling analysis, as Fig. 2 a, Fig. 2 b, Fig. 2 c,
Shown in, obtain average, pulsating wind pressure, the wind load frequency of each point and relevant index, and utilize the standard of blast spectral test model
Really property.
Step 4: FEM (finite element) model is carried out quality, rigidity, the extraction of damping equal matrix, and carries out model analysis, result
As shown in Fig. 3 a, Fig. 3 b, Fig. 3 c, Fig. 3 d.
Step 5: calculate integration, forms modal response standard deviation, as shown in fig. 4 a, it can be seen that the coupling feelings of front 2 × 2
Condition.
Step 6: calculating modal response covariance matrix, result is as shown in Figure 4 b.
Step 7: computation structure dynamic respond standard deviation, as shown in Figure 5.
For verifying the effectiveness of this invention, the most also result of calculation is contrasted with traditional method, found tradition wind
Vibration analysis method computational efficiency is relatively low, and uses the traditional frequency domain computational methods not considering modal coupling, although efficiency is higher, but
Calculate error relatively big (reaching 16.7%).The method that the present invention proposes is while ensureing computational accuracy, and computational efficiency is greatly promoted.
Illustrate to the method achieve the efficient calculating of large and complex structure wind vibration response, there is practical value.
The present invention also can have other various embodiments, in the case of without departing substantially from present invention spirit and essence thereof, and this area
Technical staff is when making various corresponding change and deformation according to the present invention, but these change accordingly and deformation all should belong to
The protection domain of appended claims of the invention.