CN105930669A - Method for calculating rigidity breathing function of non-gravity-dominated cracked rotor - Google Patents

Abstract

The invention provides a method for calculating a rigidity breathing function of a non-gravity-dominated cracked rotor. The method is high in solving precision and high in suitability. The method comprises the following steps of firstly, at a rotary vortex difference angle, determining stress tension and compression distribution of a rotor cross section at a crack in a state of assuming an initial restoring force direction to obtain a crack section closure region, wherein the rotor cross section at the crack is a crack section, and the crack section closure region is a closure section; secondly, performing calculation by utilizing a bending theory to obtain the rotor rigidity in the state; thirdly, correcting restoring force in the state by utilizing the calculated rotor rigidity, and performing iteration in sequence until the rigidity is converged; traversing all discrete values in a definition domain of the rotary vortex difference angle to obtain discrete rigidity values corresponding to the discrete values; and finally, fitting the rigidity in the discrete state into a triangular basis function continuous expression on the rotary vortex difference angle, and obtaining the rigidity of the non-gravity-dominated cracked rotor, changed according to a breathing state.

Description

Calculation method for non-gravity dominant crack rotor stiffness respiration function

Technical Field

The invention relates to the field of mechanical dynamics, in particular to a method for calculating a non-gravity dominant crack rotor stiffness respiration function.

Background

Large rotary machines such as steam turbines and compressors are the most critical and widely used mechanical equipment in the national electrical and chemical industries. Because the structure of the large-scale equipment is complex, the working condition change is large, and the fatigue crack of the equipment material is easily caused. Rotor cracking is one of the most major failures during the actual operation of the equipment. In order to monitor rotor cracks and avoid catastrophic accidents causing machine damage and human death, it is necessary to deeply research the dynamic characteristics of a crack rotor system. Accurate calculation of the respiratory stiffness of the crack rotor is a precondition for researching the dynamic characteristics of the crack rotor.

Currently, there are two types of transverse crack stiffness models: a full-open crack model and an open-close crack model. The opening and closing crack model refers to that cracks are opened sometimes, closed sometimes and half opened and half closed sometimes in the running process of the rotor, and the running condition of the actual rotor can be reflected due to comprehensive consideration of three states of the crack sections. The existing open-close crack stiffness analytic solution based on the neutral axis theory can reflect the breathing effect of the crack, but neglect the neutral axis shift caused by the asymmetry of the crack section, and assume that the neutral axis passes through the center of the crack section, so that the model has the condition of inaccurate determination of the crack opening area.

The improved rigidity calculation of the opening and closing crack based on the neutral axis theory is carried out under the assumption that the stress concentration of the crack end is not considered, when the crack is in a fully-opened state or a half-opened and half-closed state, the influence of a crack opening area on the position of the neutral axis of the section of the crack is considered, the neutral axis deviating from the centroid of the intact section by a certain distance is obtained, a part of crack surface above the neutral axis is in a compressive stress area and is in a closed state, and a part of crack surface below the neutral axis is in a tensile stress area and is in an open state. However, the improvement considers that the direction of the elastic restoring force is consistent with the axial displacement direction of the rotor, then, due to the existence of cracks, the anisotropic rigidity of the rotor is unequal, the direction of the restoring force is not consistent with the axial displacement direction, and the method still cannot accurately calculate the rigidity of the cracked rotor.

Disclosure of Invention

Aiming at the problems in the prior art, the invention provides a method for calculating the stiffness and the breathing function of a non-gravity dominated cracked rotor, which has the advantages of high solving precision and strong applicability, and can accurately solve the stiffness of rotors with different section shapes and breathing cracks in any rotating state.

The invention is realized by the following technical scheme:

a non-gravity dominated crack rotor stiffness breathing function calculation method comprises the steps that firstly, at a certain vortex difference angle, stress tension and compression distribution of a rotor cross section at a crack under the state is determined by assuming an initial restoring force direction, and a crack cross section closed area is obtained; wherein, the cross section of the rotor at the crack is the crack section; the crack section closed area is a closed section; then, calculating to obtain the rigidity of the rotor in the state by using a bending theory; correcting restoring force in the state by calculating the rigidity of the rotor, and sequentially iterating until the rigidity is converged; traversing all discrete values in the vortex difference angle definition domain to obtain a corresponding discrete stiffness value; and finally, fitting the rigidity in the discrete state into a triangular basis function continuous expression about the rotor eddy difference angle by using a data fitting method, and obtaining the rigidity of the non-gravity dominant crack rotor changing in a breathing state.

