A method for simulating the dynamics of engine complex swept blade rubbing

Abstract

The application discloses a kind of complex sweep blade's rub-impact dynamics simulation method for engine, it is related to engine technical field, comprising: step 1, using entity unit to carry out discretization operation to actual blade, obtain the dynamics matrix of blade, by fixed interface modal synthesis method and dynamics matrix fitting combined mode within the blade reduction dynamics matrix of variable speed range is constructed;Step 2, according to the relative position of blade-casing and three-dimensional motion relationship, based on multi-coordinate system, the expression of blade tip surface intrusion is established;Step 3, combined with the idea of displacement and load equivalence, use distributed gap to deduce rub-impact surface load, on the basis of finite element model, the dynamics equation of blade-casing rub-impact is established;Step 4, the rub-impact response of blade is solved using Newmark method, and the rub-impact stress of blade is calculated by the physical equation and geometric equation of constitutive relation.

Description

Technical Field

[0001] This invention relates to the field of engine technology, and in particular to a simulation method for the collision dynamics of complex swept blades in engines. Background Technology

[0002] Fans, compressors, and turbine blades are critical structural components in aero engines responsible for energy conversion, but they are highly susceptible to structural failure. Due to the harsh operating environment, blade creep, corrosion, fatigue, and fracture can lead to blade failure. Furthermore, because of the extremely high aerodynamic efficiency requirements of engines, the blade-casing design is near-zero clearance during cruise, while the rotor itself experiences complex vibrations that affect this clearance. Therefore, blade-casing rubbing has become a common phenomenon in aero engines, further deteriorating the blade's operating environment and significantly increasing the likelihood of various forms of blade failure in engineering. Although measures such as wear-resistant coatings have been implemented to reduce damage caused by rubbing, numerous catastrophic accidents due to blade rubbing have still occurred in recent years.

[0003] Essentially, blade rubbing is a complex mechanical problem encompassing multiple physics fields, multiple scales, and non-smooth characteristics. When blades and casings rub against each other, at the microscale, various physical phenomena exist, such as local contact wear, strong impact, and friction. At the macroscale, it can lead to highly complex coupled dynamic behaviors of the structural system, including the blade, casing, and even the rotor. Dynamic simulation techniques for blade rubbing began with the interaction between the rotor disk and stationary stator components, providing a preliminary understanding of the moving disk-stationary blade rubbing phenomenon. However, the actual blade structure determines that the mechanical characteristics of blade-casing rubbing are fundamentally different from those of moving disk-stationary blade rubbing. Besides this, the most common technique uses beam-plate theory to analyze blade friction. While these techniques can be used to study the basic dynamic characteristics of blade and casing structures, they still neglect the complex configuration of blades in actual aero-engines, leading to inaccurate simulations. To more realistically reflect blade rubbing characteristics, recent research has focused on three-dimensional blade rubbing simulation techniques. However, existing simulation methods often suffer from large model sizes and high computational costs. Meanwhile, existing technologies have shortcomings in modeling the three-dimensional friction force of blades, with most still using point contact models to approximate the friction load.

[0004] As modern aero-engines demand increasingly higher design precision, simulation technology must become closer to reality. By considering the actual structure of the blades in collision analysis, establishing corresponding simulation methods, analyzing the dynamic characteristics of blades under collision failure, and conducting in-depth research on the failure mechanism of blade collision, we can provide a theoretical basis for suppressing blade friction failure and ensuring the smooth operation of the engine, which has significant engineering value. Summary of the Invention

[0005] The purpose of this invention is to provide a simulation method for the collision dynamics of complex swept blades in engines, which solves the problems of large model size, large computational load, and defects in the three-dimensional collision force modeling method of existing technologies.

[0006] To achieve the above objectives, this invention provides a simulation method for the collision dynamics of complex swept blades in engines, comprising the following steps:

[0007] Step 1: Construct the blade reduction dynamics matrix; Discretize the actual blade using solid elements to obtain the blade's dynamics matrix, and construct the blade reduction dynamics matrix within the variable speed range by combining the fixed interface modal synthesis method and dynamics matrix fitting.

[0008] Step 2: Establish the expression for the blade tip intrusion amount; Based on the relative position and three-dimensional motion relationship between the blade and the casing, establish the expression for the blade tip intrusion amount based on multiple coordinate systems;

[0009] Step 3: Establish the dynamic equations for blade-casing rubbing; combining the idea of ​​displacement and load equivalence, derive the rubbing surface load using distributed clearance, and establish the dynamic equations for blade-casing rubbing based on the finite element model.

[0010] Step 4: Calculate blade rubbing stress; use the Newmark method to solve for the blade rubbing response, and calculate the blade rubbing stress using the physical and geometric equations of the constitutive relation.

[0011] Preferably, the process of constructing the blade reduction dynamics matrix in step 1 is as follows:

[0012] S11. Discretize the actual blade to obtain the solid element finite element model and dynamic matrix of the blade. Assuming that the solid element finite element model of the blade has n nodes, the blade dynamic equations are constructed as follows:

[0013] ;

[0014] In the formula, This is the quality matrix; Here is the stiffness matrix; The damping matrix obtained after applying Rayleigh damping is, i.e. ,in A coefficient proportional to mass. A coefficient proportional to stiffness; External load on the blade; It is the displacement vector; Let u be the first derivative of u with respect to time, representing the velocity vector; Let be the second derivative of u with respect to time, and denote the acceleration vector;

[0015] S12. Due to the stiffening effect generated by the high-speed rotation of the aero-engine rotor system, the stiffness matrix of the blade is affected by the engine speed. Under the influence of the rotation effect, the dynamic equation can be expressed as a function of the rotational speed. Related formats:

[0016] ;

[0017] To avoid repeated calculations of the stiffness matrix at different speeds, a polynomial is used across the speed range. Fitting the stiffness matrix internally:

[0018] ;

[0019] In the formula, This represents the 0th term of the fitted stiffness matrix. This represents the first term of the fitted stiffness matrix. This represents the second term, which is the fitted stiffness matrix. , , The expressions are:

[0020] ;

[0021] S13. Based on the modal synthesis method, the dynamic matrix is ​​reduced in dimension. Combining the above fitting process of the dynamic matrix, the reduced dynamic matrix under the stiffening effect is derived; let... j Indicates the number of degrees of freedom in the interface. i Indicates the number of internal degrees of freedom. j and i The relational expression is:

[0022] ;

[0023] In the formula, Indicates the number of nodes;

[0024] The dynamic equation is expressed as follows:

[0025] ;

[0026] In the formula, This represents the top-left mass block matrix. This represents the top-right mass block matrix. This represents the lower left mass block matrix. This represents the lower right mass block matrix. , express , The derivative with respect to time, This represents the damping block matrix in the upper left corner. This represents the damping block matrix in the upper right corner. This represents the lower left damping block matrix. This represents the lower right damping block matrix. , express , The derivative with respect to time, This represents the upper left stiffness block matrix. Top-right stiffness block matrix Lower left stiffness block matrix The lower right stiffness block matrix, This represents the displacement component vector corresponding to the internal degrees of freedom. This represents the displacement component vector corresponding to the physical degrees of freedom. Indicates external force;

