Dynamic modeling method for rotating blade with breathing effect and variable cross-section crack

Abstract

This invention provides a dynamic modeling method for a rotating airfoil with a variable cross-section cracked blade exhibiting a breathing effect, belonging to the field of mechanical dynamics. This method considers the influence of the pre-torsion angle on the potential energy of a micro-segment straight beam element, introduces a potential energy correction term, and models a healthy airfoil blade based on finite element beam element theory. Then, based on the strain energy release rate and Castingliano's theorem, a cracked beam element is introduced. To simulate crack breathing, the contact state of the crack surface is determined based on the stress value at the crack surface, establishing a dynamic model of the rotating airfoil blade considering the crack breathing effect. The effectiveness of the model is verified by comparison with rotating cracked blade models established using ANSYS solid elements and contact elements. Finally, the effects of different crack parameters and load parameters on the nonlinear vibration characteristics of the cracked blade are investigated.

Description

Technical Field

[0001] This invention relates to the field of mechanical dynamics, and more particularly to a dynamic modeling method for a rotating airfoil with a variable cross-section cracked blade exhibiting a breathing effect. Background Technology

[0002] Rotating blades, as one of the most critical components of aero-engines, directly influence the engine's normal operating performance through their structure and condition. However, due to potential material and manufacturing defects in the blades themselves, and the complex cyclic loads imposed on the blades during engine operation (frequent starting, acceleration, deceleration, and shutdown), these factors can easily lead to cracks or even fractures in the blades. Blade fracture can cause engine failure, affecting mission completion, or even serious accidents. Therefore, dynamic modeling of cracked blades and studying their vibration response characteristics are of great significance to the dynamics research and development of aero-engines.

[0003] In research on dynamic modeling of cracked single blades, some scholars have used shell elements and solid elements to model cracked blades. However, shell element modeling is relatively complex, while solid element modeling involves large computational costs and low solution efficiency. Considering breathing cracks makes the modeling process even more difficult. Therefore, few scholars have modeled rotating blades with breathing cracks based on shell elements and solid elements. The geometry of real blades is quite complex. Therefore, many scholars simplify the blade as a cantilever beam with a rectangular cross-section and then introduce cracks for research. However, there are few studies on blades with realistic shapes. Modeling real airfoil blades based on beam elements is inaccurate and cannot accurately account for the crack nonlinearity caused by crack breathing due to the airfoil's variable cross-section, as well as the efficiency of crack nonlinear modeling. Summary of the Invention

[0004] In view of this, the purpose of this invention is to propose a dynamic modeling method for a rotating airfoil variable cross-section cracked blade with breathing effect, so as to solve the problems in the existing dynamic modeling of rotating cracked blades, such as the lack of a dynamic modeling method for airfoil variable cross-section cracked blades with breathing effect and the lack of equivalent use of cracked beam elements to represent the breathing effect of cracked blades.

[0005] The technical means employed in this invention are as follows:

[0006] A dynamic modeling method for a rotating airfoil with a variable cross-section cracked blade exhibiting a breathing effect includes the following steps:

[0007] S1. Considering the effects of centrifugal stiffening, rotational softening, and Coriolis force, a dynamic model of a real blade rotating micro-segment beam element with variable cross-section is established based on finite element theory and Timoshenko beam theory.

[0008] S2. Based on the strain energy release rate and Castingliano's theorem, a cracked beam element is introduced into the dynamic model to establish a dynamic model of a rotating airfoil section blade that considers the crack breathing effect.

[0009] Furthermore, S1 specifically includes the following steps:

[0010] S11. Establish an airfoil blade model;

[0011] The airfoil blade is assumed to be a rotating Timoshenko cantilever beam, which is then divided into N elements. The k-th rotating Timoshenko straight beam element is taken as the research object, and its motion differential equation is established; where l k γ is the length of the k-th discrete finite element; β0 is the initial installation angle of the blade, and γ k Let β be the pre-torsion angle of the k-th Timoshenko beam element, and let β be the installation angle of the k-th Timoshenko beam element. k =β0+γ k α0 is the blade rotational angular displacement;

[0012] S12. Determine the blade cross-section parameters;

[0013] The parameters of each section are determined based on the geometric properties of the actual cross-section; the cross-section of the k-th element is defined, first by describing the coordinate x... k The closed curves that form the profile of this cross section are used to parametrically represent the function expressions of the blade tip and blade back through spline fitting; the expression of the cross section area, the coordinate transformation matrix of the local coordinate system and the element coordinate system, the expression of the centroid coordinates of the cross section, the moment of inertia of the cross section about the axis, the torsional stiffness of the airfoil section, and the additional torsional moment of inertia of the cross section caused by the pre-torsion angle are obtained.

[0014] S13. Considering the influence of the pre-twist angle on the potential energy of the micro-segment straight beam element, an energy correction term is introduced. After correction, the total potential energy U of the rotating variable cross-section micro-segment straight beam is obtained. e The calculation expression is obtained; then, combined with the kinetic energy of the rotating variable cross-section micro-segment straight beam, the dynamic model of the rotating micro-segment beam element is obtained.

[0015] Furthermore, in S11, the displacement of a point P on any section of the blade in the three directions in the global coordinate system can be expressed as:

[0016]

[0017] Where, θ x θ y and θ z To describe the cross section around x of an arbitrary element e axis, y e axis and ze The second kind of Euler angle is introduced by the rotational angular displacement of the axis; x, y, z are the coordinates of point P in the local coordinate system; R k The origin of the unit coordinate system is o e In the rotating coordinate system OX r Y r Z r The coordinates in the x-direction below; u, v, and w are the coordinates of the centroid of the element section along the element coordinate system o. e x e Axis, o e y e axis and o e z e Displacement of the axis.

[0018] Furthermore, in S12, the torsional stiffness J of the airfoil section... s The empirical expression is:

[0019]

[0020] Where A is the cross-sectional area, I max and I min These are the maximum and minimum moments of inertia of the cross section, respectively; the additional torsional moment of inertia J of the cross section due to the pre-torsion angle γ(L). a The formula for calculation is:

[0021]

[0022] Furthermore, in S13, after correction, the total potential energy U of the rotating variable cross-section micro-segment straight beam is obtained. e The calculation expression is:

[0023]

[0024] Among them, U bte and U s Let γ' be the bending-pendulum-axis-torsion coupling potential energy and shear potential energy of the k-th rotating variable cross-section micro-segment straight beam, respectively, and γ′ be the rate of change of the torsion angle. F b The dynamic model of the actual blade micro-segment beam element is obtained by representing the basic excitation load applied along the blade axis and simplifying the process:

[0025]

[0026] Among them, M e The mass matrix of a micro-segment beam element in the element coordinate system; G is the structural stiffness matrix of the Timoshenko beam element; e Here is the Coriolis force matrix for the beam element. Here is the centrifugal stiffening matrix of the beam element. For rotation softening matrix; The element angular acceleration stiffness matrix is ​​given by... Let angular acceleration be expressed, then it satisfies δ e and F e These represent the element node displacement and element node force vector, respectively.

[0027] Furthermore, the dynamic model of the cracked rotating airfoil blade includes healthy beam elements and cracked beam elements, and S2 specifically includes the following steps:

[0028] S21. First, a dynamic model of a healthy rotating airfoil section blade is established based on beam element theory. Then, the healthy beam elements on the model are replaced with cracked beam elements to achieve the purpose of implanting cracks.

[0029] S22. By measuring the vibration displacement of the crack surface at various times, calculate the bending moment at the crack surface, calculate the stress on the crack surface, and determine the contact state of the crack surface based on the stress value at the crack surface, thereby simulating crack breathing.

[0030] S23. The external load on the blade was simulated using aerodynamics, and a dynamic model of a rotating airfoil blade considering the crack breathing effect was established.

[0031] Furthermore, S21 specifically includes:

[0032] Total strain energy U of cracked beam element m for:

[0033] U m =U n +U c

[0034] Among them, U n The strain energy of a crack-free beam element; U c This refers to the loss of strain energy caused by cracks;

[0035] U n and U c These can be represented as follows:

[0036]

[0037]

[0038] Among them, A x and A c These represent the cross-sectional area of ​​the unit and the crack surface area, respectively; K Ir K IIr and K IIIr , are the stress intensity factors for type I, type II and type III cracks, respectively, r = 1, 2, ..., 6;

[0039] The expression for the stress intensity factor is:

[0040]

[0041]

[0042]

[0043] Among them, F1(a), F2(a), F II (a) and F III (a) is the stress intensity factor correction coefficient for the crack; a is the crack depth, h c The cross-sectional thickness is along the crack direction;

[0044] The nodal displacements of the cracked beam element are divided into the displacements of the healthy beam element and the displacements of the cracked beam element. The nodal displacements are as follows:

[0045]

[0046] When the crack opens, the structural stiffness matrix of the cracked beam element is:

[0047]

[0048] Among them, f nc and f c These are the compliance matrix of a healthy beam element and the additional compliance matrix caused by the introduction of cracks, respectively; T c e This is the transformation matrix.

[0049] Furthermore, S22 specifically includes:

[0050] A crack breathing function is established to determine the breathing state. Based on the vibration displacement of the crack surface at various times, the bending moment at the crack surface is calculated, and then the stress on the crack surface is calculated. The contact state of the crack surface is determined based on the stress value at the crack surface, and the bending deformation caused by the crack leads to the crack surface circumferential... e z e and o e y e and o e y e Bending moment M of the shaft z and M y They are represented as follows:

[0051]

[0052] in, and For the o e z e and oe y e Angular displacement of the shaft;

[0053] The structural stiffness matrix of the m-th healthy beam is replaced with the stiffness matrix of the cracked beam to obtain the structural stiffness matrix of the rotating cracked blade.

[0054] The crack depth ratio α is characterized by the proportion of the cracked region area to the total cross-sectional area, i.e., α = S1 / (S1 + S2); the crack location ratio δ is introduced, assuming the cracked blade consists of n parts. all The segment beam is composed of the nth segment. crack If the segment beam is a cracked beam, then:

[0055]

[0056] Furthermore, S23 specifically includes:

[0057] To simulate the external loads experienced by compressor blades during operation, an aerodynamic force F is introduced. p The aerodynamic force on each beam segment is The expression is as follows:

[0058]

[0059] Among them, A k Let A be the area of ​​the k-th unit that bears the aerodynamic force. p f r t and t represent the aerodynamic force amplitude, excitation frequency, and time acting on the rotating cracked blade, respectively; N v1 N v2 N v3 and N v4 These are four shape functions representing the bending displacement direction of the Timoshenko beam element.

[0060] Finally, the dynamic model of the rotating cracked blade is obtained:

[0061]

[0062] Where M is the mass matrix of the rotating cracked blade; K e G is the structural stiffness matrix of the rotating cracked blade; G is the Coriolis force matrix of the rotating cracked blade, and K is the structural stiffness matrix of the rotating cracked blade. c K is the centrifugal stiffening matrix of the rotating cracked blade. s K is the rotation softening matrix; acc The angular acceleration stiffness matrix of the rotating cracked blade is given by... Let angular acceleration be expressed, then it satisfies δ and F represent the nodal displacement and nodal force vector of the rotating cracked blade, respectively; F p The aerodynamic force on the rotating cracked blade.

[0063] Compared with the prior art, the present invention has the following advantages:

[0064] The present invention provides a realistic airfoil blade modeling method that discretizes the airfoil blade into rotating Timoshenko cantilever beams and introduces a correction term to account for the influence of the pre-torsion angle on the additional torsional moment of inertia and potential energy of the micro-segment straight beam element, thereby achieving accurate modeling of blades with complex geometries.

[0065] The method for modeling real airfoil blades using beam elements provided by this invention reduces the modeling difficulty of real airfoil blades by using Timoshenko beam theory, and achieves efficient modeling of blades with airfoil cross sections.