Preferably, the method specifically comprises the following steps:

1) dispersing a domain [0,2 pi ] of the vortex difference angle psi into N equal parts, and dividing the crack section into a plurality of discrete surface units;

2) under a certain selected rotor vortex difference angle, assuming that the restoring force direction and the rotor axis displacement direction are opposite, calculating the stress of each surface unit, and calculating to obtain a closed section according to a neutral layer theory;

3) calculating the position of an inertia main shaft of the current closed section, repeating the step 2, and recalculating the stress of the surface unit of the current closed section until a convergent closed section is obtained;

4) calculating the rigidity of the rotor in the state by using a bending theory through the convergent closed section;

5) calculating the rotor restoring force by the following formula by calculating the rotor stiffness under the rotor vortex difference angle:

$\left[\begin{array}{c}{F}_x\\ F_y\end{array}\right]=-\left[\begin{array}{cc}{k}_x& k_{xy}\\ k_{xy}& k_y\end{array}\right]\left[\begin{array}{c}{x}_{o^{\′}}\\ y_{o^{\′}}\end{array}\right];$

wherein, F_xAnd F_yComponent of restoring force along x and x axes in a fixed coordinate system o-xy, k_x，k_yAnd k_xyThree-dimensional flexural stiffness, x, of the rotor relative to a fixed coordinate system_o'And y_o'Is the axial displacement of the rotor relative to a fixed coordinate system;

6) repeating the steps 2) -5) until the restoring force in the state is converged, and obtaining the rotor rigidity under the rotor vortex difference angle;

7) selecting any non-calculated rotor vortex difference angle, and repeating the steps 2) to 6) until traversing N rotor vortex difference angles in the defined domain; obtaining the corresponding rotor stiffness under each rotational vortex difference angle;

8) and fitting the N rotor rigidities in the discrete rotor eddy difference angle state into a continuous expression of a triangular basis function, namely a rigidity breathing function, by using a data fitting method, so as to obtain the rigidity value of the breathing change.

Further, the fitting of N discrete-state rotor stiffnesses into a continuous expression of a triangular basis function with respect to a rotor vortex difference angle by using a data fitting method includes the following steps:

a step of making the axisymmetrically distributed horizontal and vertical stiffness in the form of a cos (ψ) basis function, wherein,

k_x＝a(1)cos(ψ)+a(2)cos(ψ)²+a(3)cos(ψ)³+a(4)cos(ψ)⁶+a(5)cos(ψ)⁷+a(6)cos(ψ)⁸+a(7)cos(ψ)⁹+a(8)cos(ψ)¹⁰+a(9)cos(ψ)¹⁵+a(10)cos(ψ)¹⁹+a(11)

$k_y=\underset{i=1}{\overset{7}{\Σ}}a\left(i\right)cos{\left(\mathrm{\ψ}\right)}^i+a\left(8\right);$

a step of making the cross stiffness of the centrosymmetric distribution into a Fourier series form, wherein,

$k_{xy}=\underset{i=1}{\overset{5}{\Σ}}a\left(i\right)sin\left(i\mathrm{\ψ}\right)+\underset{i=6}{\overset{10}{\Σ}}a\left(i\right)cos\left(\left(i-5\right)\mathrm{\ψ}\right)+a\left(11\right);$

in the formula 3, a (i) is a undetermined coefficient, and a vortex shedding difference angle psi belongs to [0,2 pi ];

and obtaining coefficients a (i) of the stiffness values under the N discrete vortex difference angles according to a least square data fitting principle by using the obtained stiffness values under the N discrete vortex difference angles, and determining a continuous expression of the stiffness respiratory function.

Further, the centroid and centroid principal axis of inertia position of the current closed section are calculated by the following formula,

