CN106326669A - Method for determining flutter velocity of slender ship based on transfer function - Google Patents

Abstract

The invention relates to a method for determining flutter velocity of a slender ship based on transfer function, and belongs to the technical field of aeroelesticity of flight vehicles. The method comprises the following steps: 1, obtaining a differential equation of the flutter of the slender ship through a differential equation of transverse bending vibration of beams and a model of a quasi-steady aerodynamic force of the slender ship. 2, dispersing the slender ship into slender ship units, and then proceeding Fourier conversion of the differential equation of the flutter of the slender ship units within the units, and 3, the flutter velocity of the slender ship is solved further by the transfer function. According to the method for determining flutter velocity of a slender ship based on transfer function, instead of adopting the low-order mode of the transverse vibration of the beams to describe approximately the vibration of the slender ship, differential equation of transverse bending vibration of beams is adopted to describe accurately the transverse vibration of the slender, therefore the flutter velocity of the slender ship calculated is more accurate. Compared with the prior technology, the method for solving the flutter velocity of the slender ship is much simpler.

Description

Slender bodies flutter speed based on transmission function determines method

Technical field

The present invention relates to a kind of slender bodies flutter speed based on transmission function and determine method, belong to flight vehicle aerodynamic elastic Technical field.

Background technology

For airframe, missile airframe, the plug-in rocket of wing etc., when slenderness ratio is bigger (generally higher than about 10), It is referred to as slender bodies.Slender bodies tremor refers to the dynamic instability in aircraft flight, and flight speed now is referred to as quivering Vibration velocity degree.The self-excited vibration occurred when flight reaches flutter speed, great majority all can cause catastrophic consequence.Tremor divides The condition that tremor occurs is obtained in analysis exactly, that is obtains flutter speed, and seeks to avoid quivering in the range of vehicle flight speeds Shake occur measure, and then find improve tremor critical velocity method.

On the one hand when carrying out slender bodies flutter analysis and calculate, it is to be appreciated that the dynamics of slender bodies structure, another Aspect is it is to be appreciated that the aerodynamic characteristics of slender bodies body structure surface UNSTEADY FLOW.At present, in order to avoid analysis meter in engineering Calculate excessively complicated, reduce slender bodies vibrational degrees of freedom when slender bodies Structural Dynamics calculates as far as possible, generally ignore high order mode, Lower mode merely with slender bodies vibration carrys out the bending vibration displacement of approximate representation slender bodies.Then, accurate in conjunction with slender bodies Permanent gas dynamic theory, solves slender bodies flutter speed.

It practice, slender bodies is unlimited multivariant elastomer, the most not restriction to degree of freedom in flutter analysis. Therefore, the shortcoming of above-mentioned computational methods is: decrease the degree of freedom of slender bodies vibration in order to avoid excessively complicated calculating, right Slender bodies actual vibration has made approximation, makes the accuracy of result of calculation reduce.

Summary of the invention

The technical problem to be solved is to provide one, and computational methods are simple and direct, result of calculation is accurate based on biography The slender bodies flutter speed of delivery function determines method.

The present invention solves the required technical scheme used of its technical problem:

A kind of slender bodies flutter speed based on transfer function method determines that method, described method utilize lateral deflection of beam to shake The dynamic differential equation and slender bodies quasi-steady aerodynamic force model, draw the slender bodies tremor differential equation, then to described slender bodies from Dissipate for slender bodies unit, in unit, the slender bodies unit tremor differential equation is carried out Fourier conversion, then uses transmission function Method solves slender bodies flutter speed；

One, described method needs the computing formula of foundation:

A. lateral deflection of beam oscillatory differential equation and slender bodies quasi-steady aerodynamic force model is utilized to set up slender bodies tremor The differential equation:

The present invention utilizes lateral deflection of beam oscillatory differential equation formula (1) to describe the bending vibration of slender bodies:

$\frac{{\∂}^2\left(EI\frac{{\∂}^{2}h}{\∂y^2}\right)}{\∂y^2}+m\frac{{\∂}^{2}h}{\∂t^2}+\mathrm{\Δ}p=0---\left(1\right)$

In formula, h is slender bodies bending vibration displacement, and unit is rice, and EI is slender bodies bending rigidity, units Newtons rice², M is slender bodies linear mass, and unit is kilogram, and △ p is the aerodynamic force of slender bodies unit length, and unit is Newton/meter, y For the long axial coordinate of light face type, unit is rice, and t is the time, and unit is the second,For slender bodies bending vibration displacement h to arbor to seat The second-order partial differential coefficient of mark y,For the slender bodies bending vibration displacement h second-order partial differential coefficient to time t；

According to slender bodies quasi-steady aerodynamic force model, the aerodynamic force △ p obtaining slender bodies unit length is following formula (2):

$\mathrm{\Δ}p=-\{\mathrm{\ρ}A(\frac{{\∂}^{2}h}{\∂t^2}+V\frac{{\∂}^{2}h}{\∂t\∂y})+\mathrm{\ρ}V\frac{\∂A}{\∂y}(\frac{\∂h}{\∂t}+V\frac{\∂h}{\∂y})+\mathrm{\ρ}AV(\frac{{\∂}^{2}h}{\∂t\∂y}+V\frac{{\∂}^{2}h}{\∂y^2})\}---\left(2\right)$

In formula, V is air speed, and unit is meter per second, and A is slender bodies cross-sectional area, and unit is square metre, and ρ is atmospheric density, Unit is kilograms per cubic meter；

(2) formula substitution (1) formula is obtained slender bodies tremor differential equation (3):

$\frac{{\∂}^2\left(EI\frac{{\∂}^{2}h}{{\mathrm{ay}}^2}\right)}{\∂y^2}+m\frac{{\∂}^{2}h}{\∂t^2}-\mathrm{\ρ}A(\frac{{\∂}^{2}h}{\∂t^2}+V\frac{{\∂}^{2}h}{\∂t\∂y})-\mathrm{\ρ}V\frac{\∂A}{\∂y}(\frac{\∂h}{\∂t}+V\frac{\∂h}{\∂y})-\mathrm{\ρ}AV(\frac{{\∂}^{2}h}{\∂t\∂y}+V\frac{{\∂}^{2}h}{\∂y^2})=0---\left(3\right)$