[0066] The present invention provides a method for modeling real airfoil cracked blades, which calculates the stress on the crack surface by the magnitude of the bending moment at the crack surface, and judges the contact state of the crack surface based on the stress value at the crack surface, thus realizing efficient modeling of airfoil section cracked blades that take into account the crack breathing effect.

[0067] Traditional Timoshenko beam elements are relatively accurate in calculating the bending modes of beams with arbitrary cross-sections. However, when calculating torsional modes, they only provide relatively accurate solutions for beams with regular cross-sections (rectangular, circular, etc.). For beams with irregular cross-sections such as airfoils, the torsional error between Timoshenko beam elements and solid elements is significant. This is mainly because the traditional Timoshenko beam model ignores the warping of the cross-section, resulting in inaccurate torsional calculations. To accurately predict torsion, this invention introduces a torsional correction factor. For each airfoil cross-section in the blade, it is stretched axially into a straight beam of unequal length. A corresponding solid model is established based on the finite element software ANSYS. Using the solid model as a reference, the torsional coefficient of the model is corrected. Finally, the torsional correction factor of the beam at each length is spline-fitted to approximately obtain the torsional correction factor of any micro-segment beam. Attached Figure Description

[0068] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0069] Figure 1 This is a schematic diagram of the finite element discrete model of the blade of the present invention;

[0070] Figure 2 The distance x from the leaf root in this invention k A schematic diagram of the cross-section of the k-th beam element at point k;

[0071] Figure 3 This is a comparison chart of the dynamic frequency curves of healthy blades in this invention;

[0072] Figure 4 Comparison of vibration mode cloud diagrams for the healthy blades of this invention; (a) is the vibration mode diagram of the solid element model, and (b) is the vibration mode diagram of the beam element model in this paper;

[0073] Figure 5 The following are dynamic frequency curves of the first four frequencies of the blade of the present invention under different installation angles: (a) is a curve showing the change of the first bending natural frequency of the blade with the installation angle; (b) is a curve showing the change of the first torsional natural frequency of the blade with the installation angle; (c) is a curve showing the change of the second bending natural frequency of the blade with the installation angle; and (d) is a curve showing the change of the second torsional natural frequency of the blade with the installation angle.

[0074] Figure 6 The following are dynamic frequency curves of the first four frequencies of the blade of the present invention under different pre-twist angles: (a) is the curve of the first bending natural frequency of the blade as a function of the pre-twist angle; (b) is the curve of the first torsional natural frequency of the blade as a function of the pre-twist angle; (c) is the curve of the second bending natural frequency of the blade as a function of the pre-twist angle; and (d) is the curve of the second torsional natural frequency of the blade as a function of the pre-twist angle.

[0075] Figure 7 This is a schematic diagram of a cracked beam unit containing transverse penetrating cracks according to the present invention;

[0076] Figure 8 This is a schematic diagram of the crack cross-section of the present invention;

[0077] Figure 9 The following are comparison diagrams of the blade tip response vibration response of the present invention: (a) Comparison diagram of blade tip response vibration response under superharmonic resonance; (b) Comparison diagram of blade tip response vibration response under first-order resonance.

[0078] Figure 10 (a) is a three-dimensional spectrum diagram of different aerodynamic amplitudes under superharmonic resonance; (b) is a three-dimensional spectrum diagram of different aerodynamic amplitudes under first-order resonance.

[0079] Figure 11 The following are comparison diagrams of harmonic amplitude under superharmonic resonance conditions according to the present invention: (a) is a comparison diagram of harmonic amplitude under 0 harmonic resonance conditions; (b) is a comparison diagram of harmonic amplitude under 1 harmonic resonance conditions; and (c) is a comparison diagram of harmonic amplitude under 2 harmonic resonance conditions.

[0080] Figure 12The following are comparison diagrams of the harmonic amplitude under the first harmonic resonance state of the present invention: (a) is a comparison diagram of the harmonic amplitude under the first harmonic resonance state with 0 harmonics; (b) is a comparison diagram of the harmonic amplitude under the first harmonic resonance state with 1 harmonics; and (c) is a comparison diagram of the harmonic amplitude under the first harmonic resonance state with 2 harmonics. Detailed Implementation

[0081] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0082] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0083] This invention provides a method for dynamic modeling of a rotating airfoil with a variable cross-section crack containing a breathing crack, comprising the following steps:

[0084] S1. Considering the effects of centrifugal stiffening, rotational softening, and Coriolis force, a dynamic model of the rotating micro-segment beam element was established based on finite element theory and Timoshenko beam theory. The specific process is as follows:

[0085] S11. Establish an airfoil blade model;

[0086] The airfoil blade is idealized as a rotating Timoshenko cantilever beam and divided into N units. Figure 1 The diagram illustrates a discrete finite element model of the blade. We will now focus on the k-th rotating Timoshenko straight beam element and establish its equations of motion. The length of the k-th discrete finite element is l. k ;β0 and γ kLet represent the initial installation angle of the blade and the pre-torsion angle of the k-th Timoshenko beam element, respectively, and let β represent the installation angle of the k-th Timoshenko beam element. k =β0+γ k α0 is the angular displacement of the blade rotation.

[0087] The displacement of a point P on any section of the blade in the three directions in the global coordinate system can be expressed as:

[0088]

[0089] Where, θ x θ y and θ z To describe the cross section around x of an arbitrary element e axis, y e axis and z e The second kind of Euler angle is introduced by the rotational angular displacement of the axis; x, y, z are the coordinates of point P in the local coordinate system; R k The origin of the unit coordinate system is o e In the rotating coordinate system OX r Y r Z r The coordinates in the x-direction below; u, v, and w are the coordinates of the centroid of the element section along the element coordinate system o. e x e Axis, o e y e axis and o e z e Displacement of the axis.

[0090] S12. Determine the blade cross-section parameters;