[0027] Modal synthesis with a fixed interface is used to reduce the degrees of freedom of the modal substructure. The transformation relationship between physical coordinates and modal coordinates is as follows:

[0028] ;

[0029] In the formula, The principal mode matrix is ​​normalized to its mass. For the constraint mode matrix, and ; A zero matrix; It is the identity matrix; Principal modal coordinates; Represents a reduced matrix; This represents the reduced displacement vector; i , j , k These represent the number of internal degrees of freedom, the number of interface degrees of freedom, and the number of retained modes, respectively.

[0030] Multiply both sides of the dynamic equation ,get:

[0031] ;

[0032] In the formula, This represents the second derivative of the reduced displacement vector with respect to time. This represents the derivative of the reduced displacement vector with respect to time. This represents the reduced displacement vector;

[0033] Multiply both ends by Expanding the subsequent dynamic equations yields matrix elements with the following correspondence:

[0034] ;

[0035] In the formula, , This represents the reduced mass matrix; , This represents the reduced damping matrix; , This represents the reduced stiffness matrix; Represents the identity matrix. This represents the reduced top-right mass block matrix. This represents the reduced lower left mass block matrix. This represents the reduced lower right mass block matrix. This represents the second derivative of the principal modal coordinates with respect to time. This represents the reduced top-left damping block matrix. This represents the reduced upper right damping block matrix. This represents the reduced lower left damping block matrix. This represents the reduced lower right damping block matrix. This represents the derivative of the principal modal coordinates with respect to time. This represents the reduced top-left stiffness block matrix. This represents the stiffness block matrix in the upper right corner after scaling. Represents the principal mode coordinates. This represents the reduced external force vector. All reduced dynamic matrices are of order k+j.

[0036] Preferably, the process of establishing the expression for the blade tip intrusion amount in step 2 is as follows:

[0037] Under three-dimensional friction, the tip gap and intrusion conditions vary at different points on the blade tip. The entire blade tip is miniaturized, and the distribution gap at the tip is derived kinematically using these miniature elements as units; a coordinate system is then established. XYZO , X R -Y R -Z R -O R , X r - Y r -Z r -O r Used to describe the motion of blades during the rubbing process; coordinate system XYZO It is a globally fixed coordinate system. In the globally fixed coordinate system, the origin is... O Taken from the very front end of the engine shaft. ZThe axis coincides with the engine axis; coordinate system X R -Y R -Z R -O R It is a global rotating coordinate system at the center of the wheel in a static state. X R -Y R -Z R -O R middle, Z R Axis and Global Fixed Coordinate System XYZO of Z Axis coincidence, O R The origin is located on the engine axis, and X R -O R -Y R Located at the leaf tip and Z R Global rotating coordinate system in a plane with perpendicular axes X R -Y R -Z R -O R The rotational speed is the same as the rotational speed ω in the initial state. X R shaft and X Axis parallel; coordinate system X r -Y r -Z r -O r It is a locally rotating coordinate system fixed to the center of the roulette wheel. X r -Y r -Z r -O r center, origin O r Located at the center of the wheel, each coordinate axis is relative to the global rotating coordinate system. X R -Y R -ZR -O R Corresponding parallel;

[0038] The parameters of the blade and casing structure include the casing radius. roulette radius Leaf length , Representing the axial coordinate, when the blade thickness is small and the airflow channel where the blade is located varies significantly along the axial direction, the blade length can be considered to change only with the axial position. In subsequent derivations, the independent variable representing the axial coordinate is omitted for simplicity. z However, it should be noted that the blade vibration parameters, geometric parameters, etc., will change with the change of the axial coordinate z;

[0039] S21. Based on the initial position of the blade, rotate any infinitesimal element on the blade tip in a local coordinate system from its stationary state. X r -Y r - Z r -O r The position vector in is denoted as ,in, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector. This indicates that the position vector of any infinitesimal element in a stationary state is determined by the bladed disk size parameters and is located in a locally rotating coordinate system. and and infinitesimal elements in a locally rotating coordinate system X r -Y r -Z r -O r The relative positions are determined; under the condition that the blade vibrates due to unsteady excitation load, the blade vibration displacement is transferred in a local rotating coordinate system. X r -Y r -Z r -O r The middle is recorded as , This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-component of the vector; under the condition of blade vibration, the coordinate system of any infinitesimal element at the blade tip in the locally rotating coordinate system is obtained. X r -Y r -Z r - O r The position vector below for:

[0040] ;

[0041] In the formula, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector;

[0042] S22, Based on the local rotating coordinate system X r -Y r -Z r -O r With global rotating coordinate system X R -Y R -Z R -O R The relationship is obtained so that any infinitesimal element at the leaf tip can be used in the global rotating coordinate system. X R -Y R -Z R -O R The position vector below for:

[0043] ;

[0044] In the formula, It is a locally rotated coordinate system X r -Y r -Z r -O r origin O r In the global rotating coordinate system X R -Y R-Z R -O R The position vector in the image, when there is no whirl in the impeller rotor. When the impeller rotor experiences whirl, , and It is the vortex vector of the impeller rotor center in X R shaft and Y R The projected coordinates of the axis;

[0045] S23, Leaf tip in global fixed coordinate system XYZO The position vector below Represented as:

[0046] ;

[0047] In the formula, It is the global rotating coordinate system in the initial state. X R -Y R -Z R -O R With global fixed coordinate system XYZO The axial distance between the origins Let be the rotation transformation matrix between two coordinates. The expression for the rotation transformation matrix is:

[0048] ;

[0049] In the formula, Indicates time;

[0050] S24. Establish a global fixed coordinate system. XYZO of X axis, Y axis, Z The unit vectors of the axes are respectively , k, the leaf tip position vector Recorded as:

[0051] ;

[0052] In the formula, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector;

[0053] Under conditions of rotor whirling and blade tip vibration, the amount of intrusion between the blade tip element and the casing at any axial position. for:

[0054] ;

[0055] In the formula, Indicates the radius of the casing.

[0056] Preferably, the process of establishing the dynamic equations for blade-casing rubbing in step 3 is as follows:

[0057] S31. Treating blade rubbing as an equivalent result of surface-to-surface contact between the blade and the casing, and introducing rubbing pressure to address the distribution gap at the blade tip, the normal rubbing force applies to the entire area of ​​the blade tip. A Integrating the normal contact pressure load, the normal friction force The expression is:

[0058] ;

[0059] In the formula, The normal contact pressure generated when any micro-element at the blade tip intrudes into the casing; This represents the frictional stiffness per unit area at the location of the micro-element. Let be the amount of intrusion between the blade tip and the casing at this micro-element location; here, it is assumed that the normal rubbing load direction of any micro-element is consistent; theoretically, the corresponding rubbing load at any micro-element... There are differences, but Depending on the mechanical properties of the casing structure, The change in the value is negligible, and... Treating it as a constant, we obtain the normal contact friction force. The expression is:

[0060] ;

[0061] The friction stiffness between the blade and the casing is denoted as... Friction stiffness per unit area The expression is:

[0062] ;

[0063] S32. Based on Coulomb's law of friction, the tangential contact stress corresponding to any infinitesimal element at the blade tip is obtained, and a collision load model for the collision analysis of the blade and casing in three-dimensional space is constructed:

[0064] ;