B. axially divide n slender bodies unit along slender bodies, set up the tremor differential equation of slender bodies unit:

Usually, physical parameter EI of slender bodies, m, A etc. are axially changes along it, solve for convenience, according to limited Unit's method, axially divides n unit along slender bodies, in order to make the unity of form tremor differential equation of slender bodies unit, and order:

$\mathrm{\ξ}=\frac{y-y_j}{l_j}---\left(4\right)$

In formula, y_jFor the distance of jth slender bodies unit front nodal point distance slender bodies forward terminal, l_jFor jth slender bodies list The length of unit, ξ is axial dimensionless coordinate in slender bodies unit, (4) formula understand ξ ∈ (0,1)；

When division unit is abundant, it is assumed that EI, m, A approximately linear change in each slender bodies unit can be managed, Then just like lower linear expression formula in jth slender bodies unit:

$\left\{\begin{array}{c}EI\left(\mathrm{\ξ}\right)=EI\left(0\right)+\frac{EI\left(1\right)-EI\left(0\right)}{l_j}\mathrm{\ξ}\\ m\left(\mathrm{\ξ}\right)=m\left(0\right)+\frac{m\left(1\right)-m\left(0\right)}{l_j}\mathrm{\ξ}\\ A\left(\mathrm{\ξ}\right)=A\left(0\right)+\frac{A\left(1\right)-A\left(0\right)}{l_j}\mathrm{\ξ}\end{array}\right.---\left(5\right)$

In formula, EI (0), m (0), A (0) are EI, m, A value at jth slender bodies unit front nodal point, EI (1), m (1), A (1) is EI, m, A value at jth slender bodies unit posterior nodal point, and EI (ξ), m (ξ), A (ξ) are that EI, m, A are elongated in jth Value at ξ in body unit；

(4) formula and (5) formula are substituted into (3) formula, the tremor differential equation of available slender bodies unit:

$\begin{array}{c}\frac{EI\left(\mathrm{\ξ}\right)}{l_j^4}\frac{{\∂}^{4}h}{\∂{\mathrm{\ξ}}^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\∂}^{3}h}{\∂{\mathrm{\ξ}}^3}-\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{l_j^2}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\\ -\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{l_j}\frac{{\∂}^{2}h}{\∂t\∂\mathrm{\ξ}}-{\mathrm{\ρV}}^2\frac{A\left(1\right)-A\left(0\right)}{l_j^2}\frac{\∂h}{\∂\mathrm{\ξ}}+\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]\frac{{\∂}^{2}h}{\∂t^2}-\mathrm{\ρ}V\frac{A\left(1\right)-A\left(0\right)}{l_j}\frac{\∂h}{\∂t}=0\end{array}---\left(6\right)$

C. transfer function method is utilized to solve slender bodies flutter speed

1. variable about time t in (6) formula is made Fourier conversion, available following formula (7):

$\begin{array}{c}\frac{EI\left(\mathrm{\ξ}\right)}{l_j^4}\frac{{\∂}^{4}h}{\∂{\mathrm{\ξ}}^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^3}-\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{l_j^2}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\\ -\[i\mathrm{\ω}\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{l_j}+V^2\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{l_j^2}\]\frac{\∂h}{\∂\mathrm{\ξ}}+\[{\left(i\mathrm{\ω}\right)}^2\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]-\left(i\mathrm{\ω}\right)V\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{l_j}\]h=0\end{array}---\left(7\right)$

In formula, imaginary numberω is angular frequency, and unit is radian per second；

2. the tremor differential equation of slender bodies unit is rewritten as state space equation form；

For the ease of application transmission function theory, the vectorial η of one state variable of definition_e(ξ, ω) such as following formula (8):

${\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})={\left[\begin{array}{cccc}h& \frac{\∂h}{\∂y}& \frac{{\∂}^{2}h}{\∂y^2}& \frac{{\∂}^{3}h}{\∂y^3}\end{array}\right]}^T---\left(8\right)$

In formula, T represents vector transposition；

Based on state variable vector η_e(ξ, ω), can be write as the form such as following formula (9) of state equation by (7) formula:

$\frac{\∂{\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})}{\∂\mathrm{\ξ}}=F_e(\mathrm{\ξ},\mathrm{\ω},V){\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})+g_e(\mathrm{\ξ},\mathrm{\ω})---\left(9\right)$

In formula, F_e(ξ,ω,V)、g_eThe expression formula of (ξ, ω) can draw according to (7) formula, specific as follows

$\left\{\begin{array}{c}{F}_e(\mathrm{\ξ},\mathrm{\ω},V)=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ A_4& A_3& A_2& A_1\end{array}\right]\\ g_e(y,\mathrm{\ω})=0\end{array}\right.---\left(10\right)$

In formula,

$\left\{\begin{array}{c}{A}_1=-2\frac{EI\left(1\right)-EI\left(0\right)}{EI\left(\mathrm{\ξ}\right)}\\ A_2=\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{EI\left(\mathrm{\ξ}\right)}{l}_j^2\\ A_3=i\mathrm{\ω}\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{EI\left(\mathrm{\ξ}\right)}{l}_j^2+V^2\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{EI\left(\mathrm{\ξ}\right)}{l}_j^2\\ A_4=-\frac{{\left(i\mathrm{\ω}\right)}^2\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]}{EI\left(\mathrm{\ξ}\right)}{l}_j^4+\left(i\mathrm{\ω}\right)V\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{EI\left(\mathrm{\ξ}\right)}{l}_j^3\end{array}\right.---\left(11\right)$

According to transmission function theory, the canonical solution of equation (9) formula is:

η_e(ξ, ω)=H_e(ξ,ω,V)γ_e(ω)+d_e(ξ,ω) (12)