[0091] The moment of inertia and area of ​​each section of a real blade are different, and the centroid and bending center of the airfoil section are generally not in the same position. Therefore, it is necessary to determine the parameters of each section based on the geometric properties of the actual section.

[0092] To define the cross section of the k-th element, we first need to describe the coordinates x and y. k The closed curves forming this cross-sectional profile, since the geometric input is given by a set of points, can have their functional expressions for the leaf base and leaf back parametrically represented by spline fitting. Now, using fourth-order polynomial interpolation, the leaf base and leaf back profiles are expressed as follows:

[0093]

[0094] Where a0, a1, ..., a4 and a0', a1', ..., a4' are all constants, and the blade cross-sectional profile can be approximated using the above formula:

[0095] y(z) = y1(z) - y2(z)

[0096] The cross-sectional area can be expressed as:

[0097] A = A dA=∫ A dydz

[0098] Figure 2 x is the distance from the leaf root k A schematic diagram of the cross-section of the k-th beam element at point o, with the element coordinate system o. e z e y e coordinate axis o e y e and o e z e All axes are principal axes of inertia of the cross-sectional shape, and the origin is the centroid; coordinate system o f y f z f The two coordinate axes o f y f and o f z f Relative to the coordinate axis o e y e and o e z e Parallel; coordinate axis o f y f and o e y e The distance is represented by dy, and the coordinate axis o f z f and o e z e The distance is represented by dz; the coordinate system is o. f y f z f The origin o f That is, the shear center of the cross section; Let T be the angle between the coordinate axes of the local coordinate system and the principal inertial axis (element coordinate system) of the blade. The coordinate transformation matrix T between the local coordinate system and the element coordinate system can be expressed as:

[0099]

[0100] Centroid of unit section o e Along y e axis and z e The linear displacement components v and w of the axis and the shear center o of the element section f Along y e axis and z e The relationship between the linear displacement components v1 and w1 of the shaft is as follows:

[0101]

[0102] The entire cross section for o e y e axis and o e z e The static moment of the axis, i.e., the first-order moment of inertia, can be written as:

[0103]

[0104] The expression for the centroid coordinates of the cross section can be obtained:

[0105]

[0106] Section to o e y e Axis, o e z e The second moment of inertia I of the axis y I z Inertial product I yz and the polar moment of inertia I of the cross section x They are expressed as follows:

[0107]

[0108] Except for a few regular and simple cross-sections (rectangular, circular, annular, etc.), the torsional moment of inertia of a cross-section can only be approximated or measured experimentally for other complex cross-sections. The torsional stiffness J of an airfoil cross-section... s The empirical expression is:

[0109]

[0110] Furthermore, for pre-twisted blades, the additional torsional moment of inertia J of the cross section caused by the pre-twisting angle γ(L) a The calculation formula can be expressed as follows:

[0111]

[0112] The torsional moment of inertia J of the blade was preliminarily calculated. t =J s +J a Here, an approximate algorithm is used to determine the upper and lower bounds of the torsional moment of inertia to verify the correctness of the empirical formula. First, the torsional stiffness is represented by the warping function ψ:

[0113]

[0114] Where ψ satisfies the boundary value problem of the Laplace equation:

[0115]

[0116] make:

[0117]

[0118]

[0119] Therefore, we get:

[0120]

[0121] so,

[0122] J t =I(ψ)-I1(ψ)=I(ψ)

[0123] The total potential energy U of the rotating variable cross-section micro-segment straight beam e The calculation expression is as follows:

[0124] U e =U bte +U s

[0125] Among them, U bte Let U be the bending-swinging-axis-torsional coupling potential energy of the k-th rotating variable cross-section micro-segment straight beam. s Its shear potential energy.

[0126] The bending-pendulum-axis-torsion coupling potential energy and shear potential energy of the k-th segment of the rotating variable cross-section straight beam can be expressed as follows:

[0127]

[0128]

[0129] Where u1, v1, and w1 represent the shear centers of the element cross section, respectively. f Along x e axis, y e axis and z e The linear displacement components of the axis; A represents the cross-sectional area, F c In the unit coordinate system o e x e y e z e Below, at a distance of o from the origin. e The centrifugal force experienced by the micro-segment of the straight beam at a distance x; E and G represent the elastic modulus and shear modulus of the blade, respectively; κ y and κ z Let be the shear correction factors for the beam section along the principal axes of inertia y and z, respectively; the expression for C0 is:

[0130]

[0131] S13. Introduce an energy correction term;

[0132] Furthermore, considering the influence of the pre-twist angle on the potential energy of the micro-segment straight beam element, a correction term is introduced. After correction, the total potential energy U of the rotated variable cross-section micro-segment straight beam is... e The calculation expression can be rewritten as:

[0133]

[0134] Where γ′ is the rate of change of the twist angle. F b This indicates the basic excitation load applied along the blade axis.

[0135] The kinetic energy calculation formula for a rotating micro-segment beam element is as follows:

[0136]

[0137] in,

[0138]

[0139] Finally, the following system of second-order linear differential equations with constant coefficients for the micro-segment beam element is obtained:

[0140]

[0141] Among them, M e The mass matrix of a micro-segment beam element in the element coordinate system; G is the structural stiffness matrix of the Timoshenko beam element; e Here is the Coriolis force matrix for the beam element. Here is the centrifugal stiffening matrix of the beam element. For rotation softening matrix; The element angular acceleration stiffness matrix is ​​given by... Let angular acceleration be expressed, then it satisfies δ e and F e These represent the element node displacement and element node force vector, respectively.

[0142] Traditional Timoshenko beam elements are relatively accurate in calculating the bending modes of beams with arbitrary cross-sections. However, when calculating torsional modes, they only provide relatively accurate solutions for beams with regular cross-sections (rectangular, circular, etc.). For beams with irregular cross-sections such as airfoils, the torsional error between Timoshenko beam elements and solid elements is significant. This is mainly because the traditional Timoshenko beam model ignores the warping of the cross-section, resulting in inaccurate torsional calculations. To accurately predict torsion, a torsional correction factor is introduced. For each airfoil cross-section in the blade, it is stretched axially into a straight beam of unequal length. A corresponding solid model is established using the finite element software ANSYS. The torsional coefficient of the model is corrected based on the solid model. Finally, the torsional correction factor of the beam at each length is spline-fitted to approximately obtain the torsional correction factor of any micro-segment beam.

