A method for dynamic analysis of a rotor system with blade rub

Abstract

The present application belongs to the technical field of rotor dynamics, and discloses a rotor system dynamics analysis method containing blade rub. The traversal cycle loading mode is used to simulate the rotor unbalance effect, the rotation of the blade is simulated, the vibration displacement of the blade tip node and the vibration displacement of the blade root node are calculated after the centrifugal deformation of the blade is calculated separately, the difference is calculated to obtain the blade deformation, the blade deformation coordinates are converted, the position vector of the disc and the position vector of the rotating shaft are superposed, the blade tip intrusion amount after coating wear is calculated, and the normal rub force vector of the multi-blade is constructed. Then, the beam element, the shell element and the solid element are used to construct the rotor system dynamics model containing blade rub under the motor flight condition, wherein the nonlinear force model is used to calculate the bearing to consider the influence of the nonlinearity on the rub behavior; the dynamics calculation is carried out through the reduction of the degree of freedom and the interface separation, and finally the rub dynamics response of the system is calculated through the numerical simulation analysis, and the influence of the motor analysis on the rub characteristics is analyzed.

Description

Technical Field

[0001] This invention relates to the field of rotor dynamics technology, and in particular to a method for dynamic analysis of rotor systems involving blade rubbing. Background Technology

[0002] As a core component of typical rotating machinery such as gas turbines and aero engines, the safety and stability of dual rotors are crucial. Because they consist of two rotors and their supports often utilize double-half inner ring multi-point contact ball bearings, their dynamic characteristics under internal and external excitations are more complex. However, in modeling dual rotor systems, deep groove ball bearings are commonly used, neglecting the multi-point contact characteristics of the double-half inner ring bearings. Furthermore, research on rubbing failures in dual rotor systems largely focuses on system modeling in an inertial frame and the nonlinear vibration characteristics induced by rubbing, neglecting the influence of non-inertial loads on rubbing failures. In particular, research considering coating wear in a non-inertial frame has not been published. Therefore, the vibration characteristics under the combined action of maneuvering flight loads and rubbing loads lack sufficient research, the coupling mechanism between the two is unclear, and related phenomena have not been thoroughly revealed. Moreover, for multi-blade rubbing, no published literature has proposed a normal rubbing force model considering the actual unbalance phase and blade phase relationship to calculate the load.

[0003] Blade rubbing has attracted considerable attention due to the higher linear velocity at the blade tip, the large rubbing energy, and the significant impact on rotor dynamics. Treating the blade-casing rubbing contact as rotor-stator contact is a common approach, including the penalty function method, the Lagrange multiplier method, and the augmented Lagrange multiplier method. However, these contact methods introduce iterative complexity, leading to extremely low computational efficiency. Therefore, methods based on rubbing force have garnered significant attention. Currently, the paper "Ma Hui, Guo Xumin, Zeng Jin, et al. A method for determining the blade-casing rubbing relationship: 201910878227.9 [P]. 2023-03-24" proposes a method for determining the blade-casing rubbing relationship. This method, based on a semi-analytical approach, derives the differential equations of motion between a single blade and the stator, calculates the blade tip intrusion, and then determines the rubbing force of a single blade and its induced dynamic response. The paper "Ma Hui, Guan Hong, Guo Xumin, et al. A method for calculating the normal impact force of a blade-coated casing considering thermal-wear effects: 202510938064.4 [P]. 2025-10-28" further constructs a single-blade-coated casing impact force model considering frictional heat and wear effects. This model integrates the composite stiffness characteristics of the casing coating and the supporting structure, thus enabling the impact force calculation to cover the response characteristics under different structural stiffness conditions. The paper "Yu Pingchao, Tao Xuanjun, Zhang Jia'e, et al. A simulation method for impact dynamics of complex swept blades in engines: 202510476654.X [P]. 2025-07-29" focuses on the true geometric shape of the blade and proposes a simulation method for impact dynamics of complex swept blades in engines. This method considers the true relative position and three-dimensional motion relationship between the blade and the casing, establishes an expression for the intrusion amount of the blade tip profile based on a multi-coordinate system, and aims to achieve accurate calculation of the true single-blade impact force. "Shi Tielin, Zhang Zhisong, Hu Cheng. A Simulation Method for Blade-Casing Friction Considering Structural Coupling of Aero-engines: 202310387762.0 [P]. 2023-07-07" proposes a simulation method for blade-casing friction considering structural coupling of aero-engines, focusing on multi-blade friction faults. However, it has not yet effectively considered the influence of unbalance phase and blade phase on the friction process. In summary, relatively rich mechanical models have been proposed for single-blade friction, while research on multi-blade friction remains insufficient.

[0004] Furthermore, most studies on the overall mechanical models of structures that consider rubbing effects focus on the rubbing dynamic response characteristics of dual-rotor systems under unbalanced loads while neglecting blade flexibility. A few studies also focus on the rubbing dynamic response of dual-rotor systems under maneuvering flight conditions, again neglecting blade flexibility. Therefore, to fill the gap in research on overall mechanical models of structures that consider rubbing effects, it is urgent to consider maneuvering flight and blade flexibility, and to construct a dynamic model of a complex rotating-stator system involving multi-blade rubbing.