in the formula, A_iThe area of the inner face unit of the closed section is shown, and n is the number of the inner face units of the closed section;androtating the centroid coordinate under the coordinate system o ' -x ' y ',andclosed cross-section pair mandrel x^*And y^*The moment of inertia and the product of inertia of,theta is centroid inertial principal axis and centroid axis x^*The included angle of (a).

Compared with the prior art, the invention has the following beneficial technical effects:

in the process of solving the section stress, the crack sections are dispersed into surface units, and the rigidity of the rotor with any section shape can be calculated, but the existing method is only limited to the calculation of the rigidity of the rotor with the circular section; in the stiffness calculation process, the iterative step of restoring force correction is adopted, so that the restoring force direction of the rotor in any rotor vortex difference angle state can be determined. In the existing model based on the assumption of dominant gravity, the restoring force direction is only assumed to be the reverse direction of gravity. Therefore, the method breaks through the constraint of the gravity dominance assumption and is suitable for solving the non-gravity dominance crack rotor respiratory stiffness. The obtained crack rotor respiration stiffness calculation model provides support for researching the vibration characteristics of the crack rotor and provides a theoretical basis for realizing the monitoring of the crack of the rotor of the large-scale rotating equipment.

Furthermore, rigidity fitting is carried out through a triangular basis function, so that crack rotation dynamics modeling and solving are facilitated, and the calculation accuracy and efficiency are improved.

Drawings

Fig. 1 is a flowchart of a stiffness respiratory function calculation method according to an embodiment of the present invention.

FIG. 2 is a schematic cross-sectional view of a crack described in an example of the present invention.

FIG. 3 is a schematic representation of the iterative process coordinate system and principal axes of inertia as described in the examples of the present invention.

Fig. 4a is a schematic closed cross-sectional view of an example of the present invention at a vortex slip angle ψ of 0.

Fig. 4b is a schematic closed cross-section view of an example of the invention at a transition vortex difference angle ψ pi/8.

Fig. 4c is a schematic closed cross-sectional view of an example of the present invention at a transition vortex difference angle ψ pi/4.

Fig. 4d is a schematic closed cross-section at a transition vortex difference angle ψ of 3 pi/8 in an example of the invention.

Fig. 5 is a schematic diagram of the solution of the three-way relative respiratory stiffness at a crack relative depth μ of 0.2 in the example of the present invention.

FIG. 6 is a plot of a fit function of x-direction relative stiffness versus discrete data for an example of the present invention.

FIG. 7 is a plot of a y-direction relative stiffness fit function versus discrete data for an example of the present invention.

FIG. 8 is a cross-relative stiffness fit function versus discrete data plot for an example of the present invention.

Detailed Description

The present invention will now be described in further detail with reference to specific examples, which are intended to be illustrative, but not limiting, of the invention.

The method determines the closed section of the crack section in the state according to the tension-compression stress balance principle of the crack section of the crack rotor under any rotational vortex differential angle. And calculating the crack rotor rigidity in the state according to the closed section and by combining a bending theory. And solving the corresponding restoring force by utilizing the obtained rigidity in the state, wherein the restoring force is used for correcting the tensile and compressive stress area in the state and correcting the restoring force direction in the state. And obtaining the rotor stiffness corresponding to the rotor vortex difference angle convergence in the state through iterative calculation. Obtaining the rotor stiffness converged under all the rotational vortex difference angles by adopting the same method; and finally, fitting the rotor stiffness at the discrete state into a continuous expression of a triangular basis function by using a data fitting method to obtain a stiffness respiration function and obtain the dynamic stiffness of the non-gravity dominant crack rotor. The method has high solving precision and strong applicability, can accurately solve the rigidity of rotors with different section shapes and breathing cracks in any rotating state, and is convenient for dynamic modeling and response characteristic research of the crack rotors.

The method comprises the steps of correcting the direction and the size of restoring force of a rotor in a current rotating state by using an iteration method, determining stress tension-compression distribution of a crack section in the state through the restoring force to obtain a closed section, and calculating the rigidity of the rotor in the state by using a bending theory. And (5) correcting restoring force by calculating rigidity, and sequentially iterating until convergence.

And fitting the rigidity at the discrete state into a continuous expression of the triangular basis function by using a data fitting method.