In formula, γ_e(ω) vector formed for the displacement given by boundary condition or power, can manage according to transfer function method It is as follows that opinion obtains expression formula:

${\mathrm{\γ}}_e\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}h(0,\mathrm{\ω})\\ \frac{\∂h(0,\mathrm{\ω})}{\∂\mathrm{\ξ}}\\ h(1,\mathrm{\ω})\\ \frac{\∂h(1,\mathrm{\ω})}{\∂\mathrm{\ξ}}\end{array}\right\}---\left(13\right)$

In formula, H_e(ξ,ω,V)、f_eThe expression formula of (ξ, ω) obtains also dependent on transfer function method theory, as follows:

$\left\{\begin{array}{c}{H}_e(\mathrm{\&xi;},\mathrm{\&omega;},V)={\mathrm{\&Phi;}}_F(\mathrm{\&xi;},0,\mathrm{\&omega;},V){\&lsqb;M_b{\mathrm{\&Phi;}}_F(0,0,\mathrm{\&omega;},V)+N_b{\mathrm{\&Phi;}}_F(1,0,\mathrm{\&omega;},V)\&rsqb;}^{-1}\\ d_e(\mathrm{\&xi;},\mathrm{\&omega;})={\&Integral;}_0^{1}G(\mathrm{\&xi;},\mathrm{\&zeta;},\mathrm{\&omega;},V)g_e(\mathrm{\&zeta;},\mathrm{\&omega;})d\mathrm{\&zeta;}\\ \mathrm{G(}\mathrm{\&xi;},\mathrm{\&zeta;},\mathrm{\&omega;},V\mathrm{)=}\left\{\begin{array}{cc}{H}_e(\mathrm{\&xi;},\mathrm{\&omega;},V)M_b{\mathrm{\&Phi;}}_F(0,\mathrm{\&zeta;},\mathrm{\&omega;},V)& \mathrm{\&zeta;}<\mathrm{\&xi;}\\ -H_e(\mathrm{\&xi;},\mathrm{\&omega;},V)N_b{\mathrm{\&Phi;}}_F(1,\mathrm{\&zeta;},\mathrm{\&omega;},V)& \mathrm{\&zeta;}>\mathrm{\&xi;}\end{array}\right.\\ {\mathrm{\&Phi;}}_F(\mathrm{\&xi;},\mathrm{\&zeta;},\mathrm{\&omega;},V)=e^{F_e(\mathrm{\&xi;},\mathrm{\&omega;},V)\mathrm{\&zeta;}}\end{array}\right.---\left(14\right)$

In formula, variable ζ ∈ (0,1), M_bFor slender bodies unit front end boundary condition selection matrix, N_bAfter slender bodies unit End boundary condition selection matrix, can obtain according to transfer function method theory, as follows:

$M_b=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],N_b=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]---\left(15\right)$

Owing to (10) formula providing g_e(y, ω)=0, is understood d by (14) formula_e(ξ, ω)=0, thus (12) formula can simplify For:

η_e(ξ, ω)=H_e(ξ,ω,V)γ_e(ω) (16)

3. unit equilibrium equation is set up according to slender bodies internal force relation；

From moment of flexure and the shearing formula of mechanics of materials beam, on slender bodies cross section, the expression formula of moment of flexure and shearing is writeable For:

$\left\{\begin{array}{c}M(y,\mathrm{\&omega;})=EI\frac{{\&part;}^{2}h}{\&part;y^2}\\ F(y,\mathrm{\&omega;})=\frac{\&part;\left(EI\frac{{\&part;}^{2}h}{\&part;y^2}\right)}{\&part;y}\end{array}\right.---\left(17\right)$

In slender bodies unit, and utilize EI linear change in unit it is assumed that (17) formula can be written as:

${\mathrm{\&sigma;}}_e(\mathrm{\&xi;},\mathrm{\&omega;})=\left\{\begin{array}{c}M(\mathrm{\&xi;},\mathrm{\&omega;})\\ F(\mathrm{\&xi;},\mathrm{\&omega;})\end{array}\right\}=\left\{\begin{array}{c}\frac{EI\left(\mathrm{\&xi;}\right)}{l_i^2}\frac{{\&part;}^{2}h}{\&part;{\mathrm{\&xi;}}^2}\\ \frac{EI\left(\mathrm{\&xi;}\right)}{l_i^3}\frac{{\&part;}^{3}h}{\&part;{\mathrm{\&xi;}}^3}+\frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}\frac{{\&part;}^{2}h}{\&part;{\mathrm{\&xi;}}^2}\end{array}\right\}---\left(18\right)$

Introduce state variable vector (8) formula, thus (18) can be written as matrix form:

σ_e(ξ, ω)=E_eη(ξ)η_e(ξ,ω) (19)

In formula, E_eη(ξ) its expression formula can be obtained according to (18) formula as follows:

$E_{e\mathrm{\&eta;}}\left(\mathrm{\&xi;}\right)=\left[\begin{array}{cccc}0& 0& \frac{EI\left(\mathrm{\&xi;}\right)}{l_i^2}& 0\\ 0& 0& \frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}& \frac{EI\left(\mathrm{\&xi;}\right)}{l_i^3}\end{array}\right]---\left(20\right)$

(16) formula is substituted into (19) formula, available

σ_e(ξ, ω)=E_eη(ξ)H_e(ξ,ω,V)γ_e(ω) (21)

In formula, γ_e(ω) can be considered the generalized displacement at slender bodies unit two node, if ξ takes 0 and 1, then can get thin Relation between generalized internal force and the generalized displacement of long body unit two node, thus set up unit equilibrium equation, it may be assumed that

Make ξ=0, ξ=1, then have