[0005] The main drawbacks of existing blade rubbing models and rotor system dynamics modeling methods involving rubbing are summarized as follows: (1) Several local mechanical models have been proposed to describe the blade-casing collision process. These models are based on their own simplifications and assumptions, have different applicable conditions, and each has its own advantages and disadvantages. For example, although the contact dynamics-based method can analyze the high-speed collision process well, it needs to deal with the vibration deformation of the blade-casing, complex contact judgment, and numerical iteration in the calculation, resulting in low computational efficiency and difficulty in convergence. Therefore, it is more difficult to apply it to the collision dynamics analysis of the whole machine. In contrast, the blade-casing collision model based on the quasi-static deformation assumption can explicitly express the relationship between the collision force and the structural and operating parameters, which is convenient to embed in the rotor or even the whole machine collision simulation analysis. However, this model has not effectively considered the influence of blade phase when dealing with multi-blade collisions. In addition, the speed fluctuation caused by collisions has not been effectively considered.

[0006] (2) Most studies on rubbing failures in dual-rotor systems focus on system modeling in the inertial frame and the nonlinear vibration characteristics induced by rubbing, while neglecting the influence of maneuvering flight loads and nonlinear bearing support loads on rubbing failures. In addition, for blade-casing rubbing, most studies simplify the bladed disk structure to a disk and ignore the flexibility of the blades, making it impossible to accurately analyze the vibration behavior of the blades; although some scholars retain the real bladed disk structure, they do not accurately consider the influence of rotor vibration on blade vibration during the rubbing process, resulting in distortion of the tip clearance and rubbing force calculations in the analysis of multi-blade rubbing.

[0007] To address the issues of lacking a high-fidelity multi-blade rubbing model, unclear rotor-stator vibration coupling mechanism, and ambiguous rubbing phenomena in the aforementioned studies, a multi-blade rubbing model is proposed based on vector synthesis and a "pseudo-rotation" strategy to simulate blade rotation. Considering the influence of the vibration displacement of the inner and outer rings of the bearing on the ball contact characteristics, a quasi-static model of a double-inner-ring multi-point contact ball bearing is constructed. Furthermore, based on the dynamic modeling technology of a dual-rotor system under maneuvering flight, a dynamic model of a double-inner-ring bearing-dual-rotor-disc-coated-casing system with blade-coated-casing rubbing fault in a non-inertial frame is established. The influence of different inertial loads on the rubbing-induced rotor-stator vibration response and wear under maneuvering flight conditions is analyzed, revealing the influence of rotor-stator structural parameters, speed fluctuations, coating parameters, and maneuvering flight on the rubbing response.

[0008] To quantitatively analyze the impact of rubbing on the dynamic characteristics of the engine's rotor-stator system, it is first necessary to establish a mechanical model that can fully reflect the rubbing process. This model mainly consists of two parts: a local rubbing mechanical model describing the rubbing contact process, and a structural overall mechanical model that takes into account the rubbing effect. Summary of the Invention

[0009] To address the aforementioned problems, this invention proposes a dynamic analysis method for rotor systems involving blade rubbing. This method considers the phase relationships of the actual blades and incorporates multi-blade rubbing force vectors into the dynamic model of the rotor system. First, this invention uses ergonomic cyclic loading to simulate rotor imbalance effects and blade rotation. After calculating the centrifugal deformation of the blades individually, the blade deformation is calculated by subtracting the vibration displacement of the blade tip node and the vibration of the blade root node. The blade deformation coordinates are then transformed and superimposed with the position vectors from the disk center to the blade root and the rotation axis to calculate the blade tip position vector and the blade tip intrusion amount considering coating wear, thus constructing a multi-blade normal rubbing force vector. Next, beam elements, shell elements, and solid elements are used to construct a dynamic model of the rotor system involving blade rubbing under maneuvering flight conditions. Bearings are calculated using nonlinear force models to account for their nonlinear effects on rubbing behavior. Dynamic calculations are performed by reducing degrees of freedom and interface separation. Finally, numerical simulation analysis is used to calculate the system's rubbing dynamic response and analyze the impact of maneuvering analysis on rubbing characteristics.

[0010] The technical solution of this invention is as follows: A dynamic analysis method for a rotor system with blade rubbing, which calculates the dynamic response based on the excitation load, considers the influence of frictional torque caused by rubbing on the rotational speed, obtains the instantaneous rotational speed, and calculates the blade azimuth angle; iteratively applies the excitation load to simulate the unbalanced effect of the rotor system with blade rubbing, calculates the centrifugal deformation of the blade, obtains the vibration displacement of the blade tip node and the vibration displacement of the blade root node, and uses the difference between the two as the blade deformation and calculates the blade position vector in the rotating coordinate system. The blade position vector is then projected onto the global coordinate system after coordinate transformation; considering coating wear, the blade position vector projected onto the global coordinate system is superimposed with the position vector from the disk center to the blade root and the position vector of the rotating shaft to calculate the blade tip position vector and the blade tip intrusion amount after considering coating wear, constructing a multi-blade rubbing force vector; a dynamic model of the rotor system with blade rubbing under maneuvering flight conditions is constructed using beam elements, shell elements, and solid elements, and dynamic calculations are performed by reducing the degrees of freedom and separating the interfaces; the rubbing dynamic response of the rotor system with blade rubbing is calculated through numerical simulation analysis.

[0011] The frictional torque is calculated based on the frictional force, and the angular acceleration is calculated based on the moment of inertia of the rotor system containing blade rubbing. The instantaneous rotational speed is then calculated based on the angular acceleration. exist Under the influence of the action, the formulas for calculating angular acceleration and instantaneous rotational speed are: In the formula, Angular acceleration, For the resultant torque, For rotational inertia, The radial distance from the point of contact to the axis of rotation. The rotational speed before the collision occurred. Instantaneous rotational speed t rub For the duration of the rubbing, This is the speed feedback coefficient, with units of N·m / (rad / s). For the target speed, The total frictional force at all blade tip nodes. For the first i Friction at the leaf tip node, The number of nodes at the leaf tip. For load torque, This refers to the motor drive torque. This represents the frictional torque.