[0065] In the formula, Indicates tangential friction load. Indicates the coefficient of friction at the contact point;

[0066] S33. Due to the relatively high stiffness of the casing, the dynamic characteristics of the casing are ignored here. The casing does not vibrate; therefore, the blade tip intrusion depends only on the blade vibration. The frictional normal pressure load on the blade tip element... for:

[0067] ;

[0068] In the formula, For a global fixed coordinate system XYZO The unit normal vector of the lower blade tip element and the casing at the contact point, and the direction of the friction normal force. X r -O r -Y r When the planes are parallel, the unit normal vector is obtained. The expression is:

[0069] ;

[0070] S34. Obtain the tangential friction load of the contact friction according to Coulomb's friction law. The expression is as follows:

[0071] ;

[0072] In the formula, The coefficient of friction at the contact point is... For the normal load modulus, It is the unit tangential vector of the relative velocity direction between the blade tip element and the casing contact point;

[0073] S35. Based on the expressions for the normal pressure load and the tangential friction load, the tip micro-element in the global fixed coordinate system is obtained. XYZO The three components of the impact load applied to the object are:

[0074] ;

[0075] ;

[0076] ;

[0077] In the formula, , , express , , The derivative with respect to time;

[0078] S36. Based on the components of the friction load, in the global fixed coordinate system XYZO The rubbing load of the blade tip micro-element is denoted as follows: According to the local rotating coordinate system X r -Y r -Z r -O r With global fixed coordinate system XYZO The transformation matrix between them yields the result in the locally rotated coordinate system. X r -Y r -Z r -O r Lower friction load vector for:

[0079] ;

[0080] In the formula, This represents the transpose of T. Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component;

[0081] Any element at the blade tip in a locally rotating coordinate system X r -Y r -Z r -O r The frictional force received by the lower body for:

[0082] ;

[0083] In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component, The area of ​​the blade tip element is considered to be the friction load within any element when a denser mesh is used. By treating the discrete unit as a infinitesimal element, the frictional force is obtained. for:

[0084] ;

[0085] The friction force on the discrete element is applied equally to each node of the element. The friction force is applied to each node in an average manner. ni The first discrete unit i 1 node s The number of nodes representing the discrete element is obtained as follows:

[0086] ;

[0087] In the formula, Represents the frictional force at the nodes. Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component;

[0088] ni The corresponding physical degrees of freedom in the reduced model m Group, corresponding reduction force vector Represented as:

[0089] ;

[0090] In the formula, k +3 m -2, k +3 m -1, k +3 m The number of rows representing the vectors is used to set the friction force of the entire leaf tip. :

[0091] ;

[0092] The dynamic equations for blade-casing rubbing are established as follows:

[0093] .

[0094] Preferably, the process of calculating the blade friction stress in step 4 is as follows:

[0095] Based on the governing equations, the Newmark numerical integration method is used to solve the problem. The temporal motion law of the entire blade is obtained by reducing the theory. Combined with the finite element theory, the blade friction force is calculated through physical and geometric equations. The blade structure has irregular and non-uniform characteristics. The blade model given in this invention is divided into tetrahedral elements. Hexahedral elements can be regarded as a combination of multiple tetrahedral elements, which are mathematically equivalent.

[0096] S41. Construct any tetrahedral element of the blade in a locally rotated coordinate system. Xr -Y r -Z r -O r The element linear displacement function , and The displacements of each node in the element are:

[0097] ;

[0098] In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component;

[0099] The coordinates of each node are obtained from the dynamic equations:

[0100] ;

[0101] In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component;

[0102] The element linear displacement function is expressed as:

[0103] ;

[0104] ;

[0105] ;

[0106] In the formula, The spatial x-axis position variable within the cell, The spatial y-direction position variable within the cell. This represents the spatial z-axis position variable within the cell. This represents the x-direction displacement of node 1 in the element. This represents the x-direction displacement of node 2 in the element. This represents the x-direction displacement of node 3 in the element. This represents the x-direction displacement of node 4 in the element. This represents the y-direction displacement of node 1 in the element. This represents the y-direction displacement of node 2 in the element. This represents the y-direction displacement of node 3 in the element. This represents the y-direction displacement of node 4 in the element. This represents the z-direction displacement of node 1 in the element. This represents the z-direction displacement of node 2 in the element. This represents the z-direction displacement of node 3 in the element. This represents the z-direction displacement of node 4 in the element. H It is a determinant composed of node coordinates. Let H represent the algebraic cofactor of the element in the a-th row and b-th column of the determinant. The expression for H is:

[0107] ;

[0108] In the formula, This represents the initial x-position of node 1 in the cell. This represents the initial x-position of node 2 in the element. This indicates the initial x-position of node 3 in the element. This indicates the initial x-position of node 4 in the element. This indicates the initial y-position of node 1 in the cell. This indicates the initial y-position of node 2 in the cell. This indicates the initial y-position of node 3 in the cell. This indicates the initial y-position of node 4 in the element. This represents the initial z-axis position of node 1 in the cell. This represents the initial z-axis position of node 2 in the element. This represents the initial z-axis position of node 3 in the element. This indicates the initial z-axis position of node 4 in the cell;

[0109] S42. Position coordinates are obtained through a linear displacement function. Differentiation yields the element strain :

[0110] ;

[0111] In the formula, Indicates strain in the x-direction. Indicates strain in the y-direction. Indicates strain in the z-direction. Indicates shear strain in the xy direction. Indicates the shear strain in the yz direction. Indicates the shear strain in the zx direction;

[0112] S43, through the unit constitutive matrix Establish the stress-strain relationship of the unit element to obtain the friction stress. for:

[0113] ;

[0114] In the formula, the constitutive matrix With unit elastic modulus E Compared to Poisson v Related:

[0115] .

[0116] Therefore, this invention employs the aforementioned simulation method for the collision dynamics of complex swept-over compressor blades. For compressor blades with complex swept-over characteristics, it establishes a collision dynamics model and nonlinear solution techniques to make the simulation results closer to the actual collision situation. It addresses the complex motion problem considering the spatial coupling of the blade tip after considering the three-dimensional configuration, and establishes a multi-coordinate system transformation to determine the spatial clearance distribution during the blade tip-casing collision process. It proposes the concept of collision surface load and combines it with load equivalence theory to construct the collision force in the dynamic equations. The blade modeling integrates model reduction and dynamic matrix fitting methods, ensuring that the overall scheme simultaneously considers solution efficiency and computational stability. Based on constitutive relations, physical equations, and geometric equations, it reconstructs the blade collision stress, more intuitively reflecting the weak points of the blade during the collision process.

[0117] This invention can be used to study the vibration characteristics and damage mechanisms of blade-casing rubbing in modern aero engines, enriching the relevant theoretical methods of blade rubbing dynamics, and providing theoretical and simulation method support for fault analysis and troubleshooting in engineering.

[0118] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0119] Picture 1 This is an overall flowchart of the rubbing dynamics simulation method for complex swept blades of an engine according to the present invention.

[0120] Picture 2 This is a schematic diagram of the three-dimensional spatial features of the engine compressor blades according to an embodiment of the present invention;

[0121] Picture 3 This is a schematic diagram of the finite element model of the blade according to an embodiment of the present invention;

[0122] Picture 4 These are nodes representing the physical degrees of freedom of the interface in this embodiment of the invention.