${\mathrm{\&sigma;}}_e\left(\mathrm{\&omega;}\right)=\left\{\begin{array}{c}{\mathrm{\&sigma;}}_e(0,\mathrm{\&omega;})\\ {\mathrm{\&sigma;}}_e(1,\mathrm{\&omega;})\end{array}\right\}=\left[\begin{array}{c}{E}_{e\mathrm{\&eta;}}\left(0\right)H_e(0,\mathrm{\&omega;},V)\\ E_{e\mathrm{\&eta;}}\left(1\right)H_e(1,\mathrm{\&omega;},V)\end{array}\right]{\mathrm{\&gamma;}}_e\left(\mathrm{\&omega;}\right)---\left(22\right)$

Above formula is the unit equilibrium equation characterizing slender bodies cell node generalized internal force with generalized displacement relation, thusCan be considered the generalized stiffness matrix of unit, can be designated as:

$K_e(\mathrm{\&omega;},V)=\left[\begin{array}{c}{E}_{e\mathrm{\&eta;}}\left(0\right)H_e(0,\mathrm{\&omega;},V)\\ E_{e\mathrm{\&eta;}}\left(1\right)H_e(1,\mathrm{\&omega;},V)\end{array}\right]---\left(23\right)$

4. slender bodies unit equilibrium equation is assembled and form slender bodies whole machine balancing equation；

Carrying out unit assembling according to FInite Element, available slender bodies whole machine balancing equation is

K (ω, V) γ (ω)=F (ω) (24)

In formula, K (ω, V) can be considered that Bulk stiffness matrix, γ (ω) can be considered that integral node motion vector, F (ω) are each The vector that cell node internal force is assembled into.

5. flutter speed is determined；

For slender bodies, it is not required to during tremor consider gravity, the impact of aerodynamic force suffered by slender bodies is included just by this patent Degree matrix K (ω, V),

Slender bodies is not acted on by other power in addition, thus external applied load is zero, according to cell node internal force and outer load Lotus balances, and can draw

F (ω)=0 (25)

When slender bodies tremor, γ (ω) should have a untrivialo solution, and now whole machine balancing equation must is fulfilled for condition and is

Det [K (ω, V)]=0 (26)

Owing to K (ω, V) is complex matrix, the null essential condition of its determinant be matrix determinant value real part with Imaginary part is zero, i.e.