[0012] blade azimuth angle α b The expression is: In the formula, This indicates the initial blade azimuth angle.

[0013] The simulated unbalance effect of a rotor system with blade rubbing is specifically achieved by cyclically applying the excitation load, calculating the unbalance force of the rotor system with blade rubbing based on the calculated instantaneous rotational speed, and then performing force decomposition calculations. x Unbalanced forces F unx and y Unbalanced forces F uny ; The expression for the unbalanced force of a rotor system including blade rubbing is: In the formula, m un It is unbalanced mass. e un It is the eccentricity of unbalanced mass.

[0014] The blade deformation includes the centrifugal deformation vector of the disk and the blade as a whole. and the centrifugal deformation vector of the leaf tip relative to the leaf root for: In the formula, and These are the centrifugal deformation vectors of the leaf tip and the leaf root in the local coordinate system of the leaf root, respectively. This is the inverse of the stiffness matrix of the rotor system containing blade rubbing. The centrifugal force vector of the rotor system containing blade rubbing; Under the action of frictional force, the blade moves in the local coordinate system at the leaf root. o b x b y b Dynamic deformation vector The formula for calculation is: In the formula, ( x b , y b ) represents the deformation of the leaf tip in the local coordinate system of the leaf root, ( x d , y d () represents the deformation of the leaf root in the local coordinate system of the leaf root; Calculate the position vector of the blade in the rotating coordinate system, and project the position vector of the blade in the rotating coordinate system onto the global coordinate system after coordinate transformation based on the blade azimuth angle; Leaf tip node A The position vector in the local coordinate system at the leaf root is: In the formula, Let be the position vector of the leaf tip node in the local coordinate system at the leaf root. e t This refers to the amount of tip elongation caused by thermal effects; The blade length; In a rotating coordinate system, it is denoted as Its expression is: In the formula, The angle between the local coordinate system and the rotated coordinate system at the leaf root; After transforming from the rotating coordinate system to the global coordinate system using a coordinate transformation matrix, the blade position vector in the global coordinate system... r b3 The expression is: In a rotating coordinate system, the position vector from the center of the disk to the leaf root. The expression is: In the formula, R d The radius of the disk; x di , y di() represents the deformation of the disk node in the rotating coordinate system; In the global coordinate system, the position vector from the center of the disk to the leaf root. r b2 The expression is: After considering parallel misalignment, the leaf tip node in the global coordinate system A Position vector for: , In the formula, e x and e y for x To static parallel misalignment quantity and y Towards static parallel misalignment; and For the shaft at x The vibration displacement and rotation axis in the direction y Vibration displacement in the direction, This is the position vector of the rotation axis in the global coordinate system.

[0015] Unit direction vector of frictional force and the unit direction vector of frictional force for: In the formula, , and for B The position vector of a point C The position vector of the point A The position vector of the point; the casing is uniformly and equally divided along the circumference using the finite element method. N c 1 node B Point and C The point is the leaf tip node. A Two nodes that are circumferentially adjacent, and B Point and x The included angle of the axis is greater than C Point and x Angle between axes; blade tip node A exist BC The foot of the perpendicular on the line is A 'point; , and The expression is: In the formula, ( , ), ( , )and( , )for B The coordinates of the point C The coordinates of the point and A 'The coordinates of the point; α B and α C Coated casing B The phase angle corresponding to the point and the coating on the casing C The phase angle corresponding to the point; u B and v B for B Normal displacement of the point and B Tangential displacement of the point; u C and v C for C Normal displacement of the point and C Tangential displacement of the point; The radius of the coated casing; The tip intrusion amount is calculated based on the corresponding displacement relationship, and the expression is: In the formula, For the tip of the blade without considering coating wear i The amount of intrusion into each node; All nodes at the blade tip are matched to discrete points in the coated casing based on the principle of closest phase. t The wear amount at discrete points of the coating is calculated based on the contact rubbing state at a given time. t At time +1, the wear amount at the corresponding discrete point of the coating casing is read through the phase of the blade tip node to update the blade tip intrusion amount. The blade tip intrusion amount after considering coating wear is: In the formula, e w This refers to the amount of wear on the coating. To account for the amount of blade tip intrusion after coating wear.

[0016] The process of constructing the multi-blade friction force vector is as follows: Leaf tip i Friction vector at each node and friction force vector for: In the formula, k eq For the equivalent stiffness of the coated casing, For the speed correlation coefficient, μ The coefficient of friction; After weighted averaging of the friction force vector and the rubbing force vector, the first tip of the leaf... i Averaged friction force vector at each node and leaf tip i The friction force vector after averaging at each node for: , When leaf tip node A When subjected to frictional force, the coated casing experiences a reaction force. According to static equilibrium and torque equilibrium, B Point and C The vectors of the contact force and the friction force acting on the point are: , In the formula, and for B The sum of the frictional force vectors at the points C The vector of the frictional force at the point. and for B The frictional force vector at the point and C The frictional force vector at a point, for A 'Point and C Distance between points for B Point and C The distance between points is used to calculate the tip intrusion and friction force of any blade using the same calculation method to construct a multi-blade friction force vector.

[0017] The construction of the rotor system dynamics model with blade rubbing under the maneuvering flight conditions is as follows: Considering the influence of coating wear on the blade tip clearance, the rubbing force vector and friction force vector are averaged and applied to all nodes at the blade tip and the corresponding nodes of the coating casing; the rotor system dynamics model with blade rubbing under the maneuvering flight conditions is constructed using beam elements, shell elements and solid elements. In the formula, M dr , C dr , G dr and K drThe mass matrix, damping matrix, gyroscope matrix, and stiffness matrix of the rotor system containing blade rubbing are given. The damping matrix of the rotor system containing blade rubbing is calculated using Rayleigh damping. F dr This is the updated excitation load vector for a rotor system including blade rubbing, which incorporates the maneuvering flight inertial load vector. F man Misaligned load vector F m Multi-blade rubbing force vector F rub Centrifugal force vector F ce Unbalanced force vector of the rotor F u and the nonlinear contact force vector of the bearing F b ; F rub Include , , , , , ; F u Include F unx and F uny , F b Includes the contact load of double-ring inner ring angular contact ball bearings and the contact load of deep groove ball bearings; , and Let be the updated acceleration, velocity, and displacement vector of the rotor system with blade rubbing. For the system dynamics model of a rotor system including blade rubbing under maneuvering flight conditions, Newmark- β The dynamic response is calculated by first performing LU decomposition on the stiffness matrix, and then calculating the contact load of the double-half inner ring angular contact ball bearing at each time step using the Newton-Raphson algorithm. The interface separates the stiffness matrix and damping matrix into constant and fluctuating parts: In the formula, and The constant part and the wave part of the damping matrix are respectively represented by the damping matrix. and This refers to the constant part and the fluctuating part of the stiffness matrix; for the constant part of the stiffness matrix... The LU decomposition technique is used to decompose it into an upper triangular part and a lower triangular part; The dynamic model of a rotor system with blade rubbing under maneuvering flight conditions after interface separation is as follows: The dynamic model of a rotor system with blade rubbing under reduced degrees of freedom maneuvering flight conditions is as follows: In the formula, , and The mass matrix, damping matrix, and stiffness matrix of the rotor system with blade rubbing and reduced degrees of freedom are given. , and Let be the acceleration vector, velocity vector, and displacement vector of the rotor system with blade rubbing, with reduced degrees of freedom. The external excitation load vector for a rotor system with blade rubbing and reduced degrees of freedom.

[0018] Compared with existing dynamic analysis methods for rotor systems involving blade rubbing, the method proposed in this invention has the following advantages: (1) This dynamic model can calculate the actual blade azimuth angle and effectively consider the influence of rotor unbalance phase and blade phase angle on blade rubbing behavior. (2) This dynamic model can comprehensively consider complex load environments, including but not limited to rotor unbalance excitation, additional load excitation induced by maneuvering flight, nonlinear support excitation of bearings, rubbing force excitation, etc., and the dynamic model has high calculation accuracy. (3) This dynamic model adopts reduced degrees of freedom and interface separation, resulting in high calculation efficiency. Therefore, the dynamic analysis method for rotor systems involving blade rubbing constructed in this invention can efficiently and accurately analyze the vibration behavior of the system and has strong engineering application value. Attached Figure Description

[0019] Figure 1 The diagram shows the non-point contact friction of multiple blades; (a) shows no contact friction, (b) shows non-point contact friction of multiple blades, and (c) shows the deformation and stress analysis.

[0020] Figure 2 The diagram illustrates the effects of unbalanced phase; (a) represents unbalanced phase 0°; (b) represents unbalanced phase 90°; and (c) represents unbalanced phase 180°.

[0021] Figure 3Schematic diagram for blade tip clearance calculation; (a) blade-casing; (b) blade tip contact; (c) friction force; Figure 4 A flowchart of the dynamic analysis method for a rotor system involving blade rubbing; Figure 5 A finite element model of a rotor system containing blade rubbing. Figure 6 This is a diagram illustrating a dive-climb maneuver. Figure 7 The displacement time-domain waveforms of the bladed disk nodes under the dive-climb maneuver are shown; (a) is the displacement time-domain waveform in the x-direction; (b) is the displacement time-domain waveform in the y-direction. Figure 8 (a) is the time-domain waveform of the normal friction force under the dive-climb maneuver; (b) is an enlarged view of the time-domain waveform of the friction force. Detailed Implementation

[0022] The frictional torque caused by rotor-stator contact can lead to fluctuations in rotor speed, and in extreme cases, shaft jamming. Specifically, during continuous friction between the rotor and stator, the contact surfaces generate high temperatures, causing localized melting and rotor jamming. In this situation, the rotor is subjected to enormous torsional torque, leading to torsional fracture; this phenomenon is known as jammed torsional fracture failure. Taking the Jeffcott rotor system as an example, assuming the mass of the contact disk is... m d , radius is R d The motor drive torque is T d The average load torque is T l The average speed of the rotor is ω av The instantaneous speed of the rotor is ω ins The formulas for calculating angular acceleration and instantaneous angular velocity are: In the formula, Angular acceleration, For the resultant torque, For rotational inertia, The radial distance from the point of contact to the axis of rotation. The rotational speed before the collision occurred. Instantaneous rotational speed t rub For the duration of the rubbing, This is the speed feedback coefficient, with units of N·m / (rad / s). For the target speed, The total frictional force at all blade tip nodes. For the tip friction force, The number of nodes at the leaf tip. For load torque, This refers to the motor drive torque. This represents the frictional torque.

[0023] blade azimuth angle α b The expression is updated to: In the formula, This represents the initial azimuth angle.

[0024] The expression for the unbalanced force of a rotor system including blade rubbing is: In the formula, m un It is unbalanced mass. e un It is the eccentricity of unbalanced mass.

[0025] Rotor systems involving blade rubbing are multi-blade-casing systems; multi-blade-casing systems in a stationary state, such as... Figure 1 As shown, each blade tip has an initial radial clearance with the coated casing, indicated by a solid red line; the casing is indicated by a solid green line, the coating by a solid yellow line, and the blade by a solid black line. Under external excitations such as rotor imbalance, misalignment, and maneuvering flight, multi-blade non-fixed-point rubbing may occur between the blade and the coated casing. Specifically, multiple blades may make contact with the coated casing at different circumferential positions (see...). Figure 1 (b) The rubbing area is indicated by a red dashed line. Once a rubbing failure occurs, both normal and tangential forces will act on the blade (see...). Figure 1 (c) The coated casing affects system vibration and causes speed fluctuations, while the disk deforms. Because the blades are established in the local coordinate system at the blade root, while the shaft and disk are established in the global coordinate system, the blade deformation is calculated and processed by transforming the coordinate system to the global coordinate system. Different unbalance phases affect the blade's vibration response, such as... Figure 2 As shown. Figure 2In (a), the green arrow represents the unbalanced force, the green dashed line represents the blade tip trajectory, the red dashed line represents the rotor trajectory, and the black arrow represents the rotor vortex direction. To consider the influence of the unbalanced phase and the blade phase, this invention uses a cyclic loading method to simulate the rotor unbalanced effect and the blade rotation. First, the centrifugal deformation of the blade is calculated separately. After calculating the dynamic response, the blade deformation is calculated by subtracting the vibration displacements of the blade tip and root nodes. The blade deformation is then transformed by coordinates and superimposed with the position vectors of the disk and the rotation axis to indirectly simulate the rotation of the blade model. The centrifugal deformation vectors of the disk and the blade as a whole are... r st and the centrifugal deformation vector of the leaf tip relative to the leaf root. r st-tr for: In the formula, and These are the centrifugal deformation vectors of the leaf tip and the leaf root in the local coordinate system of the leaf root, respectively. This is the inverse of the stiffness matrix of the rotor system containing blade rubbing. The centrifugal force vector of the rotor system containing blade rubbing; Under the action of frictional force, the blade moves in the local coordinate system at the leaf root. o b x b y b Dynamic deformation vector r bt The formula for calculation is: In the formula, ( x b , y b ) represents the deformation of the leaf tip in the local coordinate system of the leaf root, ( x d , y d () represents the deformation of the leaf root in the local coordinate system of the leaf root; Leaf tip node A The position vector in the local coordinate system at the leaf root is: In the formula, Let be the position vector of the leaf tip node in the local coordinate system at the leaf root. e t This refers to the amount of tip elongation caused by thermal effects; The blade length; Due to the installation angle of the blade βblade The influence of the leaf root local coordinate system o b x b y b and rotating coordinate system o r x r y r Angle exists β blade To achieve vector addition in the global coordinate system, it is necessary to transform and rotate the vectors to a rotating coordinate system. The transformation matrix is... T 1g . r stt In a rotating coordinate system, it is denoted as Its expression is: After transforming from the rotating coordinate system to the global coordinate system using a coordinate transformation matrix, the blade position vector in the global coordinate system... r b3 The expression is: In a rotating coordinate system, the position vector from the center of the disk to the leaf root. The expression is: In the formula, R d The radius of the disk; x di , y di () represents the deformation of the disk node in the rotating coordinate system; In the global coordinate system, the position vector from the center of the disk to the leaf root. r b2 The expression is: After considering parallel misalignment, the leaf tip node in the global coordinate system A Position vector r A for: , In the formula, e x and e y for x To static parallel misalignment quantity and y Towards static parallel misalignment; and For the shaft at x The vibration displacement and rotation axis in the direction y Vibration displacement in the direction, This is the position vector of the rotation axis in the global coordinate system.

[0026] Unit direction vector of frictional force and the unit direction vector of frictional force for: In the formula, r B , r C and r A' for B point, C Point and A The position vector of a point is expressed as follows: In the formula, ( x B , y B ), ( x C , y C )and( x A' , y A' )for B point, C Point and A 'The coordinates of the point, such as Figure 3 As shown. α B 、α C and α D Coated casing B, C and D The phase angles corresponding to the three points. u B and v B for B The normal and tangential displacements of a point; u C and v C for C The normal and tangential displacements of a point; The radius of the coated casing; The tip intrusion can be calculated based on the corresponding displacement relationship, and its expression is: In the formula, For the tip of the blade without considering coating wear i The amount of intrusion into each node; To account for the impact of coating wear on blade tip clearance, the coating casing is uniformly and equally dispersed circumferentially. N c A point, such as Figure 3 In (b), the solid black line represents the surface of the casing, and the dashed red line represents the blade tip; Figure 3 In (c), the casing surface is represented by a red dashed line. All nodes at the blade tip are matched to discrete points on the casing based on the closest phase principle. t The wear amount at discrete points of the coating is calculated based on the contact rubbing state at a given time. t At time +1, the wear amount of the corresponding discrete point of the casing is read through the phase of the blade tip node to update the blade tip intrusion amount. Therefore, the blade tip intrusion amount is updated to: In the formula, e w This represents the amount of wear on the coating.

[0027] Leaf tip i The tangential friction force vector of each node and normal friction force vector for: In the formula, L b For the blade length, k eq For the equivalent stiffness of the coated casing, For the speed correlation coefficient, μ is the coefficient of friction.

[0028] After weighted averaging of the friction force vector and the rubbing force vector, the first tip of the leaf... i Averaged friction force vector at each node and leaf tip i The friction force vector after averaging at each node for: , When the leaf tip A When a point is subjected to a frictional force, the coated casing will experience a reaction force. According to static equilibrium and torque equilibrium, B Point and C The normal and tangential forces acting on the point are: , In the formula, and for B The sum of the frictional force vectors at the points C The vector of the frictional force at the point. and for B The frictional force vector at the point and C The frictional force vector at a point, for A 'Point and C Distance between points for B Point and C The distance between points is used to calculate the tip intrusion and friction force of any blade using the same calculation method to construct a multi-blade friction force vector.

[0029] The system dynamics model of a rotor system with blade rubbing under maneuvering flight conditions without reduced degrees of freedom is as follows: In the formula, M dr , C dr , G dr and K dr The mass matrix, damping matrix, gyroscope matrix, and stiffness matrix of the rotor system containing blade rubbing are given. The damping matrix is ​​calculated using Rayleigh damping. F dr This is the updated excitation load vector for a rotor system including blade rubbing, which incorporates the maneuvering flight inertial load vector. F man Misaligned load vector F m Multi-blade rubbing force vector F rub Centrifugal force vector F ce Unbalanced force vector of the rotor F u and the nonlinear contact force vector of the bearing F b ; F rub Include , , , , , ; F u IncludeF unx and F uny , F b Includes the contact load of double-ring inner ring angular contact ball bearings and the contact load of deep groove ball bearings; , and Let be the updated acceleration, velocity, and displacement vector of the rotor system with blade rubbing. For the system dynamics model of a rotor system with blade rubbing under maneuvering flight conditions without reduced degrees of freedom, Newmark- β The dynamic response is calculated by first performing LU decomposition on the stiffness matrix, and then calculating the contact load of the double-half inner ring angular contact ball bearing at each time step using the Newton-Raphson algorithm. The system stiffness matrix and damping matrix are separated into constant and fluctuating components using nonlinear interface separation technology: In the formula, and The constant and fluctuating parts of the damping matrix are represented by the two components. and This represents the constant and fluctuating portions of the stiffness matrix. For a constant stiffness matrix... The equation is decomposed into an upper triangular part and a lower triangular part using the LU decomposition technique, which greatly improves the computational efficiency of the linear part of the equation.

[0030] The dynamic model of a rotor system with blade rubbing under maneuvering flight conditions after interface separation is written as follows: The dynamic model of a rotor system with blade rubbing under reduced degrees of freedom maneuvering flight conditions is as follows: In the formula, , and The mass matrix, damping matrix, and stiffness matrix of the rotor system with blade rubbing and reduced degrees of freedom are given. , and Let be the acceleration vector, velocity vector, and displacement vector of the rotor system with blade rubbing, with reduced degrees of freedom. The external excitation load vector for a rotor system with blade rubbing and reduced degrees of freedom.

[0031] Taking a dual-rotor-disc-casing system with double-half-inner ring bearings as an example, the calculation process for a rotor system involving blade rubbing is shown below. Figure 4 The dual-rotor-disc-casing system with double-half-inner ring bearings features inner and outer shafts driven by motors, with vibration coupling achieved through intermediate bearings. The inner shaft is supported on both sides by deep groove ball bearings, while the front support of the outer shaft is supported by double-half-inner ring bearings, simulating the support configuration of the high-pressure rotor's front support. In the modeling process of this invention, the inner shaft, outer shaft, and inner casing are modeled using beam elements, the inner and outer bladed disks using shell elements, and the rubbing coated outer casing is modeled using solid elements to account for local vibrations of the coated casing, effectively simulating the rubbing contact behavior between the blades and the coated casing.

[0032] Based on beam elements, shell elements, and solid elements, the dynamic model of the dual-rotor-disc-coated casing is as follows: Figure 5 As shown, the black arrow represents the excitation load; the red arrow points to the left view of the dynamic model. For this model, the mesh is refined axially in the potential rubbing area of ​​the coated casing to match the axial positions of the blade tip nodes and the coated casing nodes. As can be seen from the left view, due to the different coating thicknesses, the initial clearance at the blade tip varies, thus potentially leading to varying degrees of circumferential rubbing.

[0033] Dive-climb maneuvers result in additional parameter and load excitations, altering the rotor system's vibration characteristics and blade tip amplitude, thereby inducing blade-coated casing rubbing. A schematic diagram of a dive-climb maneuver is shown below. Figure 6 As shown, it is divided into 7 stages, with the forward direction from right to left. It can be briefly summarized as level flight, entering a dive, continuing to dive, diving at a constant speed and pulling up, continuing to climb, exiting the climb and level flight. This invention focuses on the vibration state of the rotor system in the first four stages.

[0034] The time-domain waveform of the displacement of the bladed disk node during maneuvering flight is as follows: Figure 7 As shown in the figure, a dive maneuver significantly alters the beam element nodes at the bladed disk location. x Vibrational displacement. During the dive phase, x The displacement is significantly biased towards the positive half-axis and the amplitude is increasing; during the continued dive phase... x The displacement is significantly biased towards the negative half-axis and its amplitude decreases; during the uniform dive and pull-up phase, the displacement is further biased towards the negative half-axis. This change in displacement directly affects the tip intrusion, thereby altering the tip friction force and the system's motion state. For y Vibration displacement in the direction, such as Figure 7 As shown in (b), the opposite pattern of change can be observed. Furthermore, x The displacement amplitude in the direction is significantly greater than y Displacement amplitude.

[0035] Maneuvering flight can alter the vibration patterns of the rotor shaft system due to additional excitation, thereby affecting the normal friction force between the blades and the casing, such as... Figure 8 As shown, there are 16 blades in total, each represented by a different colored dashed line. Figure 8 (a) illustrates the evolution of the normal friction force of all blades under maneuvering flight conditions, while Figure 8 (b) shows a magnified view of the blade normal friction force from 1.5s to 1.53s (corresponding to the dive entry phase). As can be seen from the figure, the normal friction force between the blades and the coated casing increases significantly during the dive entry phase due to rotor vibration displacement changes, and remains the highest for blade #15, similar to the level flight phase. However, once the dive continues and then accelerates into the pull-up phase, the maximum friction force is slightly greater than that of the level flight phase, and the blades with the highest friction force change; at this point, blades #11 and #12 have the highest friction force. This is because the maneuvering flight simultaneously alters... x Towards vibration displacement and y This is caused by vibrational displacement. Furthermore, it is due to... Figure 8 (b) It can be seen that during the dive phase, although the maximum friction force increases, the number of blades involved in friction decreases, and blades 9 to 11 no longer experience friction failures. These phenomena indicate that the blades involved in the most severe friction change depending on the different maneuvering flight conditions.

Claims

1. A method for dynamic analysis of a rotor system involving blade rubbing, characterized in that, Based on the dynamic response calculated by the excitation load, the influence of frictional torque caused by rubbing on the rotational speed is considered to obtain the instantaneous rotational speed and calculate the blade azimuth angle. The excitation load is applied cyclically to simulate the unbalanced effect of the rotor system with blade rubbing. The centrifugal deformation of the blade is calculated to obtain the vibration displacement of the blade tip node and the vibration displacement of the blade root node. The difference between the two is used as the blade deformation and the blade position vector in the rotating coordinate system is calculated. The blade position vector is then projected onto the global coordinate system after coordinate transformation. Considering coating wear, the blade position vector projected onto the global coordinate system is superimposed with the position vector from the disk center to the blade root and the position vector of the rotation axis to calculate the blade tip position vector and the blade tip intrusion amount after considering coating wear, thus constructing a multi-blade rubbing force vector. A dynamic model of the rotor system with blade rubbing under maneuvering flight conditions is constructed using beam elements, shell elements, and solid elements. Dynamic calculations are performed by reducing degrees of freedom and interface separation. The rubbing dynamic response of the rotor system with blade rubbing is calculated through numerical simulation analysis.

2. The method for dynamic analysis of a rotor system involving blade rubbing according to claim 1, characterized in that, The frictional torque is calculated based on the frictional force, and the angular acceleration is calculated based on the moment of inertia of the rotor system containing blade rubbing. The instantaneous rotational speed is then calculated based on the angular acceleration. exist Under the influence of the action, the formulas for calculating angular acceleration and instantaneous rotational speed are: In the formula, Angular acceleration, For the resultant torque, For rotational inertia, The radial distance from the point of contact to the axis of rotation. The rotational speed before the collision occurred. Instantaneous rotational speed t rub For the duration of the rubbing, This is the speed feedback coefficient, with units of N·m / (rad / s). For the target speed, The total frictional force at all blade tip nodes. For the first i Friction at the leaf tip node, The number of nodes at the leaf tip. For load torque, This refers to the motor drive torque. This represents the frictional torque.

3. The method for dynamic analysis of a rotor system involving blade rubbing according to claim 2, characterized in that, blade azimuth angle α b The expression is: In the formula, This indicates the initial blade azimuth angle.

4. The method for dynamic analysis of a rotor system involving blade rubbing according to claim 2, characterized in that, The simulated unbalance effect of a rotor system with blade rubbing is specifically achieved by cyclically applying the excitation load, calculating the unbalance force of the rotor system with blade rubbing based on the calculated instantaneous rotational speed, and then performing force decomposition calculations. x Unbalanced forces F unx and y Unbalanced forces F uny ; The expression for the unbalanced force of a rotor system including blade rubbing is: In the formula, m un It is unbalanced mass. e un It is the eccentricity of unbalanced mass.

5. The method for dynamic analysis of a rotor system involving blade rubbing according to claim 1, characterized in that, The blade deformation includes the centrifugal deformation vector of the disk and the blade as a whole. and the centrifugal deformation vector of the leaf tip relative to the leaf root for: In the formula, and These are the centrifugal deformation vectors of the leaf tip and the leaf root in the local coordinate system of the leaf root, respectively. This is the inverse of the stiffness matrix of the rotor system containing blade rubbing. The centrifugal force vector of the rotor system containing blade rubbing; Under the action of frictional force, the blade moves in the local coordinate system at the leaf root. o b x b y b Dynamic deformation vector The formula for calculation is: In the formula, ( x b , y b ) represents the deformation of the leaf tip in the local coordinate system of the leaf root, ( x d , y d () represents the deformation of the leaf root in the local coordinate system of the leaf root; Calculate the position vector of the blade in the rotating coordinate system, and project the position vector of the blade in the rotating coordinate system onto the global coordinate system after coordinate transformation based on the blade azimuth angle; Leaf tip node A The position vector in the local coordinate system at the leaf root is: In the formula, Let be the position vector of the leaf tip node in the local coordinate system at the leaf root. e t This refers to the amount of tip elongation caused by thermal effects; The blade length; In a rotating coordinate system, it is denoted as Its expression is: In the formula, The angle between the local coordinate system and the rotated coordinate system at the leaf root; After transforming from the rotating coordinate system to the global coordinate system using a coordinate transformation matrix, the blade position vector in the global coordinate system... r b3 The expression is: In a rotating coordinate system, the position vector from the center of the disk to the leaf root. The expression is: In the formula, R d The radius of the disk; x di , y di () represents the deformation of the disk node in the rotating coordinate system; In the global coordinate system, the position vector from the center of the disk to the leaf root. r b2 The expression is: After considering parallel misalignment, the leaf tip node in the global coordinate system A Position vector for: ， In the formula, e x and e y for x To static parallel misalignment quantity and y Towards static parallel misalignment; and For the shaft at x The vibration displacement and rotation axis in the direction y Vibration displacement in the direction, This is the position vector of the rotation axis in the global coordinate system.

6. The method for dynamic analysis of a rotor system involving blade rubbing according to claim 5, characterized in that, Unit direction vector of frictional force and the unit direction vector of frictional force for: In the formula, , and for B The position vector of a point C The position vector of the point A The position vector of the point; the casing is uniformly and equally divided along the circumference using the finite element method. N c 1 node B Point and C The point is the leaf tip node. A Two nodes that are circumferentially adjacent, and B Point and x The included angle of the axis is greater than C Point and x Angle between axes; blade tip node A exist BC The foot of the perpendicular on the line is A 'point; , and The expression is: In the formula, ( , ), ( , )and( , )for B The coordinates of the point C The coordinates of the point and A 'The coordinates of the point; α B and α C Coated casing B The phase angle corresponding to the point and the coating on the casing C The phase angle corresponding to the point; u B and v B for B Normal displacement of the point and B Tangential displacement of the point; u C and v C for C Normal displacement of the point and C Tangential displacement of the point; The radius of the coated casing; The tip intrusion amount is calculated based on the corresponding displacement relationship, and the expression is: In the formula, For the tip of the blade without considering coating wear i The amount of intrusion into each node; All nodes at the blade tip are matched to discrete points in the coated casing based on the principle of closest phase. t The wear amount at discrete points of the coating is calculated based on the contact rubbing state at a given time. t At time +1, the wear amount at the corresponding discrete point of the coating casing is read through the phase of the blade tip node to update the blade tip intrusion amount. The blade tip intrusion amount after considering coating wear is: In the formula, e w This refers to the amount of wear on the coating. To account for the amount of blade tip intrusion after coating wear.

7. The method for dynamic analysis of a rotor system involving blade rubbing according to claim 6, characterized in that, The process of constructing the multi-blade friction force vector is as follows: Leaf tip i Friction vector at each node and friction force vector for: In the formula, k eq For the equivalent stiffness of the coated casing, For the speed correlation coefficient, μ The coefficient of friction; After weighted averaging of the friction force vector and the rubbing force vector, the first tip of the leaf... i Averaged friction force vector at each node and leaf tip i The friction force vector after averaging at each node for: , When leaf tip node A When subjected to frictional force, the coated casing experiences a reaction force. According to static equilibrium and torque equilibrium, B Point and C The vectors of the contact force and the friction force acting on the point are: ， In the formula, and for B The sum of the frictional force vectors at the points C The vector of the frictional force at the point. and for B The frictional force vector at the point and C The frictional force vector at a point, for A 'Point and C Distance between points for B Point and C The distance between points is used to calculate the tip intrusion and friction force of any blade using the same calculation method to construct a multi-blade friction force vector.

8. The method for dynamic analysis of a rotor system involving blade rubbing according to claim 7, characterized in that, The construction of the rotor system dynamics model with blade rubbing under the maneuvering flight conditions is as follows: Considering the influence of coating wear on the blade tip clearance, the rubbing force vector and friction force vector are averaged and applied to all nodes at the blade tip and the corresponding nodes of the coating casing; the rotor system dynamics model with blade rubbing under the maneuvering flight conditions is constructed using beam elements, shell elements and solid elements. In the formula, M dr , C dr , G dr and K dr The mass matrix, damping matrix, gyroscope matrix, and stiffness matrix of the rotor system containing blade rubbing are given. The damping matrix of the rotor system containing blade rubbing is calculated using Rayleigh damping. F dr This is the updated excitation load vector for a rotor system including blade rubbing, which incorporates the maneuvering flight inertial load vector. F man Misaligned load vector F m Multi-blade rubbing force vector F rub Centrifugal force vector F ce Unbalanced force vector of the rotor F u and the nonlinear contact force vector of the bearing F b ; F rub Include , , , , , ; F u Include F unx and F uny , F b Includes the contact load of double-ring inner ring angular contact ball bearings and the contact load of deep groove ball bearings; , and Let be the updated acceleration, velocity, and displacement vector of the rotor system with blade rubbing. For the system dynamics model of a rotor system including blade rubbing under maneuvering flight conditions, Newmark- β The dynamic response is calculated by first performing LU decomposition on the stiffness matrix, and then calculating the contact load of the double-half inner ring angular contact ball bearing at each time step using the Newton-Raphson algorithm. The interface separates the stiffness matrix and damping matrix into constant and fluctuating parts: In the formula, and The constant part and the wave part of the damping matrix are respectively represented by the damping matrix. and This refers to the constant part and the fluctuating part of the stiffness matrix; for the constant part of the stiffness matrix... The LU decomposition technique is used to decompose it into an upper triangular part and a lower triangular part; The dynamic model of a rotor system with blade rubbing under maneuvering flight conditions after interface separation is as follows: The dynamic model of a rotor system with blade rubbing under reduced degrees of freedom maneuvering flight conditions is as follows: In the formula, , and The mass matrix, damping matrix, and stiffness matrix of the rotor system with blade rubbing and reduced degrees of freedom are given. , and Let be the acceleration vector, velocity vector, and displacement vector of the rotor system with blade rubbing, with reduced degrees of freedom. The external excitation load vector for a rotor system with blade rubbing and reduced degrees of freedom.

Images