[0123] Picture 5 This is a schematic diagram showing the main structural parameters and coordinate definitions of the bladed disk and casing in a static state according to an embodiment of the present invention;

[0124] Picture 6 This is a schematic diagram illustrating the kinematic relationship and coordinate definition between the impeller and the casing under rubbing conditions, according to an embodiment of the present invention.

[0125] Picture 7 This is a schematic diagram of the friction force in an embodiment of the present invention. Detailed Implementation

[0126] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0127] Please see Picture 1-Picture 7 A simulation method for the collision dynamics of complex swept blades in an engine includes the following steps:

[0128] Step 1: Construct the blade reduction dynamics matrix; Discretize the actual blade using solid elements to obtain the blade's dynamics matrix, and construct the blade reduction dynamics matrix within the variable speed range using a combination of fixed-interface modal synthesis and dynamics matrix fitting; the specific process is as follows:

[0129] S11, The compressor blade configuration of a certain engine is as follows: Picture 2 As shown, the actual blade is discretized to obtain a solid element finite element model of the blade, as follows. Picture 3 As shown, the solid element finite element model is divided using 10-node tetrahedral SOILD187 elements, totaling 3064 elements and 6279 nodes. The blade material is 1Cr11Ni2W2MoV, with a density, elastic modulus, and Poisson's ratio of 4.5 × 10⁻⁹ (t / mm³), 1.11 × 10⁵ (MPa), and 0.3, respectively. With n = 6279 blade nodes, the blade dynamic equations are constructed as follows:

[0130] ;

[0131] In the formula, It is a quality matrix; It is the stiffness matrix; It is the damping matrix obtained after applying Rayleigh damping, i.e. ; It is the external load on the blade; It is a displacement vector; It is the first derivative of u with respect to time, representing the velocity vector; It is the second derivative of u with respect to time, representing the acceleration vector; and By damping ratio and the upper and lower bounds of the frequency of concern and Decide, and The expression is as follows:

[0132] ;

[0133] ;

[0134] In this embodiment , and The values ​​are 0.1, 15000, and 100.

[0135] S12. Due to the stiffening effect generated by the high-speed rotation of the aero-engine rotor system, the stiffness matrix of the blade is affected by the engine speed. Under the influence of the rotation effect, the dynamic equation can be expressed as a function of the rotational speed. Related formats:

[0136] ;

[0137] To avoid repeated calculations of the stiffness matrix at different speeds, a polynomial is used across the speed range. Fitting the stiffness matrix internally:

[0138] ;

[0139] In the formula, This represents the 0th term of the fitted stiffness matrix. This represents the first term of the fitted stiffness matrix. This represents the second term, which is the fitted stiffness matrix. , , The expressions are:

[0140] ;

[0141] In the formula, Generally, the maximum value of the analysis rotation speed is taken; in this embodiment, it is 10000 r / min.

[0142] S13. Based on the modal synthesis method, the dynamic matrix is ​​reduced in dimension. Combining the above fitting process of the dynamic matrix, the reduced dynamic matrix under the stiffening effect is derived; let... j Indicates the number of degrees of freedom in the interface. i Indicates the number of internal degrees of freedom. j and i The relational expression is:

[0143] ;

[0144] In the formula, n represents the number of nodes.

[0145] In the embodiment, the physical degrees of freedom of the interface are selected from nodes such as Picture 4 As shown, 12 nodes are selected at the leading edge of the leaf tip. j =36, and the resulting dynamic equation is:

[0146] ;

[0147] In the formula, This represents the top-left mass block matrix. This represents the top-right mass block matrix. This represents the lower left mass block matrix. This represents the lower right mass block matrix. , express , The derivative with respect to time, This represents the damping block matrix in the upper left corner. This represents the damping block matrix in the upper right corner. This represents the lower left damping block matrix. This represents the lower right damping block matrix. , express , The derivative with respect to time, This represents the upper left stiffness block matrix. Top-right stiffness block matrix Lower left stiffness block matrix The lower right stiffness block matrix, This represents the displacement component vector corresponding to the internal degrees of freedom. This represents the displacement component vector corresponding to the physical degrees of freedom. Indicates external force;

[0148] By using fixed-interface modal synthesis to reduce the degrees of freedom of the modal substructure, the transformation relationship between physical coordinates and modal coordinates is as follows:

[0149] ;

[0150] In the formula, The principal mode matrix is ​​normalized to its mass. For the constraint mode matrix, and ; A zero matrix; It is the identity matrix; Principal modal coordinates; Represents a reduced matrix; This represents the reduced displacement vector; i , j, k These represent the number of internal degrees of freedom, the number of interface degrees of freedom, and the number of retained modes, respectively. In this example, k is 15.

[0151] Multiply both sides of the dynamic equation ,get:

[0152] ;

[0153] In the formula, This represents the second derivative of the reduced displacement vector with respect to time. This represents the derivative of the reduced displacement vector with respect to time. This represents the reduced displacement vector.

[0154] Multiply both ends by Expanding the subsequent dynamic equations yields matrix elements with the following correspondence:

[0155] ;

[0156] The unknowns in the formula have the following relationship:

[0157] ;

[0158] ;

[0159] ;

[0160] ;

[0161] ;

[0162] ;

[0163] ;

[0164] ;

[0165] In the formula, Let be the frequency of the r-th mode.

[0166] In the formula, , This represents the reduced mass matrix; , This represents the reduced damping matrix; , This represents the reduced stiffness matrix; Represents the identity matrix. This represents the reduced top-right mass block matrix. This represents the reduced lower left mass block matrix. This represents the reduced lower right mass block matrix. This represents the second derivative of the principal modal coordinates with respect to time. This represents the reduced top-left damping block matrix. This represents the reduced upper right damping block matrix. This represents the reduced lower left damping block matrix. This represents the reduced lower right damping block matrix. This represents the derivative of the principal modal coordinates with respect to time. This represents the reduced top-left stiffness block matrix. This represents the stiffness block matrix in the upper right corner after scaling. Represents the principal mode coordinates. This represents the reduced external force vector. All reduced dynamic matrices are of order k+j, i.e., 51.

[0167] When considering rotational effects The variation matrix varies with rotational speed. It also changes with rotational speed, resulting in a change in the dynamic matrix. , and In reality, all are affected by the rotational speed of the analysis, which can be expressed as , and .

[0168] The blade modal frequencies were obtained using a reduced model and compared with the results of the original finite element model, as shown in Table 1. The results show that the error in the first six modal frequencies of the reduced blade is basically within 0.05%, indicating the reliability and effectiveness of the reduction method.

[0169] Table 1 Vibration mode frequencies after blade reduction

[0170]

[0171] Step 2: Establish the expression for the blade tip intrusion amount; based on the relative position and three-dimensional motion relationship between the blade and the casing, establish the expression for the blade tip intrusion amount based on a multi-coordinate system; the specific process is as follows:

[0172] Under three-dimensional friction, the tip gap and intrusion conditions vary at different points on the blade tip. The entire blade tip is miniaturized, and the distribution gap at the tip is derived kinematically using these miniature elements as units; a coordinate system is then established. XYZO , X R -Y R -Z R -O R ,X r - Y r -Z r -O r This is used to describe the movement of blades during the rubbing process, such as Picture 5 and Picture 6 As shown; coordinate system XYZO It is a globally fixed coordinate system. In the globally fixed coordinate system, the origin is... O Taken from the very front end of the engine shaft. Z The axis coincides with the engine axis; coordinate system X R -Y R -Z R -O R It is a global rotating coordinate system at the center of the wheel in a static state. X R -Y R -Z R -O R middle, Z R Axis and Global Fixed Coordinate System XYZO of Z Axis coincidence, O R The origin is located on the engine axis, and X R -O R -Y R Located at the leaf tip and Z R Global rotating coordinate system in a plane with perpendicular axes X R -Y R -Z R -O R The rotational speed is the same as the rotational speed ω in the initial state. X R shaft and X Axis parallel; coordinate system X r -Y r -Z r -O rIt is a locally rotating coordinate system fixed to the center of the roulette wheel. X r -Y r -Z r -O r center, origin O r Located at the center of the wheel, each coordinate axis is relative to the global rotating coordinate system. X R -Y R -Z R -O R Corresponding parallels.

[0173] Blade and casing structural parameters include casing radius roulette radius and leaf length , Representing the axial coordinate, considering the small blade thickness and the significant axial variation of the airflow channel where the blade is located, the blade length can be assumed to change only with the axial position. In subsequent derivations, the independent variable z representing the axial coordinate is omitted for simplicity. It should be noted that blade vibration parameters, geometric parameters, etc., will vary with different axial coordinates z.

[0174] S21. Based on the initial position of the blade, any infinitesimal element on the blade tip in a static state is rotated in the local coordinate system X. r -Y r -Z r -O r The position vector is denoted as ,in, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector. This indicates that the position vector is located in a locally rotated coordinate system; this position vector can be determined based on the bladed disk size parameters. and and its local rotating coordinate system X r -Y r -Z r -O r The relative position is determined. If the blade is subjected to an unsteady excitation load and vibrates, its vibration displacement in the local rotating coordinate system X r -Y r -Z r -O r The Chinese character is represented as , This represents the x-component of the vector. This represents the y-component of the vector. Let the z-axis component of this vector be represented; then, considering blade vibration, any infinitesimal element at the blade tip in the locally rotated coordinate system X r -Y r -Z r -O r The position vector below is:

[0175] ;

[0176] In the formula, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector;

[0177] S22, Based on the locally rotated coordinate system X r -Y r -Z r -O r With global rotating coordinate system X R -Y R -Z R -O R The relationship can be used to obtain the position of any infinitesimal element at the leaf tip in the global rotating coordinate system. X R -Y R -Z R -O R The position vector below is:

[0178] ;

[0179] In the formula, It is a locally rotated coordinate system X r -Y r -Z r -O r origin O r In the global rotating coordinate system X R -Y R -Z R -O R The position vector in the image, when there is no whirl in the impeller rotor. When the impeller rotor experiences whirl, , and It is the vortex vector of the impeller rotor center in X R shaft and Y R The projected coordinates of the axis.

[0180] S23, Leaf tip in global fixed coordinate system XYZO The position vector below is represented as:

[0181] ;

[0182] In the formula, It is the global rotating coordinate system in the initial state. X R -Y R -Z R -O R With global fixed coordinate system XYZO The axial distance between the origins Let be the rotation transformation matrix between two coordinates. The expression for the rotation transformation matrix is:

[0183] ;

[0184] In the formula, Indicates time;

[0185] S24. Establish a global fixed coordinate system. XYZO of X axis, Y axis, Z The unit vectors of the axes are respectively , k, the leaf tip position vector Recorded as:

[0186] ;

[0187] In the formula, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector;

[0188] After considering rotor whirl and blade tip vibration, the intrusion between the blade tip element and the casing at any axial position is:

[0189] ;

[0190] In the formula, Indicates the radius of the casing

[0191] Step 3: Establish the dynamic equations for blade-casing rubbing; combining the concept of displacement and load equivalence, derive the rubbing surface load using distributed clearance, and establish the dynamic equations for blade-casing rubbing based on the finite element model; the specific process is as follows:

[0192] S31. Treating blade rubbing as an equivalent result of surface-to-surface contact between the blade and the casing, and introducing rubbing pressure to address the distribution gap at the blade tip, the normal rubbing force applies to the entire area of ​​the blade tip. A Integrating the normal contact pressure load, the expression for the normal friction force is:

[0193] ;

[0194] In the formula, The normal contact pressure generated when any micro-element at the blade tip intrudes into the casing; This represents the frictional stiffness per unit area at the location of the micro-element. Let be the amount of intrusion between the blade tip and the casing at this micro-element location; here, it is assumed that the normal rubbing load direction of any micro-element is consistent. Theoretically, the corresponding value at any micro-element is... There are differences, but Depending on the mechanical properties of the casing structure, The value of usually does not change much, therefore Treating it as a constant, the expression for the normal collision friction force is:

[0195] ;

[0196] The friction stiffness between the blade and the casing, expressed as the friction stiffness per unit area, is:

[0197] ;

[0198] S32. Based on Coulomb's law of friction, the tangential contact stress corresponding to any infinitesimal element at the blade tip is obtained, thus yielding a collision load model for the collision analysis of complex blades and casings in three-dimensional space:

[0199] ;

[0200] In the formula, Indicates tangential friction load. Indicates the coefficient of friction at the contact point;

[0201] S33. Due to the relatively high stiffness of the casing, the dynamic characteristics of the casing are ignored here. The casing does not vibrate, therefore the blade tip intrusion depends only on the blade vibration; the normal pressure load of the rubbing on the blade tip element is:

[0202] ;

[0203] In the formula, Let be the unit normal vector of the blade tip element and the casing at the contact point in the global coordinate system. Clearly, the direction of the frictional normal force at this point is... X r -O r -Y r For parallel planes, the expression is:

[0204] ;

[0205] S34. The tangential friction load of the contact friction is obtained according to Coulomb's law of friction, and the expression is as follows:

[0206] ;

[0207] In the formula, The coefficient of friction at the contact point is... For the normal load modulus, Let be the unit tangential vector of the relative velocity between the blade tip element and the casing at the contact point, such as... Picture 7 As shown.

[0208] Unit tangential vector The derivation process is as follows:

[0209] Based on the leaf tip position vector The absolute velocity of the blade tip element is obtained as follows:

[0210] ;

[0211] In the formula, , , express , , The derivative with respect to time; Let be the velocity vector of the blade tip element in the rotating coordinate system; The velocity of the origin of the rotating coordinate system relative to the fixed coordinate system is called the entrainment velocity. This is the time derivative of the coordinate transformation matrix; This represents eddy covariance, and the position of the global rotating coordinate system within the global fixed coordinate system; the expressions for each variable are as follows:

[0212] ;

[0213] ;

[0214] ;

[0215] In the formula, , , , , express , , , The derivative with respect to time.

[0216] Leaf tip micro-element in global fixed coordinate system XYZO The absolute velocity components are as follows:

[0217] ;

[0218] In the formula, Representing vectors The component in the x-direction, Representing vectors The component in the y-direction.

[0219] Based on the velocity vector of the blade tip element and the unit normal vector at the contact point The normal velocity component of the blade tip element at the contact point was obtained. and tangential velocity components for:

[0220] ;

[0221] Expanding the normal and tangential velocity components yields:

[0222] ;

[0223] Finally, the unit tangential vector is obtained. The expression is:

[0224] ;

[0225] S35. Based on the expressions for the normal pressure load and the tangential friction load, the tip micro-element in the global fixed coordinate system is obtained. XYZO The three components of the impact load applied to the object are:

[0226] ;

[0227] ;

[0228] ;

[0229] In the formula, , , express , , The derivative with respect to time.

[0230] S36. Based on the components of the friction load, in the global fixed coordinate system XYZO The rubbing load of the lower blade tip micro-element can also be denoted as: According to the local rotating coordinate system X r -Y r -Z r -O r With global fixed coordinate system XYZO The transformation matrix between the two, the vector of the friction load in the rotating coordinate system. for:

[0231] ;

[0232] In the formula, This represents the transpose of T. Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component.

[0233] The friction load on the blade tip micro-element is obtained, and the number of friction elements is p. For any blade tip element e j In a locally rotating coordinate system X r -Y r -Z r -O r The frictional force received by the lower body for:

[0234] ;

[0235] In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component, This represents the area of ​​the blade tip element. When a denser mesh is used, the friction load within any element can be considered as... It is constant, that is, if the discrete unit is regarded as a infinitesimal element, then:

[0236] ;

[0237] The friction force on this element needs to be applied equally to all nodes of the element. The friction force of the element is applied to each node in an average manner. ni Representing the unit's first i 1 node s The number of nodes representing a unit, i.e.:

[0238] ;

[0239] In the formula, Represents the frictional force at the nodes. Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component;

[0240] If ni corresponds to the m-th group of physical degrees of freedom in the reduced model, then the reduction force vector corresponding to that node is... It can be represented as:

[0241] ;

[0242] In the formula, k +3 m -2, k +3 m -1, k +3 m Represents the number of rows in the vector.

[0243] by Picture 4 Taking nodes 1, 2, and 3 of the physical degrees of freedom as an example, they all belong to the same element e1 and are nodes 1, 2, and 3 of that element. Their element friction force is:

[0244] ;

[0245] In the formula, This represents the area of ​​the leaf tip unit of unit e1;

[0246] Then the nodal forces of nodes 1, 2, and 3 are all:

[0247] ; ; ; ;

[0248] In the formula, , , This represents the components of the frictional force of unit e1 in the x, y, and z directions; , , This represents the components of the friction force at node 1 in element e1 in the x, y, and z directions. , , This represents the components of the friction force at node 2 in element e1 in the x, y, and z directions. , , This represents the components of the friction force at node 3 in element e1 in the x, y, and z directions.

[0249] Picture 4 Nodes 2, 3, and 4 of the physical degrees of freedom belong to the same element e2, which is still nodes 1, 2, and 3 of that element. Their element friction force is:

[0250] ;

[0251] In the formula, This represents the tip unit area of ​​unit e2;

[0252] Then, after the element force is evenly distributed, the nodal forces at nodes 2, 3, and 4 are:

[0253] ;

[0254] ; ; ;

[0255] In the formula, , , This represents the components of the friction force at node 1 in element e2 in the x, y, and z directions.

[0256] This allows us to collect the frictional force of the entire blade tip. :

[0257] ;

[0258] Therefore, the dynamic equations for blade-casing rubbing are established:

[0259] .

[0260] Step 4: Calculate blade friction stress; use the Newmark method to solve for the blade friction response, and calculate the blade friction stress using the physical and geometric equations of the constitutive relation; the specific process is as follows:

[0261] Based on the governing equations, the Newmark method, a numerical integration approach, is used to solve the problem. The temporal motion of the entire blade is then reconstructed using reduction theory. Furthermore, the finite element method is combined with physical and geometric equations to calculate the blade's friction stress. Due to the irregularity and non-uniformity of the blade structure, the blade model presented in this invention is divided into tetrahedral elements, which is considered the preferred implementation. Hexahedral elements can be viewed as combinations of multiple tetrahedral elements, which are mathematically equivalent.

[0262] S41. Let any tetrahedral element of a blade be in a locally rotated coordinate system. X r -Y r -Z r -O r The element linear displacement function , and The displacements of each node in this element are:

[0263] ;

[0264] In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component;

[0265] The position coordinates can be obtained by solving the dynamic equations. The coordinates of each node are:

[0266] ;

[0267] In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component;

[0268] but , and It can be represented as:

[0269] ;

[0270] ;

[0271] ;

[0272] In the formula, The spatial x-axis position variable within the cell, The spatial y-direction position variable within the cell. This represents the spatial z-axis position variable within the cell. This represents the x-direction displacement of node 1 in the element. This represents the x-direction displacement of node 2 in the element. This represents the x-direction displacement of node 3 in the element. This represents the x-direction displacement of node 4 in the element. This represents the y-direction displacement of node 1 in the element. This represents the y-direction displacement of node 2 in the element. This represents the y-direction displacement of node 3 in the element. This represents the y-direction displacement of node 4 in the element. This represents the z-direction displacement of node 1 in the element. This represents the z-direction displacement of node 2 in the element. This represents the z-direction displacement of node 3 in the element. This represents the z-direction displacement of node 4 in the element. H It is a determinant composed of node coordinates. Let H represent the algebraic cofactor of the element in the a-th row and b-th column of the determinant. The expression for H is:

[0273] ;

[0274] In the formula, This represents the initial x-position of node 1 in the cell. This represents the initial x-position of node 2 in the element. This indicates the initial x-position of node 3 in the element. This indicates the initial x-position of node 4 in the element. This indicates the initial y-position of node 1 in the cell. This indicates the initial y-position of node 2 in the cell. This indicates the initial y-position of node 3 in the cell. This indicates the initial y-position of node 4 in the element. This represents the initial z-axis position of node 1 in the cell. This represents the initial z-axis position of node 2 in the element. This represents the initial z-axis position of node 3 in the element. This indicates the initial z-axis position of node 4 in the cell.

[0275] S42, Element Strain Position coordinates can be obtained through linear displacement functions Differentiation yields:

[0276] ;

[0277] In the formula, Indicates strain in the x-direction. Indicates strain in the y-direction. Indicates strain in the z-direction. Indicates shear strain in the xy direction. Indicates the shear strain in the yz direction. This represents the shear strain in the zx direction.

[0278] S43. Establish the element stress-strain relationship from the element constitutive matrix D, then the element stress for:

[0279] ;

[0280] The constitutive matrix D is only related to the element elastic modulus. E Compared to Poisson v Related:

[0281] .

[0282] Therefore, the above-mentioned simulation method for the collision dynamics of complex swept blades in engines can be used to study the vibration characteristics and damage mechanisms of blade-casing collisions in modern aero engines, enriching the relevant theoretical methods of blade collision dynamics, and providing theoretical and simulation method support for fault analysis and troubleshooting in engineering.