$\left\{\begin{array}{c}\mathrm{Re}\{\det \&lsqb;K(\mathrm{\&omega;},V)\&rsqb;\}=0\\ \mathrm{Im}\{\det \&lsqb;K(\mathrm{\&omega;},V)\&rsqb;\}=0\end{array}\right.---\left(27\right)$

Matrix K (ω, V) has air speed V and two variablees of circular frequency ω, and (27) formula has two equations just, Ke Yiding Solve.When solving above-mentioned equation group, the multiple solution of equation group may be met, i.e. there are many groups (V, ω) and equation group can be met (27) formula.According to during slender bodies tremor when a certain air speed by being stably changed into instability, thus, one group of solution that air speed V is minimum (V, ω) should be the flutter speed of slender bodies and corresponding flutter frequency.Air speed V and circular frequency ω of equation group will be met, It is designated as V respectively_czAnd ω_cz.Wherein, V_czIt is the flutter speed of slender bodies, ω_czIt is the tremor circular frequency of slender bodies.

Two, the specifically comprising the following steps that of described method

Step (one): slender bodies is divided into n unit, the following physical parameter of measurement slender bodies each unit node:

Slender bodies each unit axial length [l₁ l₂ … l_j … l_n], unit is rice；

Cross-sectional area [A at slender bodies cell node₁ A₂ … A_j … A_n+1], unit is rice²

Linear mass [m at slender bodies unit-node₁ m₂ … m_j … m_n+1], unit is kg/m；

Bending rigidity [EI at slender bodies unit-node₁ EI₂ … EI_j … EI_n+1], unit is Newton meter²；

Atmospheric density ρ, unit is kg/m³；

Step (two): determine air speed V and the approximate range of circular frequency ω that slender bodies flies, it is assumed that for

Step (three):In the range of divide suitable step-length △ V and △ ω, and carry out discrete, air speed V and circle The value of frequencies omega is:

Step (four): take air speed V=V₀, circular frequency ω takes ω successively₀+ j △ ω (j=0,1,2,3 ...), by air speed V and circle The value of frequencies omega substitutes into (11) formula with the slender bodies physical parameter in step (), calculates each unit coefficient A₁、A₂、A₃、A₄'s Value；

Step (five): by coefficient A₁、A₂、A₃、A₄Value substitute into (10) formula, obtain slender bodies each unit matrix F_e(ξ,ω,V) Value；

Step (six): by matrix F_eThe value of (ξ, ω, V) substitutes into (14) formula, in conjunction with (15) formula, calculates H_eThe value of (ξ, ω, V)；

Step (seven): by H_eThe value of (ξ, ω, V) substitutes into (23) formula, in conjunction with (20) formula, computing unit stiffness matrix K_e(ω, V) value；

Step (eight): utilize K_e(ω, V) assembles Bulk stiffness matrix K (ω, V)；

Step (nine): calculate Re{det [K (ω, V)] and Im{det [K (ω, V)] value；

Step (ten): take air speed V=V the most successively₀+ j △ V (j=1,2,3 ...), repeat step (four) to step (ten), meter Calculate Re{det [K (ω, V)] and Im{det [K (ω, V)] value；

Step (11): determine the value of air speed V of satisfied (27) formula, i.e. can determine that flutter speed V of slender bodies_cz。

Beneficial effects of the present invention is as follows:

(1) lateral deflection of beam oscillatory differential equation is utilized to carry out the oscillation crosswise of accurate description slender bodies due to the present invention, And do not use the lower mode of beam oscillation crosswise to carry out the vibration of approximate description slender bodies, so the slender bodies flutter speed calculated is more Add accurately.

(2) to solve the method for slender bodies flutter speed the simplest and the most direct for the present invention.

Accompanying drawing explanation

Fig. 1 is slender bodies front view；

Fig. 2 is slender bodies sectional view；

Fig. 3 is slender bodies structural representation；

Fig. 4 is slender bodies element subdivision schematic diagram；

Fig. 5 be the Re [detA] and Im [detA] of embodiment 1 isogram (V ∈ (0,50), ω ∈ (0,200 π)),

In the accompanying drawings: 1 slender bodies, 2 i-th slender bodiess, 3 front nodal points, 4 posterior nodal points.

Detailed description of the invention

As shown in accompanying drawing 1-5, embodiments of the invention 1 are as follows:

A length of 1.8 meters of certain slender bodies, weight 22 kilograms, a diameter of 0.025 square metre of maximum cross-section.

The concrete calculation procedure of the present embodiment 1 is as follows:

Step (one): slender bodies is divided into 10 unit, the following physical parameter of measurement slender bodies each unit node:

Each unit axial length: l₁=l₂=...=l₁₀=0.18, unit is rice；

Cross-sectional area at each unit node: A₁=0, A₂=A₃=...=A₁₀=0.025, A₁₁=0.02, unit is rice²

Linear mass m at each unit node₁=10, m₂=m₃...=m₁₀=20, unit is kg/m；

Bending rigidity EI at each unit node₁=0.5 × 10³,EI₂=...=EI₁₁=1.1 × 10³, unit is cattle Pause rice²；

Atmospheric density ρ=1.225, unit is kg/m³；

Step (two): determine that air speed V of aircraft wing and the approximate range of circular frequency ω areDivide step-lengthAccording to step (three) in Summary to step (11), calculate Re [detA] and the value of Im [detA].Figure 3 giveIn the range of Re [detA] and the isogram of Im [detA].It is found that giving from Fig. 3 ScopeInterior existence meets the some A of (26) formula, and its coordinate is (14.7,3020), and it represents slender bodies flutter speed V_cz=3020m/s, flutter frequency is 14.7Hz (tremor circular frequency ω_cz=92.32rad/s).

Claims (1)

1. a slender bodies flutter speed based on transmission function determines method, it is characterised in that described method utilizes the horizontal of beam The bending vibration differential equation and slender bodies quasi-steady aerodynamic force model, draw the slender bodies tremor differential equation, then to described carefully The long body tremor differential equation carries out Fourier conversion, then uses transfer function method to solve slender bodies flutter speed；

One, described method needs the computing formula of foundation:

A. lateral deflection of beam oscillatory differential equation and slender bodies quasi-steady aerodynamic force model is utilized to set up slender bodies tremor differential Equation:

The present invention utilizes lateral deflection of beam oscillatory differential equation formula (1) to describe the bending vibration of slender bodies:

$\frac{{\&part;}^2\left(EI\frac{{\&part;}^{2}h}{\&part;y^2}\right)}{\&part;y^2}+m\frac{{\&part;}^{2}h}{\&part;t^2}+\mathrm{\&Delta;}p=0---\left(1\right)$

In formula, h is slender bodies bending vibration displacement, and unit is rice, and EI is slender bodies bending rigidity, units Newtons rice², m is thin Long body linear mass, unit is kilogram, and △ p is the aerodynamic force of slender bodies unit length, and unit is Newton/meter, and y is light face type Long axial coordinate, unit is rice, and t is the time, and unit is the second,For slender bodies bending vibration displacement h to machine axial coordinate y Second-order partial differential coefficient,For the slender bodies bending vibration displacement h second-order partial differential coefficient to time t；

According to slender bodies quasi-steady aerodynamic force model, the aerodynamic force △ p obtaining slender bodies unit length is following formula (2):

$\mathrm{\&Delta;}p=-\{\mathrm{\&rho;}A(\frac{{\&part;}^{2}h}{\&part;t^2}+V\frac{{\&part;}^{2}h}{\&part;t\&part;y})+\mathrm{\&rho;}V\frac{\&part;A}{\&part;y}(\frac{\&part;h}{\&part;t}+V\frac{\&part;h}{\&part;y})+\mathrm{\&rho;}AV(\frac{{\&part;}^{2}h}{\&part;t\&part;y}+V\frac{{\&part;}^{2}h}{\&part;y^2})\}---\left(2\right)$

In formula, V is air speed, and unit is meter per second, and A is slender bodies cross-sectional area, and unit is square metre, and ρ is atmospheric density, unit For kilograms per cubic meter；

(2) formula substitution (1) formula is obtained slender bodies tremor differential equation (3):

$\frac{{\&part;}^2\left(EI\frac{{\&part;}^{2}h}{{\mathrm{ay}}^2}\right)}{{\mathrm{ay}}^2}+m\frac{{\&part;}^{2}h}{\&part;t^2}-\mathrm{\&rho;}A(\frac{{\&part;}^{2}h}{\&part;t^2}+V\frac{{\&part;}^{2}h}{\&part;t\&part;y})-\mathrm{\&rho;}V\frac{\&part;A}{ay}(\frac{\&part;h}{\&part;t}+V\frac{\&part;h}{ay})-\mathrm{\&rho;}AV(\frac{{\&part;}^{2}h}{\&part;t\&part;y}+V\frac{{\&part;}^{2}h}{{\mathrm{ay}}^2})=0---\left(3\right)$

B. axially divide n slender bodies unit along slender bodies, set up the tremor differential equation of slender bodies unit:

Usually, physical parameter EI of slender bodies, m, A etc. are axially changes along it, solve for convenience, according to finite element Method, axially divides n unit along slender bodies, in order to make the unity of form tremor differential equation of slender bodies unit, and order:

$\mathrm{\&xi;}=\frac{y-y_j}{l_j}---\left(4\right)$

In formula, y_jFor the distance of jth slender bodies unit front nodal point distance slender bodies forward terminal, l_jFor jth slender bodies unit Length, ξ is axial dimensionless coordinate in slender bodies unit, (4) formula understand ξ ∈ (0,1)；

When division unit is abundant, it is assumed that EI, m, A approximately linear change in each slender bodies unit can be managed, then exist Just like lower linear expression formula in jth slender bodies unit:

$\left\{\begin{array}{c}EI\left(\mathrm{\&xi;}\right)=EI\left(0\right)+\frac{EI\left(1\right)-EI\left(0\right)}{l_j}\mathrm{\&xi;}\\ m\left(\mathrm{\&xi;}\right)=m\left(0\right)+\frac{m\left(1\right)-m\left(0\right)}{l_j}\mathrm{\&xi;}\\ A\left(\mathrm{\&xi;}\right)=A\left(0\right)+\frac{A\left(1\right)-A\left(0\right)}{l_j}\mathrm{\&xi;}\end{array}\right.---\left(5\right)$

In formula, EI (0), m (0), A (0) are EI, m, A value at jth slender bodies unit front nodal point, EI (1), m (1), A (1) For EI, m, A value at jth slender bodies unit posterior nodal point, EI (ξ), m (ξ), A (ξ) are that EI, m, A are at jth slender bodies list Value at the interior ξ of unit；

(4) formula and (5) formula are substituted into (3) formula, the tremor differential equation of available slender bodies unit:

$\begin{array}{c}\frac{EI\left(\mathrm{\&xi;}\right)}{l_j^4}\frac{{\&part;}^{4}h}{\&part;{\mathrm{\&xi;}}^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\&part;}^{3}h}{\&part;{\mathrm{\&xi;}}^3}-\frac{\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)V^2}{l_j^2}\frac{{\&part;}^{2}h}{\&part;{\mathrm{\&xi;}}^2}\\ -\frac{2\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)V}{l_j}\frac{{\&part;}^{2}h}{\&part;t\&part;\mathrm{\&xi;}}-{\mathrm{\&rho;V}}^2\frac{A\left(1\right)-A\left(0\right)}{l_j^2}\frac{\&part;h}{\&part;\mathrm{\&xi;}}+\&lsqb;m\left(\mathrm{\&xi;}\right)-\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)\&rsqb;\frac{{\&part;}^{2}h}{\&part;{\mathrm{\&xi;}}^2}-\mathrm{\&rho;}V\frac{A\left(1\right)-A\left(0\right)}{l_j}\frac{\&part;h}{\&part;t}=0\end{array}---\left(6\right)$

C. transfer function method is utilized to solve slender bodies flutter speed

1. variable about time t in (6) formula is made Fourier conversion, available following formula (7):

$\begin{array}{c}\frac{EI\left(\mathrm{\&xi;}\right)}{l_j^4}\frac{{\&part;}^{4}h}{\&part;{\mathrm{\&xi;}}^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\&part;}^{3}h}{\&part;{\mathrm{\&xi;}}^3}-\frac{\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)V^2}{l_j^2}\frac{{\&part;}^{2}h}{\&part;{\mathrm{\&xi;}}^2}\\ -\&lsqb;i\mathrm{\&omega;}\frac{2\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)V}{l_j}+V^2\frac{\mathrm{\&sigma;}A\left(1\right)-\mathrm{\&sigma;}A\left(0\right)}{l_j^2}\&rsqb;\frac{\&part;h}{\&part;\mathrm{\&xi;}}+\&lsqb;{\left(i\mathrm{\&omega;}\right)}^2\&lsqb;m\left(\mathrm{\&xi;}\right)-\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)\&rsqb;-\left(i\mathrm{\&omega;}\right)V\frac{\mathrm{\&rho;}A\left(1\right)-\mathrm{\&rho;}A\left(0\right)}{l_j}\&rsqb;h=0\end{array}---\left(7\right)$