The basic flow of the method in the implementation of the invention is shown in fig. 1, the position of the principal axis of inertia and the direction of the restoring force are solved in sequence through three iterations, and all the vortex differential angles are traversed to obtain the three-way rigidity discrete value of any vortex differential angle, and the specific steps are as follows.

Firstly, as shown in FIG. 2, the crack section is composed of ① crack-free zone, ② crack-closed zone and ③ crack-open zone, V is normal direction of the crack, d is crack depth, o-xy is fixed coordinate system, o ' -x ' y ' is rotating coordinate system, theta is corner,is the swirl angle, psi is the rotational swirl difference angleAnd F is the restoring force. The crack cross section is divided into a plurality of discrete rectangular face units. At a certain vortex difference angle ψ, it is assumed that the elastic restoring force F is in the opposite direction of x'. FIG. 3, o_IIs the centroid, x, of a closed section (region ① + ②)^*And y^*Is a shaped axis parallel to x 'and y', x_IAnd y_IThe centroid is the main axis of inertia,is x_IAnd x^*The included angle of (a). Assuming main axis of inertia x_I,y_ICoinciding with the x 'y' axis, in this state, the closed section, i.e. the region ① + ② in fig. 2, is obtained by balancing the tensile and compressive stresses.

Then by the following formula

The centroid and centroid inertial axis position of the closed section, x in FIG. 3, is calculated_I,y_IAs shown. The calculation is iterated according to the flow shown in fig. 1 until the position of the principal axis of inertia converges. After the convergent centroid inertia axis is obtained, the stiffness of the centroid inertia axis in each axial direction under the coordinate system o ' -x ' y ' is calculatedAnd

by the following transformation

Obtaining the rigidity k of the material in each axial direction under a fixed coordinate system_x,k_yAnd k_xy. Wherein,is the rotor's angle of swirl in a fixed coordinate system o-xy, as shown in fig. 3.

The restoring force is corrected by the stiffness in this rotor state:and iterating until the restoring force converges. Thus, the rotor stiffness at the rotor differential angle is obtained. And calculating any other rotational vortex difference angle, and repeating the process to obtain the crack rotor stiffness at each discrete rotational vortex difference angle.

In the data fitting process, firstly, according to the symmetry of data, the expressions of the three-directional rigidity are respectively

k_x＝a(1)cos(ψ)+a(2)cos(ψ)²+a(3)cos(ψ)³+a(4)cos(ψ)⁶+a(5)cos(ψ)⁷+a(6)cos(ψ)⁸+a(7)cos(ψ)⁹+a(8)cos(ψ)¹⁰+a(9)cos(ψ)¹⁵+a(10)cos(ψ)¹⁹+a(11)

$k_y=\underset{i=1}{\overset{7}{\Σ}}a\left(i\right)cos{\left(\mathrm{\ψ}\right)}^i+a\left(8\right)$

$k_{xy}=\underset{i=1}{\overset{5}{\Σ}}a\left(i\right)sin\left(i\mathrm{\ψ}\right)+\underset{i=6}{\overset{10}{\Σ}}a\left(i\right)cos\left(\left(i-5\right)\mathrm{\ψ}\right)+a\left(11\right)$

And then solving by using a least square data fitting method to obtain a fitting coefficient a (i), namely a rigidity breathing function, so as to obtain a rigidity value of the breathing change.

In the preferred embodiment, the crack cross section shown in fig. 2 is taken as an example, where the relative depth of the crack is μ ═ D/D ═ 0.2, D is the crack depth, and D is the cross-sectional diameter. The rotational vortex difference angle was equally divided into 128 parts in the domain [0,2 π ], and the crack cross-section was discretized into a face unit with a side length of 1/100. According to the basic flow shown in fig. 1, the closed cross section of the crack section at different vortex difference angles is calculated by writing a MATLAB program, and fig. 4 shows the calculation results of the closed cross section at 4 vortex difference angles ψ 0, ψ π/8, ψ π/4 and ψ 3 π/8. Fig. 5 shows the result of the settlement of the discrete values of the 3-way stiffness at different rotor eddy difference angles, i.e. the stiffness-respiration function of the crack rotor, the stiffness value of which is respiration-variable. After data fitting, obtaining a continuous respiration function expression described by a triangular basis function as follows:

$\begin{array}{c}{k}_{x^{*}}=-0.0896\cos \left(\mathrm{\ψ}\right)+0.121\cos {\left(\mathrm{\ψ}\right)}^2-0.0295\cos {\left(\mathrm{\ψ}\right)}^3-0.145\cos {\left(\mathrm{\ψ}\right)}^6+0.026\cos {\left(\mathrm{\ψ}\right)}^7+0.175\cos {\left(\mathrm{\ψ}\right)}^8\\ -0.41\cos {\left(\mathrm{\ψ}\right)}^9-0.028\cos {\left(\mathrm{\ψ}\right)}^{10}+0.941\cos {\left(\mathrm{\ψ}\right)}^{15}-0.535\cos {\left(\mathrm{\ψ}\right)}^{19}+0.785\end{array}$

$\begin{array}{c}{k}_{y^{*}}=-0.26\cos \left(\mathrm{\ψ}\right)-0.17\cos {\left(\mathrm{\ψ}\right)}^2-0.006\cos {\left(\mathrm{\ψ}\right)}^3-0.034\cos {\left(\mathrm{\ψ}\right)}^4-0.04\cos {\left(\mathrm{\ψ}\right)}^5\\ +0.01\cos {\left(\mathrm{\ψ}\right)}^6+0.01\cos {\left(\mathrm{\ψ}\right)}^7+0.9\end{array}$

$\begin{array}{c}{k}_{x^{*}{y}^{*}}=-0.14\sin \left(\mathrm{\ψ}\right)-0.09\sin \left(2\mathrm{\ψ}\right)-0.015\sin \left(3\mathrm{\ψ}\right)-0.006\sin \left(4\mathrm{\ψ}\right)-0.005\sin \left(5\mathrm{\ψ}\right)\\ +0.00006\cos \left(\mathrm{\ψ}\right)-0.000017\cos \left(2\mathrm{\ψ}\right)-0.00008\cos \left(3\mathrm{\ψ}\right)-0.00006\cos \left(4\mathrm{\ψ}\right)\\ +0.000006\cos \left(5\mathrm{\ψ}\right)+0.000045\end{array}$

to illustrate the accuracy of the data fitting, fig. 6-8 compare the 3-way stiffness raw discrete data and the fitted data difference using relative stiffness, where K_maxRepresenting the bending stiffness of a normal rotor. As can be seen from the figure, the curve of the fitting function almost coincides with the discrete stiffness points, indicating that the error of the fitting function is within the allowable range.

Claims (4)

1. A non-gravity dominant crack rotor stiffness respiration function calculation method is characterized in that firstly, at a certain vortex difference angle, stress tension and compression distribution of a rotor cross section at a crack under the state is determined by assuming an initial restoring force direction, and a crack section closed area is obtained; wherein, the cross section of the rotor at the crack is the crack section; the crack section closed area is a closed section; then, calculating to obtain the rigidity of the rotor in the state by using a bending theory; correcting restoring force in the state by calculating the rigidity of the rotor, and sequentially iterating until the rigidity is converged; traversing all discrete values in the vortex difference angle definition domain to obtain a corresponding discrete stiffness value; and finally, fitting the rigidity in the discrete state into a triangular basis function continuous expression about the rotor eddy difference angle by using a data fitting method, and obtaining the rigidity of the non-gravity dominant crack rotor changing in a breathing state.

2. The method for calculating the rotor stiffness breathing function of the non-gravity dominant crack according to claim 1, is characterized by comprising the following steps:

1) dispersing a domain [0,2 pi ] of the vortex difference angle psi into N equal parts, and dividing the crack section into a plurality of discrete surface units;

2) under a certain selected rotor vortex difference angle, assuming that the restoring force direction and the rotor axis displacement direction are opposite, calculating the stress of each surface unit, and calculating to obtain a closed section according to a neutral layer theory;

3) calculating the position of an inertia main shaft of the current closed section, repeating the step 2, and recalculating the stress of the surface unit of the current closed section until a convergent closed section is obtained;