[0143] The effectiveness of the realistic blade model built using beam elements in this chapter is verified by comparing its natural frequencies with those of a finite element model established using the commercial finite element software ANSYS. The material parameters are as follows: density ρ = 4370 kg / m³, elastic modulus E = 125 GPa, and Poisson's ratio μ = 0.3. In ANSYS modeling, Solid185 elements were used to extract data from points on each cross-section, fit a contour curve, and then generate a solid model. When solving for the natural frequencies, for the beam element model, six degrees of freedom of the corresponding nodes on the bottom surface were constrained; for the solid element model, all degrees of freedom of all nodes on the bottom surface were constrained.

[0144] The errors of the first four natural frequencies between the solid element model and the beam element model in this paper were finally obtained, as shown in Table 1. Table 1 shows that the errors between the model in this chapter and the ANSYS solid model are small, all within 5%, proving the effectiveness of the model in this chapter.

[0145] Table 1 Comparison of Inherent Frequency Errors

[0146]

[0147] Because the natural frequency of the blades is affected by effects such as centrifugal stiffening and rotational softening, the rotor blade rotational speed is a key element that needs to be simulated. Rotational speeds are applied to the blades in three cases, with the speed range set from 0 r / min to 10000 r / min and a speed increment Δn = 100 r / min. Figure 3 This is a comparison chart of the final dynamic frequency curves. (From...) Figure 3 It can be seen that the trends of the natural frequencies of the solid model of the blade and the beam model are in good agreement under different rotational speeds, and the natural frequencies of each order increase slightly with the increase of rotational speed.

[0148] The eigenvectors corresponding to the first four natural frequencies (static frequencies) are extracted to plot the blade's mode shape. The blade is modeled using 41 nodes and 40 elements, therefore, representing the blade profile requires 41 sets of profile point data, each set consisting of 42 points. In this paper, each node of the beam element has 6 degrees of freedom, so the mode shape vector corresponding to each natural frequency consists of 246 elements. Let the mode shape vector be u... m The x, y, and z coordinates of the points that constitute the contour of each cross section are x ij y ij z ij (i = 1, 2, 3, ..., 41; j = 1, 2, 3, ..., 42). Then the coordinates (x, y) of the points on the deformed cross-section profile are... ijd ,y ijd ,z ijd This can be expressed as follows:

[0149]

[0150] Let the desired mode shape order be m (m = 1, 2, 3, ...), then u is the mode shape vector u. m The m-th column vector.

[0151] The displacements Δu, Δv, and Δw of each point in the x, y, and z directions can be expressed as follows:

[0152]

[0153] The total deformation at each point is then...

[0154] Based on the above theory, the first four mode shapes of the airfoil blade based on beam elements are plotted in this chapter. Figure 4 This is a comparison of the contour plots of the first four vibration modes. Figure 4 It can be seen that the first four mode shape contour plots of the blade model built based on beam elements in this chapter are basically consistent with those of the blade model built using ANSYS solid elements, further proving the effectiveness of the model in this chapter and the correctness of the derivation of the formulas representing the magnitude of deformation at each node. Furthermore, based on the magnitude of the deformation, Figure 5 The airfoil blade model established for this chapter includes a comparison diagram of its shape before and after vibration deformation, combined with... Figure 4 It can be clearly seen that the first four vibration modes of the airfoil variable cross section blade model in this paper are the first bending mode, the first torsion mode, the second bending mode, and the second torsion mode.

[0155] When the blade rotates, it introduces Coriolis force, centrifugal stiffening, and rotational softening effects. The effects of the installation angle β and pre-twist angle γ on the first-order bending natural frequency Y of the blade were discussed at different rotational speeds. I The first-order torsional natural frequency θX,I The second-order bending natural frequency Y II and the second-order torsional natural frequency θ X,II The impact.

[0156] To discuss the influence of the installation angle on the first four dynamic frequencies, this section selects blade installation angles β = 0°, 10°, ..., 70°; the rotational speed n range is set to [0, 20000] r / min, and the rotational speed increment Δn = 5000 r / min. The pre-twist angle remains constant, and based on the angle between the principal axes of inertia of each blade cross-section profile during modeling, the pre-twist angle γ = 24° is selected. Figure 5 The dynamic frequency curves of the blade at different installation angles are shown for the first four frequencies. Figure 5 It can be seen that when the rotational speed is 0, changes in the installation angle do not affect the magnitude of the blade's natural frequency. When the rotational speed is non-zero, as the installation angle increases, the first and second order bending natural frequencies Y... I Y II and the first and second order torsional natural frequencies θ X,I θ X,II This also increases accordingly, and the larger the installation angle β, the greater the increment of each natural frequency. Figure 5 (a) and Figure 5 (b) It can be seen that the increments of each natural frequency are related to both the rotational speed and the mounting angle. Under the conditions of high rotational speed and large mounting angle, the first-order bending natural frequency Y... I and the first-order torsional natural frequency θ X,I The increment is particularly large. This is mainly due to the increase in the installation angle β, which leads to changes in both the bending direction (Y-direction) and the torsional direction (θ-direction). X The reduction of the rotational softening effect on the (upward) side.

[0157] To discuss the influence of the pre-twist angle on the first four dynamic frequencies, this invention selects blade pre-twist angles γ = 0°, 10°, ..., 70°; the rotational speed n range is set to [0, 20000] r / min, and the rotational speed increment Δn = 5000 r / min. The installation angle remains constant, and based on the angle between the principal axis of inertia of the blade bottom surface profile and the Y-axis during modeling, the installation angle β = 37.7° is taken. Figure 6 This shows the dynamic frequency curves of the first four frequencies of the blade at different pre-twist angles. From... Figure 6 It can be seen that increasing the pre-torsion angle γ(L) leads to an increase in the first-order bending natural frequency Y. I The first-order torsional natural frequency θ X,I and the second-order torsional natural frequency θ X,II The increase of , and the second-order bending natural frequency Y II It will decrease accordingly.