In formula, imaginary numberω is angular frequency, and unit is radian per second；

2. the tremor differential equation of slender bodies unit is rewritten as state space equation form；

For the ease of application transmission function theory, the vectorial η of one state variable of definition_e(ξ, ω) such as following formula (8):

${\mathrm{\&eta;}}_e(\mathrm{\&xi;},\mathrm{\&omega;})={\left[\begin{array}{cccc}h& \frac{\&part;h}{\&part;y}& \frac{{\&part;}^{2}h}{\&part;y^2}& \frac{{\&part;}^{3}h}{\&part;y^3}\end{array}\right]}^T---\left(8\right)$

In formula, T represents vector transposition；

Based on state variable vector η_e(ξ, ω), can be write as the form such as following formula (9) of state equation by (7) formula:

$\frac{\&part;{\mathrm{\&eta;}}_e(\mathrm{\&xi;},\mathrm{\&omega;})}{\&part;\mathrm{\&xi;}}=F_e(\mathrm{\&xi;},\mathrm{\&omega;},V){\mathrm{\&eta;}}_e(\mathrm{\&xi;},\mathrm{\&omega;})+g_e(\mathrm{\&xi;},\mathrm{\&omega;})---\left(9\right)$

In formula, F_e(ξ,ω,V)、g_eThe expression formula of (ξ, ω) can draw according to (7) formula, specific as follows

$\left\{\begin{array}{c}{F}_e(\mathrm{\&xi;},\mathrm{\&omega;},V)=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ A_4& A_3& A_2& A_1\end{array}\right]\\ g_e(y,\mathrm{\&omega;})=0\end{array}\right.---\left(10\right)$

In formula,

$\left\{\begin{array}{c}{A}_1=-2\frac{EI\left(1\right)-EI\left(0\right)}{EI\left(\mathrm{\&xi;}\right)}\\ A_2=\frac{\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)V^2}{EI\left(\mathrm{\&xi;}\right)}{l}_j^2\\ A_3=i\mathrm{\&omega;}\frac{2\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)V}{EI\left(\mathrm{\&xi;}\right)}{l}_j^3+V^2\frac{\mathrm{\&rho;}A\left(1\right)-\mathrm{\&rho;}A\left(0\right)}{EI\left(\mathrm{\&xi;}\right)}{l}_j^2\\ A_4=-\frac{{\left(i\mathrm{\&omega;}\right)}^2\&lsqb;m\left(\mathrm{\&xi;}\right)-\mathrm{\&rho;}A\left(\mathrm{\&xi;}\right)\&rsqb;}{EI\left(\mathrm{\&xi;}\right)}{l}_j^4+\left(i\mathrm{\&omega;}\right)V\frac{\mathrm{\&rho;}A\left(1\right)-\mathrm{\&rho;}A\left(0\right)}{EI\left(\mathrm{\&xi;}\right)}{l}_j^3\end{array}\right.---\left(11\right)$