[0283] 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 preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A simulation method for the collision dynamics of complex swept blades in an engine, characterized in that, Includes the following steps: Step 1: Construct the blade reduction dynamics matrix; Discretize the actual blade using solid elements to obtain the blade's dynamics matrix, and construct the blade reduction dynamics matrix within the variable speed range by combining the fixed interface modal synthesis method and dynamics matrix fitting. Step 2: Establish the expression for the blade tip intrusion amount; based on the relative position and three-dimensional motion relationship between the blade and the casing, establish the expression for the blade tip intrusion amount based on a multi-coordinate system. The specific expression is as follows: ； In the formula, Indicates the radius of the casing. This indicates the amount of intrusion between the blade tip micro-element and the casing at any axial position. This represents the x-component of the leaf tip position vector. This represents the y-component of the leaf tip position vector. This represents the z-component of the leaf tip position vector; Step 3: Establish the dynamic equations for blade-casing rubbing; combining the concept of displacement and load equivalence, derive the rubbing surface load using distributed clearance, and establish the dynamic equations for blade-casing rubbing based on the finite element model. The specific expressions are as follows: ； in, This represents the frictional force across the entire blade tip. This represents the second derivative of the reduced displacement vector with respect to time. This represents the derivative of the reduced displacement vector with respect to time. This represents the reduced displacement vector. This represents the reduced mass matrix. This represents the reduced damping matrix. This represents the reduced stiffness matrix; Step 4: Calculate blade rubbing stress; use the Newmark method to solve for the blade rubbing response, and calculate the blade rubbing stress using the physical and geometric equations of the constitutive relation.

2. The method for simulating the collision dynamics of complex swept blades in an engine according to claim 1, characterized in that, The process of constructing the blade reduction dynamics matrix in step 1 is as follows: S11. Discretize the actual blade to obtain the solid element finite element model and dynamic matrix of the blade. Assume the solid element finite element model of the blade has... n The blade dynamics equations are constructed based on the given nodes as follows: ； In the formula, This is the quality matrix; Here is the stiffness matrix; The damping matrix is ​​obtained after applying Rayleigh damping. ,in A coefficient proportional to mass. A coefficient proportional to stiffness; External load on the blade; It is the displacement vector; Let u be the first derivative of u with respect to time, representing the velocity vector; Let be the second derivative of u with respect to time, and denote the acceleration vector; S12. Express the dynamic equations in relation to engine speed. Related formats: ； Using polynomials in the speed range Fitting the stiffness matrix internally: ； In the formula, This represents the 0th term of the fitted stiffness matrix. This represents the first term of the fitted stiffness matrix. This represents the second term, which is the fitted stiffness matrix. , , The expressions are: ； S13, Let j Indicates the number of degrees of freedom in the interface. i Indicates the number of internal degrees of freedom. j and i The relational expression is: ； In the formula, Indicates the number of nodes; The dynamic equation is expressed as follows: ； In the formula, This represents the top-left mass block matrix. This represents the top-right mass block matrix. This represents the lower left mass block matrix. This represents the lower right mass block matrix. , express , The derivative with respect to time, This represents the damping block matrix in the upper left corner. This represents the damping block matrix in the upper right corner. This represents the lower left damping block matrix. This represents the lower right damping block matrix. , express , The derivative with respect to time, This represents the upper left stiffness block matrix. Top-right stiffness block matrix Lower left stiffness block matrix The lower right stiffness block matrix, This represents the displacement component vector corresponding to the internal degrees of freedom. This represents the displacement component vector corresponding to the physical degrees of freedom. Indicates external force; Modal synthesis with a fixed interface is used to reduce the degrees of freedom of the modal substructure. The transformation relationship between physical coordinates and modal coordinates is as follows: ； In the formula, The principal mode matrix is ​​normalized to its mass. For the constraint mode matrix, and ; A zero matrix; It is the identity matrix; Principal modal coordinates; Represents a reduced matrix; This represents the reduced displacement vector; i , j , k These represent the number of internal degrees of freedom, the number of interface degrees of freedom, and the number of retained modes, respectively. Multiply both sides of the dynamic equation ,get: ； In the formula, This represents the second derivative of the reduced displacement vector with respect to time. This represents the derivative of the reduced displacement vector with respect to time. This represents the reduced displacement vector; Multiply both ends by Expanding the subsequent dynamic equations yields matrix elements with the following correspondence: ； In the formula, , This represents the reduced mass matrix; , This represents the reduced damping matrix; , This represents the reduced stiffness matrix; Represents the identity matrix. This represents the reduced top-right mass block matrix. This represents the reduced lower left mass block matrix. This represents the reduced lower right mass block matrix. This represents the second derivative of the principal modal coordinates with respect to time. This represents the reduced top-left damping block matrix. This represents the reduced upper right damping block matrix. This represents the reduced lower left damping block matrix. This represents the reduced lower right damping block matrix. This represents the derivative of the principal modal coordinates with respect to time. This represents the reduced top-left stiffness block matrix. This represents the stiffness block matrix in the upper right corner after scaling. Represents the principal mode coordinates. This represents the reduced external force vector.

3. The method for simulating the collision dynamics of complex swept blades in an engine according to claim 2, characterized in that, The process of establishing the expression for the tip surface intrusion amount in step 2 is as follows: Establish coordinate system XYZO , X R -Y R -Z R -O R , X r -Y r -Z r -O r Used to describe the motion of blades during the rubbing process; coordinate system XYZO It is a globally fixed coordinate system. In the globally fixed coordinate system, the origin is... O Taken from the very front end of the engine shaft. Z The shaft coincides with the engine shaft; coordinate system X R -Y R -Z R -O R It is a global rotating coordinate system at the center of the wheel in a static state. X R -Y R -Z R -O R middle, Z R Axis and Global Fixed Coordinate System XYZO of Z Axis coincidence, O R The origin is located on the engine axis, and X R -O R -Y R Located at the leaf tip and Z R Global rotating coordinate system in a plane with perpendicular axes X R -Y R -Z R -O R The rotational speed is the same as the rotational speed ω in the initial state. X R shaft and X Axis parallel; coordinate system X r -Y r -Z r -O r It is a locally rotating coordinate system fixed to the center of the roulette wheel. X r -Y r -Z r -O r center, origin O r Located at the center of the wheel, each coordinate axis is relative to the global rotating coordinate system. X R -Y R -Z R -O R Corresponding parallel; The parameters of the blade and casing structure include the casing radius. roulette radius Leaf length , Represents the axial coordinate; S21. Based on the initial position of the blade, rotate any infinitesimal element on the blade tip in a local coordinate system from its stationary state. X r -Y r -Z r -O r The position vector in is denoted as ,in, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector. This indicates that the position vector of any infinitesimal element in a stationary state is determined by the bladed disk size parameters and is located in a locally rotating coordinate system. and and infinitesimal elements in a locally rotating coordinate system X r -Y r -Z r -O r The relative positions are determined; under the condition that the blade vibrates due to unsteady excitation load, the blade vibration displacement is transferred in a local rotating coordinate system. X r - Y r -Z r -O r The middle is recorded as , This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-component of the vector; under the condition of blade vibration, the coordinate system of any infinitesimal element at the blade tip in the locally rotating coordinate system is obtained. X r -Y r -Z r -O r The position vector below for: ； In the formula, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector; S22, Based on the local rotating coordinate system X r -Y r -Z r -O r With global rotating coordinate system X R -Y R -Z R -O R The relationship is obtained so that any infinitesimal element at the leaf tip can be used in the global rotating coordinate system. X R -Y R -Z R -O R The position vector below for: ； In the formula, It is a locally rotated coordinate system X r -Y r -Z r -O r origin O r In the global rotating coordinate system X R -Y R -Z R -O R The position vector in the image, when there is no whirl in the impeller rotor. When the impeller rotor experiences whirl, , and It is the vortex vector of the impeller rotor center in X R shaft and Y R The projected coordinates of the axis; S23, Leaf tip in global fixed coordinate system XYZO The position vector below Represented as: ； In the formula, It is the global rotating coordinate system in the initial state. X R -Y R -Z R -O R With global fixed coordinate system XYZO The axial distance between the origins Let be the rotation transformation matrix between two coordinates. The expression for the rotation transformation matrix is: ； In the formula, Indicates time; S24. Establish a global fixed coordinate system. XYZO of X axis, Y axis, Z The unit vectors of the axes are respectively , k, the leaf tip position vector Recorded as: ； In the formula, This represents the x-component of the vector. This represents the y-component of the vector. This represents the z-axis component of the vector; Under conditions of rotor whirling and blade tip vibration, the amount of intrusion between the blade tip element and the casing at any axial position. for: ； In the formula, Indicates the radius of the casing.