[0158] like Figure 6 (a)- Figure 6As shown in (c), the larger the pre-twist angle γ, the higher the first-order bending natural frequency Y. I The first-order torsional natural frequency θ X,I and the second-order bending natural frequency Y II The larger the increment, the greater; conversely, such as Figure 6 As shown in (d), the second-order torsional natural frequency θ X,II The increment will decrease as the pre-twist angle increases.

[0159] S2. Using the dynamic model of the airfoil variable cross-section blade established in S1, and based on the strain energy release rate and Castingliano's theorem, a cracked beam element was introduced to establish a dynamic model of the rotating airfoil cross-section blade considering the crack breathing effect. The specific process is as follows:

[0160] S21. Establish the cracked beam element model;

[0161] Now assume that the m-th beam element is a cracked beam element. Figure 7 This is a schematic diagram of a cracked beam element containing a transverse penetrating crack. The two end nodes of this beam element are i and j, respectively, and each node is subjected to forces in six directions. Among them, P1 and P7 are axial tensile forces, and P4 and P... 10 P2 represents the axial torque, P8, P3, and P9 represent the shear forces in the y and z directions, respectively, and P5, P... 11 P6 and P 12 These are the bending moments in the y and z directions, respectively.

[0162] Total strain energy U of cracked beam element m It can be written as:

[0163] U m =U n +U c

[0164] Among them, U n The strain energy of a crack-free beam element; U c This refers to the strain energy loss caused by the crack. n and U c These can be represented as follows:

[0165]

[0166]

[0167] Among them, A x and A c These represent the cross-sectional area of ​​the unit and the crack surface area, respectively; K Ir K IIr and K IIIr(r = 1, 2, ..., 6) represent the stress intensity factors for Type I, Type II, and Type III cracks, respectively. The expression for the stress intensity factor can be written as:

[0168]

[0169]

[0170]

[0171] F1(a), F2(a), F II (a) and F III (a) The stress intensity factor correction coefficient for the crack; a is the crack depth, h c The thickness of the section in the direction of the crack is denoted as .

[0172] Similarly, the nodal displacements of a cracked beam element can be divided into the displacements of a healthy beam element and the displacements of a cracked beam element. The nodal displacements can be expressed as follows:

[0173]

[0174] When the crack opens, the structural stiffness matrix of the cracked beam element can be written as:

[0175]

[0176] Among them, f nc and f c These are the compliance matrix of a healthy beam element and the additional compliance matrix caused by the introduction of cracks, respectively; T c e This is the transformation matrix.

[0177] S22, Simulated crack breathing method;

[0178] Furthermore, to simulate crack breathing, a crack breathing function needs to be established to determine the breathing state. Based on the vibration displacement of the crack surface at various times, the bending moment at the crack surface is calculated, and then the stress on the crack surface is calculated. The contact state of the crack surface is determined based on the stress value at the crack surface. The bending deformation caused by the crack leads to the crack surface circling around... e z e and o e y e Bending moment M of the shaft z and M y They can be represented as:

[0179]

[0180] in, and For the o e z e and oe y e Angular displacement of the shaft.

[0181] The structural stiffness matrix of the m-th healthy beam is replaced with the stiffness matrix of the cracked beam to obtain the structural stiffness matrix of the rotating cracked blade.

[0182] Figure 8 This is a schematic diagram of the crack cross-section. The crack depth ratio α is characterized by the proportion of the cracked region area to the total cross-sectional area, i.e., α = S1 / (S1 + S2). Secondly, the crack location ratio δ is introduced, assuming the cracked blade consists of n... all The segment beam is composed of the nth segment. crack If the segment beam is a cracked beam, then:

[0183]

[0184] S23. Simulate external loads on the blade and establish a dynamic model of a rotating airfoil section blade considering crack breathing effect;

[0185] To simulate the external loads experienced by compressor blades during operation, an aerodynamic force F is introduced here. p The aerodynamic force on each beam segment is The expression is as follows:

[0186]

[0187] Among them, A k Let A be the area of ​​the k-th unit that bears the aerodynamic force. p f r t and t represent the aerodynamic force amplitude, excitation frequency, and time acting on the rotating cracked blade, respectively; N v1 N v2 N v3 and N v4 These are four shape functions representing the bending displacement direction of the Timoshenko beam element.

[0188] Finally, the dynamic model of the rotating cracked blade is obtained:

[0189]

[0190] Where M is the mass matrix of the rotating cracked blade; K e G is the structural stiffness matrix of the rotating cracked blade; G is the Coriolis force matrix of the rotating cracked blade, and K is the structural stiffness matrix of the rotating cracked blade. c K is the centrifugal stiffening matrix of the rotating cracked blade. s K is the rotation softening matrix; acc The angular acceleration stiffness matrix of the rotating cracked blade is given by... Let angular acceleration be expressed, then it satisfies δ and F represent the nodal displacement and nodal force vector of the rotating cracked blade, respectively; F p The aerodynamic force on the rotating cracked blade.

[0191] S3. Using the dynamic model of the rotating airfoil blade considering the crack breathing effect verified in S2, the influence of different crack parameters and load parameters on the nonlinear response of the cracked blade is studied. The specific process is as follows:

[0192] The nonlinear vibration response of the rotating cracked blade considering the breathing effect was verified. Based on the finite element software ANSYS, a rotating breathing cracked airfoil section blade was constructed using Solid 185 solid elements, Conta174 contact elements, and Targe170 elements. The material parameters were set as follows: elastic modulus E = 125 GPa, Poisson's ratio μ = 0.3, and density ρ = 4370 kg / m³. 3 The crack parameters and load parameters for the dynamic model of the cracked blade are as follows: crack depth ratio α = 0.5, crack location ratio δ = 0.5, blade rotation speed n = 1000 r / min, and aerodynamic amplitude A. p =0.01MPa. During the solution process, each solution cycle is taken as 128 load steps, and the first-order resonant frequency f... c The calculation formula is:

[0193]

[0194] Among them, f open and f close These are the first-order bending natural frequencies of the open crack and the closed crack, respectively.

[0195] Figure 9 (a) and Figure 9 (b) Comparison of the time and frequency domains of the tip vibration response of the cracked blade using the beam model and solid model under superharmonic resonance and first-order resonance states, respectively. The figures clearly show that the response results of the beam model in this paper agree well with the results of the ANSYS solid model, proving the effectiveness of the model presented in this paper. Furthermore, the solution efficiency of the method presented in this paper is more than 20 times higher than that of ANSYS, significantly saving solution time.

[0196] In engineering, aero-engine blades exhibit a variety of crack types, and the loads they bear are also variable. It is necessary to study the influence of different parameters such as crack depth, crack location, and aerodynamic amplitude on the nonlinear vibration response characteristics of cracked blades. This invention only presents the influence of aerodynamic amplitude on the nonlinear vibration response characteristics of cracked blades. The simulation uses the Newmark numerical integration method to solve the vibration response signal of the finite element model and extracts the y-axis at the tip node of the cracked blade. bThe vibration characteristics of cracks are analyzed by directional vibration displacement.

[0197] With fixed parameters such as crack depth ratio, crack location ratio, and rotational speed, and crack depth ratio α = 0.5, crack location ratio δ = 0.5, and rotational speed n = 1000 r / min, the tip vibration response of a cracked airfoil under different aerodynamic amplitudes was studied. Aerodynamic amplitude A p The change increment ΔA is taken from 0.01 MPa to 0.1 MPa. p =0.01MPa.

[0198] Considering both superharmonic resonance and first-order resonance vibration states of the cracked blade, the resulting three-dimensional spectrum is as follows: Figure 10 As shown. By Figure 10 It can be seen that the amplitude of the nonlinear vibration response of the cracked airfoil gradually increases with the increase of the aerodynamic amplitude. The vibration amplitude of the blade in the first-order resonance state is significantly greater than that in the superharmonic resonance state, and in the superharmonic resonance state, the amplitude of the vibration response at the first harmonic frequency (f) is significantly greater. r and 2 times the frequency 2f r The amplitudes are quite similar. The amplitude variations of each harmonic of the cracked blade under the two vibration states are compared as follows: Figure 11 and 12 As shown in the figure, under both vibration states, as the aerodynamic amplitude increases linearly, the amplitudes of each harmonic also show a regular linear increase.

[0199] The sequence numbers of the above embodiments of the present invention are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments.

[0200] In the above embodiments of the present invention, the descriptions of each embodiment have different focuses. For parts not described in detail in a certain embodiment, please refer to the relevant descriptions of other embodiments.

[0201] In the several embodiments provided in this application, it should be understood that the disclosed technical content can be implemented in other ways. The device embodiments described above are merely illustrative; for example, the division of units can be a logical functional division, and in actual implementation, there may be other division methods. For instance, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the displayed or discussed mutual coupling, direct coupling, or communication connection may be through some interfaces; the indirect coupling or communication connection between units or modules may be electrical or other forms.

[0202] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0203] Furthermore, the functional units in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit.

[0204] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, read-only memory (ROM), random access memory (RAM), portable hard drives, magnetic disks, or optical disks.

[0205] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method of modeling the dynamics of a breathing-effect rotating-foil variable- cross-section crack, characterized by, Includes the following steps: S1. Considering the effects of centrifugal stiffening, rotational softening, and Coriolis force, a dynamic model of a real blade rotating micro-segment beam element with variable cross-section is established based on finite element theory and Timoshenko beam theory. S2. Based on the strain energy release rate and Castingliano's theorem, a cracked beam element is introduced into the dynamic model to establish a dynamic model of a rotating airfoil section blade that considers the crack breathing effect. The dynamic model of a rotating airfoil blade with cracks includes healthy beam elements and cracked beam elements. S2 specifically includes the following steps: S21. First, a dynamic model of a healthy rotating airfoil section blade is established based on beam element theory. Then, the healthy beam elements on the model are replaced with cracked beam elements to achieve the purpose of implanting cracks. S22. By measuring the vibration displacement of the crack surface at various times, calculate the bending moment at the crack surface, calculate the stress on the crack surface, and determine the contact state of the crack surface based on the stress value at the crack surface, thereby simulating crack breathing. S22 specifically includes: A crack breathing function is established to determine the breathing state. According to the vibration displacement of the crack surface at each time, the bending moment at the crack surface is calculated, and then the stress at the crack surface is calculated. According to the stress value at the crack surface, the contact state of the crack surface is determined. The bending deformation caused by the crack causes the crack surface to rotate around the axis o e z e And o e y e The bending moment of the shaft M z And M y are respectively represented as: wherein and is the rotational angular displacement about o e z e and o e y e the axis of rotation; The structural stiffness matrix of the healthy beam is replaced by the stiffness matrix of the cracked beam to obtain the structural stiffness matrix of the rotating cracked blade. m segment health beam structure stiffness matrix is replaced by the stiffness matrix of the cracked beam to obtain the structural stiffness matrix of the rotating cracked blade. The crack depth ratio is represented by the proportion of the crack area to the total cross-sectional area α That is α =S1 / (S1+S2); the crack position ratio is introduced δ Suppose that the crack blade is composed of n all The first n crack The crack beam is the crack beam. S23. The external load on the blade was simulated using aerodynamics, and a dynamic model of a rotating airfoil blade considering the crack breathing effect was established. S23 specifically includes: Simulates the external loads borne by compressor blades during operation, introducing aerodynamic forces. F p The aerodynamic force on each beam segment is F ep k The expression is as follows: in, A k For the first k The area of ​​each unit that bears aerodynamic forces; A p , f r and t These represent the aerodynamic force amplitude, excitation frequency, and time experienced by the rotating cracked blade, respectively. N v1 , N v2 , N v3 and N v4 These are four shape functions representing the bending displacement direction of the Timoshenko beam element; Finally, the dynamic model of the rotating cracked blade is obtained: in, M The mass matrix of the rotating cracked blade; K e is the structural stiffness matrix of the rotating cracked blade; G The Coriolis force matrix of the rotating cracked blade. K c is the centrifugal stiffening matrix of the rotating cracked blade. K s is the rotation softening matrix; K acc is the angular acceleration stiffness matrix of the rotating cracked blade, using Let angular acceleration be expressed, then it satisfies ; δ and F These represent the nodal displacement and nodal force vector of the rotating cracked blade, respectively. F p The aerodynamic force on the rotating cracked blade.