4) calculating the rigidity of the rotor in the state by using a bending theory through the convergent closed section;

5) calculating the rotor restoring force by the following formula by calculating the rotor stiffness under the rotor vortex difference angle:

$\left[\begin{array}{c}{F}_x\\ F_y\end{array}\right]=-\left[\begin{array}{cc}{k}_x& k_{xy}\\ k_{xy}& k_y\end{array}\right]\left[\begin{array}{c}{x}_{o^{\′}}\\ y_{o^{\′}}\end{array}\right];$

wherein, F_xAnd F_yComponent of restoring force along x and x axes in a fixed coordinate system o-xy, k_x，k_yAnd k_xyThree-dimensional flexural stiffness, x, of the rotor relative to a fixed coordinate system_o'And y_o'Is the axial displacement of the rotor relative to a fixed coordinate system;

6) repeating the steps 2) -5) until the restoring force in the state is converged, and obtaining the rotor rigidity under the rotor vortex difference angle;

7) selecting any non-calculated rotor vortex difference angle, and repeating the steps 2) to 6) until traversing N rotor vortex difference angles in the defined domain; obtaining the corresponding rotor stiffness under each rotational vortex difference angle;

8) and fitting the N rotor rigidities in the discrete rotor eddy difference angle state into a continuous expression of a triangular basis function, namely a rigidity breathing function, by using a data fitting method, so as to obtain the rigidity value of the breathing change.

3. The method for calculating the rotor stiffness respiration function of the non-gravity dominant crack according to the claim 2, wherein the method for fitting the N rotor stiffnesses in the discrete states into the triangular basis function continuous expression about the rotor vortex difference angle by using the data fitting method comprises the following steps:

a step of making the axisymmetrically distributed horizontal and vertical stiffness in the form of a cos (ψ) basis function, wherein,

k_x＝a(1)cos(ψ)+a(2)cos(ψ)²+a(3)cos(ψ)³+a(4)cos(ψ)⁶+a(5)cos(ψ)⁷+a(6)cos(ψ)⁸+a(7)cos(ψ)⁹+a(8)cos(ψ)¹⁰+a(9)cos(ψ)¹⁵+a(10)cos(ψ)¹⁹+a(11)

$k_y=\underset{i=1}{\overset{7}{\Σ}}a\left(i\right)cos{\left(\mathrm{\ψ}\right)}^i+a\left(8\right);$

a step of making the cross stiffness of the centrosymmetric distribution into a Fourier series form, wherein,

$k_{xy}=\underset{i=1}{\overset{5}{\Σ}}a\left(i\right)sin\left(i\mathrm{\ψ}\right)+\underset{i=6}{\overset{10}{\Σ}}a\left(i\right)cos\left(\right(i-5\left)\mathrm{\ψ}\right)+a\left(11\right);$

in the formula 3, a (i) is a undetermined coefficient, and a vortex shedding difference angle psi belongs to [0,2 pi ];

and obtaining coefficients a (i) of the stiffness values under the N discrete vortex difference angles according to a least square data fitting principle by using the obtained stiffness values under the N discrete vortex difference angles, and determining a continuous expression of the stiffness respiratory function.

4. The method for calculating the rigidity breathing function of the non-gravity dominant crack rotor according to the claim 2, characterized in that the centroid and the centroid inertia principal axis position of the current closed section are calculated by the following formula,

$x_{o_I}=\frac{\underset{i=1}{\overset{n}{\Σ}}{x}_i{A}_i}{\underset{i=1}{\overset{n}{\Σ}}{A}_i},y_{o_I}=\frac{\underset{i=1}{\overset{n}{\Σ}}{y}_i{A}_i}{\underset{i=1}{\overset{n}{\Σ}}{A}_i},\tan \left(2\mathrm{\θ}\right)=-\frac{2I_{x^{*}{y}^{*}}}{I_{y^{*}}-I_{x^{*}}};$

in the formula, A_iThe area of the inner face unit of the closed section is shown, and n is the number of the inner face units of the closed section;androtating the centroid coordinate under the coordinate system o ' -x ' y ',andclosed cross-section pair mandrel x^*And y^*Theta is the centroid principal axis of inertia and the centroid axis x^*The included angle of (a).