According to transmission function theory, the canonical solution of equation (9) formula is:

η_e(ξ, ω)=H_e(ξ,ω,V)γ_e(ω)+d_e(ξ,ω) (12)

In formula, γ_e(ω) vector formed for the displacement given by boundary condition or power, can obtain according to transfer function method is theoretical As follows to expression formula:

${\mathrm{\&gamma;}}_e\left(\mathrm{\&omega;}\right)=\left\{\begin{array}{c}h(0,\mathrm{\&omega;})\\ \frac{\&part;h(0,\mathrm{\&omega;})}{\&part;\mathrm{\&xi;}}\\ h(1,\mathrm{\&omega;})\\ \frac{\&part;h(1,\mathrm{\&omega;})}{\&part;\mathrm{\&xi;}}\end{array}\right\}---\left(13\right)$

In formula, H_e(ξ,ω,V)、f_eThe expression formula of (ξ, ω) obtains also dependent on transfer function method theory, as follows:

$\left\{\begin{array}{c}{H}_e(\mathrm{\&xi;},\mathrm{\&omega;},V)={\mathrm{\&Phi;}}_F(\mathrm{\&xi;},0,\mathrm{\&omega;},V){\&lsqb;M_b{\mathrm{\&Phi;}}_F(0,0,\mathrm{\&omega;},V)+N_b{\mathrm{\&Phi;}}_F(1,0,\mathrm{\&omega;},V)\&rsqb;}^{-1}\\ d_e(\mathrm{\&xi;},\mathrm{\&omega;})={\&Integral;}_0^1(\mathrm{\&xi;},\mathrm{\&zeta;},\mathrm{\&omega;},V)g_e(\mathrm{\&zeta;},\mathrm{\&omega;})d\mathrm{\&zeta;}\\ G(\mathrm{\&xi;},\mathrm{\&zeta;},\mathrm{\&omega;},V)=\left\{\begin{array}{cc}{H}_e(\mathrm{\&xi;},\mathrm{\&omega;},V)M_b{\mathrm{\&Phi;}}_F(0,\mathrm{\&zeta;},\mathrm{\&omega;},V)& \mathrm{\&zeta;}<\mathrm{\&xi;}\\ -H_e(\mathrm{\&xi;},\mathrm{\&omega;},V)N_b{\mathrm{\&Phi;}}_F(1,\mathrm{\&zeta;},\mathrm{\&omega;},V)& \mathrm{\&zeta;}>\mathrm{\&xi;}\end{array}\right.\\ {\mathrm{\&Phi;}}_F(\mathrm{\&xi;},\mathrm{\&zeta;},\mathrm{\&omega;},V)=e^{F_e(\mathrm{\&xi;},\mathrm{\&omega;},V)\mathrm{\&zeta;}}\end{array}\right.---\left(14\right)$

In formula, variable ζ ∈ (0,1), M_bFor slender bodies unit front end boundary condition selection matrix, N_bFor slender bodies unit back end edge Boundary's condition selection matrix, can obtain according to transfer function method theory, as follows:

$M_b=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],N_b=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]---\left(15\right)$

Owing to (10) formula providing g_e(y, ω)=0, is understood d by (14) formula_e(ξ, ω)=0, thus (12) formula can be reduced to:

η_e(ξ, ω)=H_e(ξ,ω,V)γ_e(ω) (16)

3. unit equilibrium equation is set up according to slender bodies internal force relation；

From moment of flexure and the shearing formula of mechanics of materials beam, on slender bodies cross section, the expression formula of moment of flexure and shearing can be written as:

$\left\{\begin{array}{c}M(y,\mathrm{\&omega;})=EI\frac{{\&part;}^{2}h}{\&part;y^2}\\ F(y,\mathrm{\&omega;})=\frac{\&part;\left(EI\frac{{\&part;}^{2}h}{\&part;y^2}\right)}{\&part;y}\end{array}\right.---\left(17\right)$

In slender bodies unit, and utilize EI linear change in unit it is assumed that (17) formula can be written as:

${\mathrm{\&sigma;}}_e(\mathrm{\&xi;},\mathrm{\&omega;})=\left\{\begin{array}{c}M(\mathrm{\&xi;},\mathrm{\&omega;})\\ F(\mathrm{\&xi;},\mathrm{\&omega;})\end{array}\right\}=\left\{\begin{array}{c}\frac{EI\left(\mathrm{\&xi;}\right)}{l_i^2}\frac{{\&part;}^{2}h}{\&part;{\mathrm{\&xi;}}^2}\\ \frac{EI\left(\mathrm{\&xi;}\right)}{l_i^3}\frac{{\&part;}^{3}h}{\&part;{\mathrm{\&xi;}}^3}+\frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}\frac{{\&part;}^{2}h}{\&part;{\mathrm{\&xi;}}^2}\end{array}\right\}---\left(18\right)$

Introduce state variable vector (8) formula, thus (18) can be written as matrix form:

σ_e(ξ, ω)=E_eη(ξ)η_e(ξ,ω) (19)

In formula, E_eη(ξ) its expression formula can be obtained according to (18) formula as follows:

$E_{e\mathrm{\&eta;}}\left(\mathrm{\&xi;}\right)=\left[\begin{array}{cccc}0& 0& \frac{EI\left(\mathrm{\&xi;}\right)}{l_i^2}& 0\\ 0& 0& \frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}& \frac{EI\left(\mathrm{\&xi;}\right)}{l_i^3}\end{array}\right]---\left(20\right)$

(16) formula is substituted into (19) formula, available