2. The method for dynamic modeling of a rotating airfoil with a variable cross-section cracked blade containing a breathing effect according to claim 1, characterized in that, S1 specifically includes the following steps: S11. Establish an airfoil blade model; Imagine the airfoil as a rotating Timoshenko cantilever beam, and divide the Timoshenko cantilever beam into... N Unit; With the first k Taking a rotating Timoshenko straight beam element as the research object, the differential equations of motion are established; among them... l k For the first k The length of a discrete finite element; β 0 represents the initial installation angle of the blade. γ k For the first k The pre-torsion angle of the Timoshenko beam element, the first k The installation angle of a Timoshenko beam element is defined as follows: β k = β 0+ γ k ; α 0 represents the blade's rotational angular displacement; S12. Determine the blade cross-section parameters; The parameters of each section are determined based on the geometric properties of the actual section; the first... k The cross-section of each element is first described in terms of coordinates. x k The closed curves that form the profile of this cross section are used to parametrically represent the function expressions of the blade tip and blade back through spline fitting; the expression of the cross section area, the coordinate transformation matrix of the local coordinate system and the element coordinate system, the expression of the centroid coordinates of the cross section, the moment of inertia of the cross section about the axis, the torsional stiffness of the airfoil section, and the additional torsional moment of inertia of the cross section caused by the pre-torsion angle are obtained. S13. Considering the influence of the pre-twist angle on the potential energy of the micro-segment straight beam element, an energy correction term is introduced. After correction, the total potential energy of the rotating variable cross-section micro-segment straight beam is obtained. U e The calculation expression is obtained; then, combined with the kinetic energy of the rotating variable cross-section micro-segment straight beam, the dynamic model of the rotating micro-segment beam element is obtained.

3. The method for dynamic modeling of a rotating airfoil with a variable cross-section cracked blade containing a breathing effect according to claim 2, characterized in that, In S11, a point on any cross section of the blade P The displacements in the three directions in the global coordinate system can be expressed as: in, θ x , θ y and θ z For describing the cross-section of an arbitrary element x e axis, y e shaft and z e The second kind of Euler angle is introduced by the rotational angular displacement of the shaft; x , y , z for P The coordinates of the point in the local coordinate system; R k Origin of unit coordinates o e In rotating coordinate system OX r Y r Z r Below x Coordinates of direction; u , v and w These are the centroids of the element sections along the element coordinate system. o e x e axis, o e y e shaft and o e z e Displacement of the axis.

4. The method for dynamic modeling of a rotating airfoil with a variable cross-section cracked blade containing a breathing effect according to claim 2, characterized in that, In S12, the torsional stiffness of the airfoil section J s The empirical expression is: in, A For cross-sectional area, I max and I min These are the maximum and minimum moments of inertia of the cross section, respectively; due to the pre-twist angle γ ( L Additional torsional moment of inertia of the cross section caused by ) J a The formula for calculation is: 。 5. The method for dynamic modeling of a rotating airfoil with a variable cross-section cracked blade containing a breathing effect according to claim 2, characterized in that, In S13, after correction, the total potential energy of the rotating variable cross-section micro-segment straight beam is obtained. U e The calculation expression is: in, U bte and U s The first k The bending-pendulum-axis-torsion coupling potential energy and shear potential energy of a micro-segment straight beam with variable cross-section by rotation. The rate of change of the twist angle, ; F b The dynamic model of the actual blade micro-segment beam element is obtained by representing the basic excitation load applied along the blade axis and simplifying the process: in, M e The mass matrix of a micro-segment beam element in the element coordinate system; K ee is the structural stiffness matrix of the Timoshenko beam element; G e Here is the Coriolis force matrix for the beam element. K ec is the centrifugal stiffening matrix of the beam element. K es is the rotation softening matrix; K eacc is the element angular acceleration stiffness matrix, used Let angular acceleration be expressed, then it satisfies ; δ e and F e These represent the element node displacement and element node force vector, respectively.

6. The method for dynamic modeling of a rotating airfoil with a variable cross-section cracked blade containing a breathing effect according to claim 1, characterized in that, S21 specifically includes: Total strain energy of cracked beam element U m for: in, U n The strain energy of a crack-free beam element; U c This refers to the loss of strain energy caused by cracks; U n and U c These can be represented as follows: in, A x and A c These are the cross-sectional area of ​​the unit and the crack surface area, respectively. K r , K r and K r These are the stress intensity factors for Type I, Type II, and Type III cracks, respectively. r =1, 2, ..., 6; The expression for the stress intensity factor is: in, F 1 ( a ), F 2 ( a ), F II ( a )and F III ( a ) is the correction factor for the stress intensity factor of the crack; a The crack depth. h c The cross-sectional thickness is along the crack direction; The nodal displacements of the cracked beam element are divided into the displacements of the healthy beam element and the displacements of the cracked beam element. The nodal displacements are as follows: When the crack opens, the structural stiffness matrix of the cracked beam element is: in, f nc and f c These are the compliance matrix of a healthy beam element and the additional compliance matrix caused by the introduction of cracks, respectively. T c e This is the transformation matrix.

Images