4. The method for simulating the collision dynamics of complex swept blades in an engine according to claim 3, characterized in that, The process of establishing the dynamic equations for blade-casing rubbing in step 3 is as follows: S31. Treating blade rubbing as an equivalent result of surface-to-surface contact between the blade and the casing, rubbing pressure is introduced to address the distribution gap at the blade tip, with normal rubbing force... The expression is: ； In the formula, The normal contact pressure generated when any micro-element at the blade tip intrudes into the casing; This represents the frictional stiffness per unit area at the location of the micro-element. The area of ​​the entire region at the leaf tip; Let be the amount of intrusion between the blade tip and the casing at this micro-element location; here, the normal rubbing load direction of any micro-element is considered to be consistent; theoretically, the corresponding rubbing load at any micro-element... There are differences, but Depending on the mechanical properties of the casing structure, The change in the value is negligible. Treating it as a constant, we obtain the normal contact friction force. The expression is: ； The friction stiffness between the blade and the casing is denoted as... Friction stiffness per unit area The expression is: ； S32. Based on Coulomb's law of friction, the tangential contact stress corresponding to any infinitesimal element at the blade tip is obtained, and a collision load model for the collision analysis of the blade and casing in three-dimensional space is constructed: ； In the formula, Indicates tangential friction load. Indicates the coefficient of friction at the contact point; S33. Calculate the frictional normal pressure load on the blade tip element. The expression is as follows: ； In the formula, For a global fixed coordinate system XYZO The unit normal vector of the lower blade tip element and the casing at the contact point, and the direction of the friction normal force. X r -O r -Y r When the planes are parallel, the unit normal vector is obtained. The expression is: ； S34. Obtain the tangential friction load of the contact friction according to Coulomb's friction law. The expression is as follows: ； In the formula, The coefficient of friction at the contact point is... For the normal load modulus, It is the unit tangential vector of the relative velocity direction between the blade tip element and the casing contact point; S35. Based on the expressions for the normal pressure load and the tangential friction load, the tip micro-element in the global fixed coordinate system is obtained. XYZO The three components of the impact load applied to the object are: ； ； ； In the formula, , , express , , The derivative with respect to time; S36. Based on the components of the friction load, in the global fixed coordinate system XYZO The rubbing load of the blade tip micro-element is denoted as follows: According to the local rotating coordinate system X r -Y r -Z r -O r With global fixed coordinate system XYZO The transformation matrix between them yields the result in the locally rotated coordinate system. X r -Y r -Z r -O r Lower friction load vector for: ； In the formula, This represents the transpose of T. Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component; Any element at the blade tip in a locally rotating coordinate system X r -Y r -Z r -O r The frictional force received by the lower body for: ； In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component, The area of ​​the blade tip element is the friction load within any element. By treating the discrete unit as a infinitesimal element, the frictional force is obtained. for: ； set up ni The first discrete unit i 1 node s The number of nodes representing the discrete element is obtained as follows: ； In the formula, Represents the frictional force at the nodes. Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component; ni The corresponding physical degrees of freedom in the reduced model m Group, corresponding reduction force vector Represented as: ； In the formula, k +3 m -2, k +3 m -1, k +3 m The number of rows representing the vectors is used to set the friction force of the entire leaf tip. : ； The dynamic equations for blade-casing rubbing are established as follows: 。 5. The method for simulating the collision dynamics of complex swept blades in an engine according to claim 4, characterized in that, The process of calculating the blade friction stress in step 4 is as follows: Based on the governing equations, the Newmark numerical integration method is used to solve the problem. The temporal motion law of the entire blade is obtained by reducing the theory. Combined with the finite element theory, the blade friction force is calculated through physical and geometric equations. S41. Construct any tetrahedral element of the blade in a locally rotated coordinate system. X r -Y r -Z r -O r The element linear displacement function , and The displacements of each node in the element are: ； In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component; The coordinates of each node are obtained from the dynamic equations: ； In the formula, Representing vectors The x-component, Representing vectors The y-component, Representing vectors The z-axis component; The element linear displacement function is expressed as: ； ； ； In the formula, The spatial x-axis position variable within the cell, The spatial y-direction position variable within the cell. This represents the spatial z-axis position variable within the cell. This represents the x-direction displacement of node 1 in the element. This represents the x-direction displacement of node 2 in the element. This represents the x-direction displacement of node 3 in the element. This represents the x-direction displacement of node 4 in the element. This represents the y-direction displacement of node 1 in the element. This represents the y-direction displacement of node 2 in the element. This represents the y-direction displacement of node 3 in the element. This represents the y-direction displacement of node 4 in the element. This represents the z-direction displacement of node 1 in the element. This represents the z-direction displacement of node 2 in the element. This represents the z-direction displacement of node 3 in the element. This represents the z-direction displacement of node 4 in the element. H It is a determinant composed of node coordinates. Let H represent the algebraic cofactor of the element in the a-th row and b-th column of the determinant. The expression for H is: ； In the formula, This represents the initial x-position of node 1 in the cell. This represents the initial x-position of node 2 in the element. This indicates the initial x-position of node 3 in the element. This indicates the initial x-position of node 4 in the element. This indicates the initial y-position of node 1 in the cell. This indicates the initial y-position of node 2 in the cell. This indicates the initial y-position of node 3 in the cell. This indicates the initial y-position of node 4 in the element. This represents the initial z-axis position of node 1 in the cell. This represents the initial z-axis position of node 2 in the element. This represents the initial z-axis position of node 3 in the element. This indicates the initial z-axis position of node 4 in the cell; S42. Position coordinates are obtained through a linear displacement function. Differentiation yields the element strain : ； In the formula, Indicates strain in the x-direction. Indicates strain in the y-direction. Indicates strain in the z-direction. Indicates shear strain in the xy direction. Indicates the shear strain in the yz direction. Indicates the shear strain in the zx direction; S43, through the unit constitutive matrix Establish the stress-strain relationship of the unit element to obtain the friction stress. for: ； In the formula, the constitutive matrix With unit elastic modulus E Compared to Poisson v Related, the expression is as follows: 。

Images