σ_e(ξ, ω)=E_eη(ξ)H_e(ξ,ω,V)γ_e(ω) (21)

In formula, γ_e(ω) can be considered the generalized displacement at slender bodies unit two node, if ξ takes 0 and 1, then can get slender bodies Relation between generalized internal force and the generalized displacement of unit two node, thus set up unit equilibrium equation, it may be assumed that

Make ξ=0, ξ=1, then have

${\mathrm{\&sigma;}}_e\left(\mathrm{\&omega;}\right)=\left\{\begin{array}{c}{\mathrm{\&sigma;}}_e(0,\mathrm{\&omega;})\\ {\mathrm{\&sigma;}}_e(1,\mathrm{\&omega;})\end{array}\right\}=\left[\begin{array}{c}{E}_{e\mathrm{\&eta;}}\left(0\right)H_e(0,\mathrm{\&omega;},V)\\ E_{e\mathrm{\&eta;}}\left(1\right)H_e(1,\mathrm{\&omega;},V)\end{array}\right]{\mathrm{\&gamma;}}_e\left(\mathrm{\&omega;}\right)---\left(22\right)$

Above formula is the unit equilibrium equation characterizing slender bodies cell node generalized internal force with generalized displacement relation, thusCan be considered the generalized stiffness matrix of unit, can be designated as:

$K_e(\mathrm{\&omega;},V)=\left[\begin{array}{c}{E}_{e\mathrm{\&eta;}}\left(0\right)H_e(0,\mathrm{\&omega;},V)\\ E_{e\mathrm{\&eta;}}\left(1\right)H_e(1,\mathrm{\&omega;},V)\end{array}\right]---\left(23\right)$

4. slender bodies unit equilibrium equation is assembled and form slender bodies whole machine balancing equation；

Carrying out unit assembling according to FInite Element, available slender bodies whole machine balancing equation is

K (ω, V) γ (ω)=F (ω) (24)

In formula, K (ω, V) can be considered that Bulk stiffness matrix, γ (ω) can be considered that integral node motion vector, F (ω) are each unit The vector that node reaction forces is assembled into.

5. flutter speed is determined；

For slender bodies, it is not required to during tremor consider gravity, the impact of aerodynamic force suffered by slender bodies is included at rigidity square by this patent Battle array K (ω, V), slender bodies is not acted on by other power in addition, thus external applied load is zero, according to cell node internal force with outer Counterweight balance, can draw

F (ω)=0 (25)

When slender bodies tremor, γ (ω) should have a untrivialo solution, and now whole machine balancing equation must is fulfilled for condition and is

Det [K (ω, V)]=0 (26)

Owing to K (ω, V) is complex matrix, the null essential condition of its determinant is real part and the imaginary part of matrix determinant value It is zero, i.e.

$\left\{\begin{array}{c}\mathrm{Re}\{\det \&lsqb;K(\mathrm{\&omega;},V)\&rsqb;\}=0\\ \mathrm{Im}\{\det \&lsqb;K(\mathrm{\&omega;},V)\&rsqb;\}=0\end{array}\right.---\left(27\right)$

Matrix K (ω, V) has air speed V and two variablees of circular frequency ω, and (27) formula has two equations just, can surely solve.Ask When solving above-mentioned equation group, the multiple solution of equation group may be met, i.e. there are many groups (V, ω) and equation group (27) formula can be met. According to during slender bodies tremor when a certain air speed by being stably changed into instability, thus, the minimum one group of solution (V, ω) of air speed V should Flutter speed and corresponding flutter frequency for slender bodies.Air speed V and circular frequency ω of equation group will be met, be designated as respectively V_czAnd ω_cz.Wherein, V_czIt is the flutter speed of slender bodies, ω_czIt is the tremor circular frequency of slender bodies.

Two, the specifically comprising the following steps that of described method

Step (one): slender bodies is divided into n unit, the following physical parameter of measurement slender bodies each unit node:

Slender bodies each unit axial length [l₁ l₂ … l_j … l_n], unit is rice；

Cross-sectional area [A at slender bodies cell node₁ A₂ … A_j … A_n+1], unit is rice²

Linear mass [m at slender bodies unit-node₁ m₂ … m_j … m_n+1], unit is kg/m；

Bending rigidity [EI at slender bodies unit-node₁ EI₂ … EI_j … EI_n+1], unit is Newton meter²；

Atmospheric density ρ, unit is kg/m³；

Step (two): determine air speed V and the approximate range of circular frequency ω that slender bodies flies, it is assumed that for

Step (three):In the range of divide suitable step-length △ V and △ ω, and carry out discrete, air speed V and circular frequency The value of ω is:

Step (four): take air speed V=V₀, circular frequency ω takes ω successively₀+ j △ ω (j=0,1,2,3 ...), by air speed V and circular frequency The value of ω substitutes into (11) formula with the slender bodies physical parameter in step (), calculates each unit coefficient A₁、A₂、A₃、A₄Value；

Step (five): by coefficient A₁、A₂、A₃、A₄Value substitute into (10) formula, obtain slender bodies each unit matrix F_e(ξ, ω, V's) Value；

Step (six): by matrix F_eThe value of (ξ, ω, V) substitutes into (14) formula, in conjunction with (15) formula, calculates H_eThe value of (ξ, ω, V)；

Step (seven): by H_eThe value of (ξ, ω, V) substitutes into (23) formula, in conjunction with (20) formula, computing unit stiffness matrix K_e(ω, V's) Value；

Step (eight): utilize K_e(ω, V) assembles Bulk stiffness matrix K (ω, V)；

Step (nine): calculate Re{det [K (ω, V)] and Im{det [K (ω, V)] value；

Step (ten): take air speed V=V the most successively₀+ j △ V (j=1,2,3 ...), repeat step (four) to step (ten), calculate Re { det [K (ω, V)] } and Im{det [K (ω, V)] } value；

Step (11): determine the value of air speed V of satisfied (27) formula, i.e. can determine that flutter speed V of slender bodies_cz。