Rubbing dynamics simulation method for complex curved-swept blades of engine

Abstract

The invention discloses a rub-impact dynamics simulation method for complex curved and swept blades of an engine, which relates to the technical field of engines and comprises the following steps of: 1, performing discretization operation on an actual blade by using an entity unit to obtain a dynamics matrix of the blade, constructing a blade reduction dynamic matrix in a variable rotating speed range in a mode of combining a fixed interface modal synthesis method and dynamic matrix fitting; 2, according to the relation between the blade-casing relative position and the three-dimensional motion, establishing an intrusion amount expression of the blade tip molded surface of the blade based on a multi-coordinate system; step 3, in combination with a displacement and load equivalence idea, deducing a rub-impact surface load by using a distribution gap, and establishing a kinetic equation of blade-casing rub-impact on the basis of the finite element model; and 4, solving the rub-impact response of the blade by using a Newmark method, and calculating the rub-impact stress of the blade through the physical equation and the geometric equation of the constitutive relationship.

Description

Technical Field

[0001] The present invention relates to the technical field of engines, and particularly to a rubbing dynamics simulation method for complex swept and bowed engine blades. Background Art

[0002] Fan, compressor and turbine blades are key structural components responsible for energy conversion in an aeroengine, but they are extremely prone to structural failures. Since the blades are in a harsh working environment, creep, corrosion, fatigue and fracture of the blades will cause blade failures. In addition, due to the extremely high requirements of the engine for aerodynamic efficiency, at cruise conditions, the blades and the casing are designed with a near-zero clearance, while at the same time, the rotor itself has complex vibration problems that affect the clearance. Therefore, blade-casing rubbing has become a common phenomenon in aeroengines, which will further deteriorate the blade working environment and greatly increase the possibility of various forms of blade failures in engineering. Although some measures such as wear-resistant coatings have been taken to reduce the damage caused by rubbing, many catastrophic accidents caused by blade rubbing have still occurred in recent years.

[0003] Essentially, rubbing is a complex mechanical problem covering multiple physical fields, multiple scales and non-smooth characteristics. When the blade and the casing rub against each other, there are various physical phenomena such as local contact wear, strong impact and friction at the mesoscopic scale, and macroscopically, it may lead to highly complex coupled dynamic behaviors of the structural system including the blades, the casing and even the rotor. The dynamics simulation technology of blade rubbing originated from the interaction problem between the rotating disk and the fixed stator element, which provided a preliminary understanding of the rubbing phenomenon between the rotating disk and the stationary blade. However, the actual structure of the blade determines that the mechanical characteristics of blade-casing rubbing are fundamentally different from those of rotating disk-stationary blade rubbing. In addition, the most common technical solution is to use the beam plate theory to analyze the blade friction problem. Although these technologies can be used to study the basic dynamic characteristics of the blade and casing structures, they still ignore the complex configuration of the blades in an actual aeroengine, resulting in inaccurate simulation phenomena. In order to more accurately reflect the rubbing characteristics of the blades, many research works on the simulation technology of three-dimensional blade rubbing have been carried out recently. However, the existing simulation methods often have problems of large model size and large computational amount. At the same time, there are defects in the modeling method of the three-dimensional blade rubbing force in the existing technology, and most still use the point contact model to approximate the rubbing load.

[0004] With the increasing requirements of modern aeroengines for design precision, the simulation technology must be closer to the actual situation. Considering the actual structure of the blades in the rubbing analysis, establishing the corresponding simulation method, analyzing the dynamic characteristics of the blades under rubbing faults, and deeply studying the failure mechanism of blade rubbing can provide a theoretical basis for suppressing blade friction faults and the stable operation of the engine, and has important engineering value. Summary of the Invention

[0005] The object of the present invention is to provide a rubbing dynamics simulation method for complex swept and bowed blades of an engine, so as to solve the problems existing in the prior art, such as large model scale, large amount of calculation, and defects in the three-dimensional rubbing force modeling method.

[0006] To achieve the above object, the present invention provides a rubbing dynamics simulation method for complex swept and bowed blades of an engine, including the following steps:

[0007] Step 1: Construct a reduced dynamics matrix of the blade; discretize the actual blade using solid elements to obtain the dynamics matrix of the blade, and construct a reduced dynamics matrix of the blade within a variable rotational speed range by combining the fixed interface modal synthesis method and dynamics matrix fitting.

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

[0009] Step 3: Establish a dynamics equation for blade-casing rubbing; combine the idea of displacement and load equivalence, deduce the rubbing surface load using distributed gaps, and establish a dynamics equation for blade-casing rubbing based on the finite element model.

[0010] Step 4: Calculate the rubbing stress of the blade; use the Newmark method to solve the rubbing response of the blade, and calculate the rubbing stress of the blade through the physical equation and geometric equation of the constitutive relationship.

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

[0012] S11: Discretize the actual blade to obtain the solid element finite element model and dynamics matrix of the blade. Assume that the solid element finite element model of the blade has n nodes, and construct the blade dynamics equation as follows:

[0013]

[0014] In the formula, M ∈ R 3n×3n is the mass matrix; K ∈ R 3n×3n is the stiffness matrix; C ∈ R 3n×3n is the damping matrix obtained after being applied in the form of Rayleigh damping, that is, C = αM + βK, where α is the coefficient proportional to the mass, and β is the coefficient proportional to the stiffness; F ∈ R 3n×1 is the external load of the blade; u = [x1, y1, z1,..., x n , y n , z n ∈ R 3n×1 is the displacement vector; is the first derivative of u with respect to time, representing the velocity vector; is the second derivative of u with respect to time, representing the acceleration vector;

[0015] S12. Due to the stiffening effect generated by the blades under the high-speed rotation of the aero-engine rotor system, the stiffness matrix of the blades is affected by the engine speed. Under the influence of the rotational effect, the dynamic equation is expressed in a form related to the rotational speed ω:

[0016]

[0017] To avoid repeated calculations of the stiffness matrix at different rotational speeds, a polynomial is used to fit the stiffness matrix within the rotational speed range [0, ω max ]:

[0018] K(ω) = K0 + ω 2 K1 + ω 4 K2;

[0019] In the formula, K0 represents the 0th term fitting stiffness matrix, K1 represents the 1st term fitting stiffness matrix, K2 represents the 2nd term fitting stiffness matrix, and the expressions of K0, K1, and K2 are respectively:

[0020]

[0021] S13. Based on the modal synthesis method, the dynamic matrix is reduced in dimension, and combined with the above fitting process of the dynamic matrix to derive the reduced dynamic matrix under the stiffening effect; let j represent the number of interface degrees of freedom and i represent the number of internal degrees of freedom, and the relational expression between j and i is:

[0022] i = 3n - j;

[0023] In the formula, n represents the number of nodes;

[0024] The expression of the dynamic equation is obtained as:

[0025]

[0026] In the formula, M ii represents the upper left mass sub-block matrix, M ij represents the upper right mass sub-block matrix, M ji represents the lower left mass sub-block matrix, M jj represents the lower right mass sub-block matrix, represents the reciprocal of time, C ii represents the upper left damping sub-block matrix, C ij represents the upper right damping sub-block matrix, C ji represents the lower left damping sub-block matrix, C jj represents the lower right damping sub-block matrix, denote u i , the derivative of u j with respect to time, K ii denote the upper left stiffness submatrix, K ij the upper right stiffness submatrix, K ji the lower left stiffness submatrix, K jj the lower right stiffness submatrix, u i denote the displacement subvector corresponding to the internal degrees of freedom, u j denote the displacement subvector corresponding to the physical degrees of freedom, F j denote the external force;

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

[0028]

[0029] where, Φ ik ∈R i×k is the mass-normalized main modal matrix; Ψ ij ∈R i×j is the constraint modal matrix, and Ψ ij = -(K ii ) -1 K ij ; 0 jk is the 0 matrix; is the identity matrix; q k ∈R k×1 is the main modal coordinate; Γ denotes the reduction matrix; q denotes the reduced displacement vector; i, j, k respectively denote the number of internal degrees of freedom, the number of interface degrees of freedom, and the number of retained modes;

[0030] Premultiply Γ T on both sides of the dynamic equation, and we get:

[0031]

[0032] where, denotes the second derivative of the reduced displacement vector with respect to time, denotes the derivative of the reduced displacement vector with respect to time, q denotes the reduced displacement vector;

[0033] Expand the dynamic equation after premultiplying Γ T on both sides, and we get the matrix elements with the following corresponding relationships:

[0034]

[0035] where, denotes the reduced mass matrix; Denotes the reduced damping matrix; Denotes the reduced stiffness matrix; I kk Denotes the identity matrix, Denotes the reduced upper - right mass sub - matrix, Denotes the reduced lower - left mass sub - matrix, Denotes the reduced lower - right mass sub - matrix, Denotes the second - order derivative of the principal modal coordinate with respect to time, Denotes the reduced upper - left damping sub - matrix, Denotes the reduced upper - right damping sub - matrix, Denotes the reduced lower - left damping sub - matrix, Denotes the reduced lower - right damping sub - matrix, Denotes the derivative of the principal modal coordinate with respect to time, Denotes the reduced upper - left stiffness sub - matrix, Denotes the reduced upper - right stiffness sub - matrix, q k Denotes the principal modal coordinate, Denotes the reduced external force vector. The reduced dynamic matrices are all of order k + j.

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

[0037] Under three - dimensional rubbing, the tip clearances and intrusion conditions at each point of the blade tip are different. The tip is micro - elementized as a whole, and the distributed clearances of the tip are derived kinematically with the tip micro - element as the unit; establish coordinate systems X - Y - Z - O, X R -Y R -Z R -O R 、X r -Y r -Z r -O r to describe the motion of the blade during rubbing; the coordinate system X - Y - Z - O is a globally fixed coordinate system. In the globally fixed coordinate system, the origin O is taken from the front - most end of the engine axis, and the Z - axis coincides with the engine axis; the coordinate system X R -Y R -Z R -O R is the globally rotating coordinate system of the disk center in the static state. In the globally rotating coordinate system X R -Y R -Z R -O R , the Z R axis coincides with the Z - axis of the globally fixed coordinate system X - Y - Z - O, the O R origin is located on the engine axis, and X R -OR -Y R is located in the plane perpendicular to the Z axis where the blade tip is located, and the global rotation coordinate system X R -Y R -Y R -Z R -O R has the same rotational speed as the rotational speed ω. In the initial state, the X R axis is parallel to the X axis; the coordinate system X r -Y r -Z r -O r is a local rotation coordinate system fixedly connected to the center of the disk. In the local rotation coordinate system X r -Y r -Z r -O r the origin O r is located at the axis center of the disk, and each coordinate axis is parallel to the global rotation coordinate system X R -Y R -Z R -O R correspondingly;

[0038] The parameters of the blade and the casing structure include the casing radius R c (z), the disk radius R d (z) and the blade length L(z), where z represents the axial coordinate. When the blade profile thickness is small and the air flow passage where the blade is located has a large change along the axial direction, it can be considered that the blade length only changes with the axial position; in the subsequent derivation process, to simplify the expression, the independent variable z representing the axial coordinate is omitted, but it should be noted that the blade vibration parameters, geometric parameters, etc. will all change with the change of the axial coordinate z;

[0039] S21. According to the initial position of the blade, any microelement on the blade tip in the stationary state is recorded as the position vector in the local rotation coordinate system X r -Y r -Z r -O r as where, represents the x-component of the vector, represents the y-component of the vector, represents the z-component of the vector, T represents being in the local rotation coordinate system; the position vector of any microelement in the stationary state is based on the disk size parameters R d and L and the microelement in the local rotation coordinate system X r -Y r -Z r -O rDetermination of the relative position; under the condition that the blade vibrates under unsteady excitation load, the blade vibration displacement is in the local rotating coordinate system X r -Y r -Z r -O r is denoted as indicating the x-component of the vector, indicating the y-component of the vector, indicating the z-component of the vector; under the condition of blade vibration, the position vector of any micro-element at the blade tip in the local rotating coordinate system X r -Y r -Z r -O r is as follows:

[0040]

[0041] In the formula, indicating the x-component of the vector, indicating the y-component of the vector, indicating the z-component of the vector;

[0042] S22. According to the relationship between the local rotating coordinate system X r -Y r -Z r -O r and the global rotating coordinate system X R -Y R -Z R -O R the position vector of any micro-element at the blade tip in the global rotating coordinate system X R -Y R -Z R -O R is as follows:

[0043]

[0044] In the formula, Δp R is the position vector of the origin O r -Y r -Z r -O r of the local rotating coordinate system X r in the global rotating coordinate system X R -Y R -Z R -O R When there is no whirl of the disk rotor, Δp R = 0; when there is whirl of the disk rotor, and is the projection coordinates of the disk center whirling vector of the blisk rotor on the X R axis and Y R axis;

[0045] S23. The position vector p of the blade tip in the global fixed coordinate system X-Y-Z-O b1 is expressed as:

[0046]

[0047] In the formula, Δz is the axial distance between the origins of the global rotating coordinate system X R -Y R -Z R -O R and the global fixed coordinate system X-Y-Z-O, and T is the rotation transformation matrix between the two coordinate systems. The expression of the rotation transformation matrix is:

[0048]

[0049] In the formula, t represents time;

[0050] S24. Let the unit vectors of the X-axis, Y-axis, and Z-axis of the global fixed coordinate system X-Y-Z-O be i, j, and k respectively. The blade tip position vector p b1 is denoted as:

[0051]

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

[0053] Under the conditions of rotor whirling and blade tip vibration, the intrusion amount δ between the blade tip microelement and the casing at any axial position is:

[0054]

[0055] In the formula, R c represents the casing radius.

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

[0057] S31. Regarding the blade rubbing as an equivalent result after the blade and the casing are in surface-to-surface contact, the distributed clearance at the blade tip is processed by introducing the rubbing pressure. Then, the normal rubbing force is the integral of the normal contact pressure load over the entire area A on the blade tip. The expression of the normal rubbing force F n is:

[0058] F n = ∫A P n dA = ∫ A k c δdA;

[0059] In the formula, P n = k c δ is the normal contact pressure generated by the intrusion of any micro - element at the blade tip and the casing; k c is the rubbing stiffness per unit area corresponding to the position of this micro - element; δ is the intrusion amount of the blade tip at the position of this micro - element and the casing; it is considered that the directions of the normal rubbing loads of any micro - element are the same; theoretically, there are differences in k c corresponding to any micro - element, but k c depends on the mechanical properties of the casing structure, and the change law of the value of k c is ignored. Regarding k c as a constant, the normal rubbing force F n expression is:

[0060] F n = k c ∫ A δdA;

[0061] Denote the rubbing stiffness between the blade and the casing as K c , and the rubbing stiffness per unit area k c expression is:

[0062]

[0063] S32. Obtain the tangential contact stress corresponding to any micro - element at the blade tip according to Coulomb's friction law, and construct a rubbing load model for the rubbing analysis between the blade and the casing in three - dimensional space:

[0064]

[0065] In the formula, P t represents the tangential friction load, and μ represents the friction coefficient at the contact point;

[0066] S33. Since the relative stiffness of the casing is relatively large, the dynamic characteristics of the casing are ignored here, and the casing does not vibrate. Therefore, the intrusion amount of the blade tip only depends on the vibration of the blade; the normal pressure load P n received by the micro - element at the blade tip is:

[0067]

[0068] In the formula, n is the unit normal vector at the contact point between the micro - element at the blade tip and the casing in the global fixed coordinate system X - Y - Z - O, and the direction of the normal rubbing force is the same as that of X r - O r - Y rParallel to the plane, the expression for the unit normal vector n is obtained as follows:

[0069]

[0070] S34. According to Coulomb's friction law, the tangential frictional load P of the rubbing is obtained t , and the expression is as follows:

[0071] P t =-μP n t;

[0072] In the formula, μ is the friction coefficient at the contact point, and P n =||P n || is the modulus of the normal load, and t is the unit tangential vector in the direction of the tangential relative velocity at the contact point between the blade tip microelement and the casing;

[0073] S35. According to the expressions of the normal pressure load and the tangential frictional load of the rubbing, the three components of the rubbing load received by the blade tip microelement in the global fixed coordinate system X-Y-Z-O are as follows:

[0074]

[0075] In the formula, represents the derivative with respect to time;

[0076] S36. According to the components of the rubbing load, the rubbing load of the blade tip microelement is denoted as P = [P x P y P z T ; According to the transformation matrix between the local rotating coordinate system X r -Y r -Z r -O r and the global fixed coordinate system X-Y-Z-O, the rubbing load vector P r -Y r -Z r -O r in the local rotating coordinate system X r is as follows:

[0077]

[0078] In the formula, T T represents the transpose of T, represents the x-component of the vector P r , represents the y-component of the vector P r , represents the z-component of the vector P r ;

[0079] The rubbing force on any unit at the blade tip in the local rotation coordinate system X r -Y r -Z r -O r is as follows: is:

[0080]

[0081] In the formula, represents the x-component of the vector , represents the y-component of the vector , represents the z-component of the vector , A e is the area of the blade tip unit. When the adopted grid is relatively dense, it is considered that the rubbing load P in any unit r is constant. Regarding the discrete unit as a microelement body, the rubbing force is:

[0082]

[0083] The rubbing force on the discrete unit is equivalently applied to each node of the unit. The unit rubbing force is applied to each node in an average manner. Let ni represent the i-th node of the discrete unit, and s represent the number of nodes of the discrete unit, then we get:

[0084]

[0085] In the formula, represents the node rubbing force, represents the x-component of the vector , represents the y-component of the vector , represents the z-component of the vector ;

[0086] ni corresponds to the m-th group in the reduced model physical degrees of freedom, and the corresponding reduced force vector is expressed as:

[0087]

[0088] In the formula, k + 3m - 2, k + 3m - 1, k + 3m represent the number of rows of the vector. The rubbing force of the entire blade tip is obtained by assembling

[0089]

[0090] The dynamic equation for the blade-casing rubbing is established as:

[0091]

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

[0093] Based on the control equations, the Newmark numerical integration method is used for solution. According to the reduction theory, the time-series motion law of the entire blade is restored. Combining with the finite element theory, the blade rubbing force is calculated through physical equations and geometric equations. The blade structure has the characteristics of irregularity and non-uniformity. The blade model given in the present invention is divided into tetrahedral elements, and hexahedral elements can be regarded as a combination of multiple tetrahedral elements, which are mathematically equivalent.

[0094] S41. Construct the unit linear displacement function of any blade tetrahedral element in the local rotation coordinate system X r -Y r -Z r -O r and and [[ID=2,2]] The displacements of each node of the unit are:

[0095]

[0096] In the formula, represents the x-component of the vector ; represents the y-component of the vector ; represents the z-component of the vector ;

[0097] The coordinates of each node are obtained through the dynamic equation:

[0098]

[0099] In the formula, represents the x-component of the vector ; represents the y-component of the vector ; represents the z-component of the vector ;

[0100] The unit linear displacement function is expressed as:

[0101]

[0102] In the formula, x r is the spatial x-direction position variable inside the unit, y r is the spatial y-direction position variable inside the unit, z r represents the spatial z-direction position variable inside the unit. Denotes the x - direction displacement of Node 1 in the element, Denotes the x - direction displacement of Node 2 in the element, Denotes the x - direction displacement of Node 3 in the element, Denotes the x - direction displacement of Node 4 in the element, Denotes the y - direction displacement of Node 1 in the element, Denotes the y - direction displacement of Node 2 in the element, Denotes the y - direction displacement of Node 3 in the element, Denotes the y - direction displacement of Node 4 in the element, Denotes the z - direction displacement of Node 1 in the element, Denotes the z - direction displacement of Node 2 in the element, Denotes the z - direction displacement of Node 3 in the element, Denotes the z - direction displacement of Node 4 in the element, where H is the determinant composed of node coordinates, and h ab Represents the algebraic cofactor of the element in the a - th row and b - th column of the determinant. The expression of H is:

[0103]

[0104] In the formula, Denotes the initial x - direction position of Node 1 in the element, Denotes the initial x - direction position of Node 2 in the element, Denotes the initial x - direction position of Node 3 in the element, Denotes the initial x - direction position of Node 4 in the element, Denotes the initial y - direction position of Node 1 in the element, Denotes the initial y - direction position of Node 2 in the element, Denotes the initial y - direction position of Node 3 in the element, Denotes the initial y - direction position of Node 4 in the element, Denotes the initial z - direction position of Node 1 in the element, Denotes the initial z - direction position of Node 2 in the element, Denotes the initial z - direction position of Node 3 in the element, Denotes the initial z - direction position of Node 4 in the element;

[0105] S42. By taking the derivative of the position coordinates [x r y r z r T through the linear displacement function to obtain the element strain ε:

[0106]

[0107] In the formula, ε​x Denote the strain in the x direction as ε y Denote the strain in the y direction as ε z Denote the strain in the z direction as γ xy Denote the shear strain in the xy direction as γ yz Denote the shear strain in the yz direction as γ zx Denote the shear strain in the zx direction;

[0108] S43. Establish the element stress-strain relationship through the element constitutive matrix D, and obtain the rubbing stress σ as follows:

[0109] σ = Dε;

[0110] In the formula, the constitutive matrix D is related to the element elastic modulus E and Poisson's ratio v:

[0111]

[0112] Therefore, the present invention adopts the above-mentioned rubbing dynamics simulation method for complex bowed and swept blades of an engine. For compressor blades with complex bowed and swept characteristics, a rubbing dynamics model and a non-linear solution technology for such complex blades are established, making the simulation results closer to the actual rubbing situation, realizing the complex motion problem of tip space coupling after considering the three-dimensional configuration, and establishing the spatial clearance distribution during the rubbing process between the tip and the casing through multi-coordinate system transformation; the concept of rubbing surface load is proposed and the rubbing force in the dynamic equation is constructed by combining the load equivalence theory; the blade modeling integrates the model reduction and dynamic matrix fitting methods, making the overall solution take into account both the solution efficiency and the calculation stability; the rubbing stress of the blade is restored according to the constitutive relationship, physical equation and geometric equation, more intuitively reflecting the weak parts of the blade during the rubbing process.

[0113] The present invention can be used to study the vibration characteristics and damage failure mechanisms during the rubbing of blades and casings in modern aero-engines, enriches the related theoretical methods of blade rubbing dynamics, and provides theoretical and simulation method support for fault analysis and troubleshooting in engineering.

[0114] The technical solution of the present invention will be further described in detail below with reference to the drawings and embodiments. Brief Description of the Drawings

[0115] Figure 1 is the overall flowchart of a rubbing dynamics simulation method for complex bowed and swept blades of an engine according to the present invention;

[0116] Figure 2 is the schematic diagram of the three-dimensional spatial characteristics of the compressor blade of the engine in the embodiment of the present invention;

[0117] Figure 3 is the schematic diagram of the finite element model of the blade in the embodiment of the present invention;

[0118] Figure 4 The node of the interface physical degrees of freedom according to the embodiment of the present invention;

[0119] Figure 5 Schematic diagram of the main structural parameters and coordinate definitions of the bladed disk and the casing in the static state according to the embodiment of the present invention;

[0120] Figure 6 Schematic diagram of the motion relationship and coordinate definitions of the bladed disk and the casing in the rubbing state according to the embodiment of the present invention;

[0121] Figure 7 Schematic diagram of the rubbing friction force according to the embodiment of the present invention. Specific embodiments

[0122] The following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely represents selected embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention.

[0123] Please refer to Figures 1 - 7 , a rubbing dynamics simulation method for complex bowed and swept blades of an engine, comprising the following steps:

[0124] Step 1, construct a reduced dynamics matrix of the blade; discretize the actual blade using solid elements to obtain the dynamics matrix of the blade, and construct a reduced dynamics matrix of the blade within a variable speed range by combining the fixed interface modal synthesis method and dynamics matrix fitting; the specific process is as follows:

[0125] S11. The configuration of the compressor blade of a certain engine is as Figure 2 shown. Discretize the actual blade to obtain a finite element model of the solid elements of the blade, as Figure 3 shown. This finite element model of the solid elements is divided by 10-node tetrahedron SOILD187 elements, with a total of 3064 elements and 6279 nodes. The blade material is 1Cr11Ni2W2MoV, and the density, elastic modulus, and Poisson's ratio of this material are 4.5×10-9 (t / mm3), 1.11×105 (MPa), and 0.3 respectively. The number of blade nodes n = 6279, and the blade dynamics equation is constructed as:

[0126]

[0127] In the formula, M ∈ R 3n×3n is the mass matrix; K ∈ R 3n×3n is the stiffness matrix; C ∈ R 3n×3n is the damping matrix obtained after being applied in the form of Rayleigh damping, that is, C = αM + βK; F ∈ R3n×1 is the external load on the blade; u = [x1, y1, z1, …, x n , y n , z n ∈ R 3n×1 is the displacement vector; is the first-order derivative of u with respect to time, representing the velocity vector; is the second-order derivative of u with respect to time, representing the acceleration vector; α and β are determined by the damping ratio ξ and the upper and lower bounds of the concerned frequency ω i and ω j and the expressions of α and β are as follows:

[0128]

[0129] In this embodiment, ξ, ω i and ω j take the values of 0.1, 15000, and 100 respectively.

[0130] S12. Due to the stiffening effect generated by the high-speed rotation of the blade in 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 is expressed in a form related to the rotational speed ω:

[0131]

[0132] To avoid repeated calculations of the stiffness matrix at different rotational speeds, a polynomial is used to fit the stiffness matrix within the rotational speed range [0, ω max :

[0133] K(ω) = K0 + ω 2 K1 + ω 4 K2;

[0134] In the formula, K0 represents the 0th-term fitted stiffness matrix, K1 represents the 1st-term fitted stiffness matrix, K2 represents the 2nd-term fitted stiffness matrix, and the expressions of K0, K1, and K2 are respectively:

[0135]

[0136] In the formula, ω max generally takes the maximum value of the analysis rotational speed. In this embodiment, it takes 10000 r / min.

[0137] S13. Based on the modal synthesis method, the dynamic matrix is dimensionally reduced, and combined with the above fitting process of the dynamic matrix to derive the reduced dynamic matrix under the stiffening effect; let j represent the number of interface degrees of freedom and i represent the number of internal degrees of freedom, and the relational expression between j and i is:

[0138] i = 3n - j;

[0139] Wherein, n represents the number of nodes.

[0140] In the embodiment, the nodes selected for the physical degrees of freedom of the interface are as Figure 4 shown. Twelve nodes at the leading edge of the blade tip are selected. At this time, j = 36, and the obtained dynamic equation is:

[0141]

[0142] Wherein, M ii represents the upper left mass sub-block matrix, M ij represents the upper right mass sub-block matrix, M ji represents the lower left mass sub-block matrix, M jj represents the lower right mass sub-block matrix, represents the reciprocal of time, C ii represents the upper left damping sub-block matrix, C ij represents the upper right damping sub-block matrix, C ji represents the lower left damping sub-block matrix, C jj represents the lower right damping sub-block matrix, represents u i , u j the derivative with respect to time, K ii represents the upper left stiffness sub-block matrix, K ij the upper right stiffness sub-block matrix, K ji the lower left stiffness sub-block matrix, K jj the lower right stiffness sub-block matrix, u i represents the displacement sub-vector corresponding to the internal degrees of freedom, u j represents the displacement sub-vector corresponding to the physical degrees of freedom, F j represents the external force;

[0143] The fixed interface modal synthesis is adopted to reduce the degrees of freedom of the modal substructure. Then, the transformation relationship between the physical coordinates and the modal coordinates is:

[0144]

[0145] Wherein, Φ ik ∈R i×k is the mass-normalized main modal matrix; Ψ ij ∈R i×j is the constraint modal matrix, and Ψ ij = -(K ii ) -1 K ij ; 0 jk is the 0 matrix; is the identity matrix; q k ∈R k×1is the main modal coordinate; Γ represents the reduction matrix; q represents the reduced displacement vector; i, j, and k respectively represent the number of internal degrees of freedom, the number of interface degrees of freedom, and the number of retained modes. In the embodiment, k is 15.

[0146] Premultiply Γ on both ends of the dynamic equation T , and obtain:

[0147]

[0148] In the formula, represents the second derivative of the reduced displacement vector with respect to time, represents the derivative of the reduced displacement vector with respect to time, and q represents the reduced displacement vector.

[0149] Expand the dynamic equation after premultiplying Γ on both ends T , and obtain the matrix elements with the following corresponding relationships:

[0150]

[0151] The relationships of each unknown quantity in the formula are as follows:

[0152]

[0153] In the formula, λ r is the r-th modal frequency.

[0154] In the formula, represents the reduced mass matrix;

[0155] represents the reduced damping matrix;

[0156] represents the reduced stiffness matrix; I kk represents the identity matrix, represents the reduced upper-right mass sub-block matrix, represents the reduced lower-left mass sub-block matrix, represents the reduced lower-right mass sub-block matrix, represents the second derivative of the main modal coordinate with respect to time, represents the reduced upper-left damping sub-block matrix, represents the reduced upper-right damping sub-block matrix, represents the reduced lower-left damping sub-block matrix, represents the reduced lower-right damping sub-block matrix, represents the derivative of the main modal coordinate with respect to time, represents the reduced upper-left stiffness sub-block matrix, represents the reduced upper-right stiffness sub-block matrix, q kdenotes the main modal coordinates, denotes the reduced external force vector. The reduced dynamic matrices are all of order k + j, i.e., 51.

[0157] When considering the rotational effect, K(ω) changes with the rotational speed. Therefore, the transformation matrix Γ also changes with the rotational speed, resulting in the dynamic matrices and both being actually affected by the analyzed rotational speed and can be expressed as and

[0158] The reduced model is adopted to obtain the blade modal frequencies, and a comparison is made with the results of the original finite element model, as shown in Table 1. It can be seen from the results that the errors of the first six modal frequencies of the reduced blade are basically within 0.05%, indicating the reliability and effectiveness of the reduction method.

[0159] Table 1 Vibration modal frequencies of the reduced blade

[0160]

[0161]

[0162] Step 2: Establish the expression of the tip surface intrusion amount; based on the relative position and three-dimensional motion relationship between the blade and the casing, establish the expression of the tip surface intrusion amount of the blade based on multiple coordinate systems; the specific process is as follows:

[0163] The tip clearances and intrusion conditions at each tip of the blade under three-dimensional rubbing are different. The tip is discretized into micro-elements, and the distributed clearances at the tip are derived through kinematics with the tip micro-elements as the unit; establish the coordinate systems X-Y-Z-O, X R -Y R -Z R -O R 、X r -Y r -Z r -O r to describe the motion of the blade during rubbing, as shown in Figure 5 and Figure 6 ; the coordinate system X-Y-Z-O is a globally fixed coordinate system. In the globally fixed coordinate system, the origin O is taken from the front end of the engine axis, and the Z-axis coincides with the engine axis; the coordinate system X R -Y R -Z R -O R is the globally rotating coordinate system at the center of the disk in the static state. In the globally rotating coordinate system X R -Y R -Z R -O R Z RThe axis coincides with the Z-axis of the global fixed coordinate system X-Y-Z-O, and O R The origin is located on the engine axis, and X R -O R -Y R is located in the plane perpendicular to the Z R axis where the tip of the blade is located. The global rotating coordinate system X R -Y R -Z R -O R has the same rotational speed as the rotational speed ω. In the initial state, the X R axis is parallel to the X-axis; the coordinate system X r -Y r -Z r -O r is a local rotating coordinate system fixedly connected to the center of the disk. In the local rotating coordinate system X r -Y r -Z r -O r the origin O r is located at the axis center of the disk, and each coordinate axis is parallel to the global rotating coordinate system X R -Y R -Z R -O R correspondingly.

[0164] The structural parameters of the blade and the casing include the casing radius R c (z), the disk radius R d (z), and the blade length L(z). z represents the axial coordinate. Considering that the blade profile thickness is small and the air flow passage where the blade is located has a large change along the axial direction, it can be considered that the blade length only changes with the axial position. In the subsequent derivation, the independent variable z representing the axial coordinate is omitted for simplicity of expression. It should be noted that the blade vibration parameters, geometric parameters, etc. will all be different with different axial coordinates z.

[0165] S21. According to the initial position of the blade, any microelement on the tip of the blade in the static state in the local rotating coordinate system X r -Y r -Z r -O r The position vector is denoted as where represents the x-component of the vector, represents the y-component of the vector, represents the z-component of the vector, T represents being in the local rotating coordinate system; this position vector can be based on the disk size parameters R d and L and its position in the local rotating coordinate system X r -Y r -Z r -Or is determined by the relative position. If the blade vibrates under unsteady excitation loads, its vibration displacement in the local rotating coordinate system X r -Y r -Z r -O r is expressed as represents the x-component of the vector, represents the y-component of the vector, represents the z-component of the vector; then after considering the blade vibration, the position vector of any microelement at the blade tip in the local rotating coordinate system X r -Y r -Z r -O r is:

[0166]

[0167] In the formula, represents the x-component of the vector, represents the y-component of the vector, represents the z-component of the vector;

[0168] S22. According to the relationship between the local rotating coordinate system X r -Y r -Z r -O r and the global rotating coordinate system X R -Y R -Z R -O R the position vector of any microelement at the blade tip in the global rotating coordinate system X R -Y R -Z R -O R is:

[0169]

[0170] In the formula, Δp R is the position vector of the origin O r -Y r -Z r -O r of the local rotating coordinate system X r in the global rotating coordinate system X R -Y R -Z R -O R When there is no whirl of the disk rotor, Δp R = 0; when there is whirl of the disk rotor, and are the whirl vectors of the disk center of the disk rotor in the X R axis and YR The projected coordinates of the shaft.

[0171] S23. The position vector of the blade tip in the global fixed coordinate system X-Y-Z-O is expressed as:

[0172]

[0173] In the formula, Δz is the axial distance between the origins of the global rotation coordinate system X R -Y R -Z R -O R and the global fixed coordinate system X-Y-Z-O, and T is the rotation transformation matrix between the two coordinate systems. The expression of the rotation transformation matrix is:

[0174]

[0175] In the formula, t represents time;

[0176] S24. Let the unit vectors of the X-axis, Y-axis, and Z-axis of the global fixed coordinate system X-Y-Z-O be i, j, and k respectively. Denote the blade tip position vector p b1 as:

[0177]

[0178] In the formula, represents the x-component of the vector, represents the y-component of the vector, represents the z-component of the vector;

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

[0180]

[0181] In the formula, R c represents the casing radius

[0182] Step 3. Establish the dynamic equation of blade-casing rubbing; combining the idea of displacement and load equivalence, use the distributed clearance to deduce the rubbing surface load, and establish the dynamic equation of blade-casing rubbing on the basis of the finite element model. The specific process is as follows:

[0183] S31. Regarding the blade rubbing as an equivalent result after the blade and the casing are in surface-to-surface contact, introduce the rubbing pressure to process the distributed clearance at the blade tip. Then the normal rubbing force is the integral of the normal contact pressure load over the entire area A on the blade tip. The expression of the normal rubbing force is:

[0184] F n =∫ AP n dA = ∫ A k c δdA;

[0185] In the formula, P n = k c δ is the normal contact pressure generated by the intrusion of any microelement at the blade tip and the casing; k c is the rubbing stiffness per unit area corresponding to the position of this microelement; δ is the intrusion amount of the blade tip at the position of this microelement and the casing; it is considered that the directions of the normal rubbing loads of any microelement are the same. Theoretically, there are differences between the k c corresponding to any microelement, but k c depends on the structural mechanical properties of the casing, and the value of k c usually does not change much. Therefore, k c is regarded as a constant, and the expression of the normal rubbing force is obtained as:

[0186] F n = k c ∫ A δdA;

[0187] K c is the rubbing stiffness between the blade and the casing, and the expression of the rubbing stiffness per unit area is:

[0188]

[0189] S32. According to Coulomb's friction law, the tangential contact stress corresponding to any microelement at the blade tip is obtained, and thus a rubbing load model for the rubbing analysis of complex blades and casings in three-dimensional space is obtained:

[0190]

[0191] In the formula, P t represents the tangential friction load, and μ represents the friction coefficient at the contact point;

[0192] S33. Since the relative stiffness of the casing is relatively large, the dynamic characteristics of the casing are ignored here, and the casing does not vibrate. Therefore, the intrusion amount of the blade tip only depends on the vibration of the blade; the normal pressure load of the rubbing force received by the microelement at the blade tip is:

[0193]

[0194] In the formula, n is the unit normal vector of the microelement at the blade tip and the casing at the contact point in the global coordinate system. Obviously, the direction of the normal rubbing force is parallel to the X r -O r -Y r plane at this time, and its expression is:

[0195]

[0196] S34. Obtain the tangential friction load of the rub-impact according to Coulomb's friction law, and the expression is as follows:

[0197] P t = -μP n t;

[0198] where μ is the friction coefficient at the contact point, P n = ||P n || is the modulus of the normal load, and t is the unit tangential vector in the direction of the tangential relative velocity at the contact point between the blade tip element and the casing, as Figure 7 shown.

[0199] The derivation process of the unit tangential vector t is as follows:

[0200] According to the blade tip position vector p b1 , the absolute velocity of the blade tip element is obtained as:

[0201]

[0202] where represents the derivative with respect to time; is the velocity vector of the blade tip element in the rotating coordinate system; is the velocity of the origin of the rotating coordinate system relative to the fixed coordinate system, i.e., the transport velocity; is the derivative of the coordinate transformation matrix with respect to time; Δp represents the whirling motion, which is the position of the global rotating coordinate system in the global fixed coordinate system; the expressions of each variable are as follows:

[0203]

[0204] where represents x d the derivative with respect to time.

[0205] The absolute velocity components of the blade tip element in the global fixed coordinate system X - Y - Z - O are respectively:

[0206]

[0207] where x d represents the x - component of the vector Δp, and y d represents the y - component of the vector Δp.

[0208] According to the velocity vector of the blade tip element and the unit normal vector n at the contact point, the normal velocity component and the tangential velocity component of the blade tip element at the contact point are obtained as:

[0209]

[0210] The normal velocity component and the tangential velocity component are expanded to obtain:

[0211]

[0212] Finally, the expression for the unit tangential vector t is obtained as:

[0213]

[0214] S35. According to the expressions of the rubbing normal pressure load and the tangential friction load, the three components of the rubbing load received by the tip microelement in the global fixed coordinate system X-Y-Z-O are:

[0215]

[0216]

[0217] In the formula, represents the derivative with respect to time.

[0218] S36. According to the components of the rubbing load, the rubbing load of the tip microelement in the global fixed coordinate system X-Y-Z-O can also be denoted as P = [P x P y P z T . According to the transformation matrix between the local rotating coordinate system X r -Y r -Z r -O r and the global fixed coordinate system X-Y-Z-O, the vector P r of the rubbing load in the rotating coordinate system is:

[0219]

[0220] In the formula, T T represents the transpose of T, represents the x-direction component of the vector P r , represents the y-direction component of the vector P r , represents the z-direction component of the vector P r .

[0221] The rubbing load received by the tip microelement is obtained. Let the number of rubbing units be p. For any unit e j of the tip, in the local rotating coordinate system X r -Y r -Z r -O r ​The rubbing force under is:

[0222]

[0223] In the formula, represents the x - component of the vector , represents the y - component of the vector , represents the z - component of the vector , A ej is the area of the tip element. When the adopted mesh is relatively dense, it can be considered that the rubbing load P r in any element is constant, that is, regarding this discrete element as a micro - element body, so:

[0224]

[0225] The rubbing force received by this element needs to be equivalently applied to each node of the element. The element rubbing force is applied to each node in an average way. ni represents the i - th node of this element, and s represents the number of nodes of the element, that is:

[0226]

[0227] In the formula, represents the node rubbing force, represents the x - component of the vector , represents the y - component of the vector , represents the z - component of the vector ;

[0228] If ni corresponds to the m - th group in the physical degrees of freedom of the reduced model, then the reduced force vector of this node can be expressed as:

[0229]

[0230] In the formula, k + 3m - 2, k + 3m - 1, k + 3m represent the number of rows of the vector.

[0231] Taking nodes 1, 2, and 3 among the nodes to which the Figure 4 physical degrees of freedom belong as an example, they all belong to the same element e1 and are the 1st, 2nd, and 3rd nodes of this element respectively. The element rubbing force of it is:

[0232]

[0233] In the formula, A e1 represents the tip element area of element e1;

[0234] The nodal forces of nodes 1, 2, and 3 are then:

[0235]

[0236] where represents the components of the rubbing force of element e1 in the x, y, and z directions; represents the components of the rubbing force of node 1 in element e1 in the x, y, and z directions; represents the components of the rubbing force of node 2 in element e1 in the x, y, and z directions; represents the components of the rubbing force of node 3 in element e1 in the x, y, and z directions.

[0237] Figure 4 Nodes 2, 3, and 4 among the nodes belonging to the physical degrees of freedom of ,

[0237] , Figure 4 all belong to the same element e2 and are still nodes 1, 2, and 3 of this element. The element rubbing force is:

[0238]

[0239] where A e2 represents the tip element area of element e2;

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

[0241]

[0242] where represents the components of the rubbing force of node 1 in element e2 in the x, y, and z directions;

[0243] Furthermore, the rubbing force of the entire blade tip can be assembled to obtain

[0244]

[0245] Thus, the dynamic equation of blade-casing rubbing is established:

[0246]

[0247] Step 4: Calculate the blade rubbing stress; use the Newmark method to solve the rubbing response of the blade, and calculate the blade rubbing stress through the physical equation and geometric equation of the constitutive relationship; the specific process is as follows:

[0248] Based on the governing equations, the Newmark method, a numerical integration method, is used for solution, and the time-series motion law of the entire blade is restored according to the reduction theory. On this basis, combined with the finite element theory, the rubbing stress of the blade is calculated through physical equations and geometric equations. Due to the irregularity and non-uniformity of the blade structure, the blade model given in the present invention is divided into tetrahedral elements, and this is taken as the preferred implementation mode. For hexahedral elements, they can be regarded as a combination of multiple tetrahedral elements, which are mathematically equivalent.

[0249] S41. Denote the unit linear displacement function of any blade tetrahedral element in the local rotation coordinate system X r -Y r -Z r -O r as and the displacements of each node of this element are:

[0250]

[0251] wherein, represents the x-component of the vector , represents the y-component of the vector , represents the z-component of the vector ;

[0252] can be solved from the dynamic equation, and the position coordinates are [x r y r z r , T and the coordinates of each node are:

[0253]

[0254] wherein, represents the x-component of the vector , represents the y-component of the vector , represents the z-component of the vector ;

[0255] then and can be expressed as:

[0256]

[0257]

[0258] wherein, x r is the spatial x-direction position variable inside the element, y rThe spatial y-direction position variable inside the element, z r Represents the spatial z-direction position variable inside the element, Represents the x-direction displacement of Node 1 in the element, Represents the x-direction displacement of Node 2 in the element, Represents the x-direction displacement of Node 3 in the element, Represents the x-direction displacement of Node 4 in the element, Represents the y-direction displacement of Node 1 in the element, Represents the y-direction displacement of Node 2 in the element, Represents the y-direction displacement of Node 3 in the element, Represents the y-direction displacement of Node 4 in the element, Represents the z-direction displacement of Node 1 in the element, Represents the z-direction displacement of Node 2 in the element, Represents the z-direction displacement of Node 3 in the element, Represents the z-direction displacement of Node 4 in the element, H is the determinant composed of node coordinates, h ab Represents the algebraic cofactor of the element in the a-th row and b-th column of the determinant, and the expression of H is:

[0259]

[0260] In the formula, Represents the initial x-direction position of Node 1 in the element, Represents the initial x-direction position of Node 2 in the element, Represents the initial x-direction position of Node 3 in the element, Represents the initial x-direction position of Node 4 in the element, Represents the initial y-direction position of Node 1 in the element, Represents the initial y-direction position of Node 2 in the element, Represents the initial y-direction position of Node 3 in the element, Represents the initial y-direction position of Node 4 in the element, Represents the initial z-direction position of Node 1 in the element, Represents the initial z-direction position of Node 2 in the element, Represents the initial z-direction position of Node 3 in the element, Represents the initial z-direction position of Node 4 in the element.

[0261] S42. The element strain ε can be obtained by differentiating the linear displacement function with respect to the position coordinates [x r y r z r T as follows:

[0262] ​

[0263] where ε x represents the strain in the x - direction, ε y represents the strain in the y - direction, ε z represents the strain in the z - direction, γ xy represents the shear strain in the xy - direction, γ yz represents the shear strain in the yz - direction, γ zx represents the shear strain in the zx - direction.

[0264] S43. Establish the element stress - strain relationship from the element constitutive matrix D, then the element stress σ is:

[0265] σ = Dε;

[0266] where the constitutive matrix D is only related to the element elastic modulus E and Poisson's ratio v:

[0267]

[0268] Therefore, the present invention adopts the above - mentioned friction - induced vibration dynamics simulation method for complex swept - and - curved blades of an engine, which can be used to study the vibration characteristics and damage mechanisms during the blade - casing friction in modern aero - engines, enriches the theoretical methods related to blade friction - induced vibration dynamics, and provides theoretical and simulation method support for fault analysis and troubleshooting in engineering.

[0269] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions of the present invention or make equivalent substitutions, and these modifications or equivalent substitutions do not make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A rubbing dynamics simulation method for complex swept and curved blades of an engine, characterized in that, It includes the following steps: Step 1: Construct the reduced dynamics matrix of the blade; discretize the actual blade using solid elements to obtain the dynamics matrix of the blade, and construct the reduced dynamics matrix of the blade within the variable speed range by combining the fixed interface modal synthesis method and dynamics matrix fitting; Step 2: Establish the expression for the tip surface intrusion amount; Based on the relative position and three-dimensional motion relationship between the blade and the casing, establish the expression for the tip surface intrusion amount of the blade based on multiple coordinate systems; Step 3: Establish the dynamic equation of blade-casing rubbing; combining the idea of displacement and load equivalence, deduce the rubbing surface load using the distributed clearance, and establish the dynamic equation of blade-casing rubbing based on the finite element model; Step 4: Calculate the blade rubbing stress; Use the Newmark method to solve the rubbing response of the blade, and calculate the blade rubbing stress through the physical equation and geometric equation of the constitutive relationship.

2. The friction-impact dynamics simulation method for complex bowed and swept blades of an engine according to claim 1, wherein The process of constructing the reduced dynamics matrix of the blade in Step 1 is as follows: S11: Discretize the actual blade to obtain the finite element model and dynamics matrix of the solid elements of the blade. Assume that the finite element model of the blade solid elements has n nodes, and construct the blade dynamics equation as follows: where \(M\in\mathbb{R}\) 3n×3n is the mass matrix; \(K\in\mathbb{R}\) 3n×3n is the stiffness matrix; \(C\in\mathbb{R}\) 3n×3n is the damping matrix obtained after being applied in the form of Rayleigh damping, \(C = \alpha M+\beta K\), where \(\alpha\) is the coefficient proportional to the mass and \(\beta\) is the coefficient proportional to the stiffness; \(F\in\mathbb{R}\) 3n×1 is the external load on the blade; \(u = [x_1,y_1,z_1,\ldots,x\) n ,y n ,z n \(\in\mathbb{R}\) 3n×1 is the displacement vector; is the first derivative of \(u\) with respect to time, representing the velocity vector; is the second derivative of \(u\) with respect to time, representing the acceleration vector; S12: Express the dynamics equation in a form related to the engine speed ω: Using a polynomial to fit the stiffness matrix within the rotational speed range [0, ω max : K(ω) = K0 + ω 2 K1 + ω 4 K2; In the formula, K0 represents the 0th fitting stiffness matrix, K1 represents the 1st fitting stiffness matrix, K2 represents the 2nd fitting stiffness matrix, and the expressions of K0, K1, and K2 are respectively: S13: Assume that j represents the number of interface degrees of freedom and i represents the number of internal degrees of freedom. The relational expression between j and i is: i = 3n - j; In the formula, n represents the number of nodes; The obtained dynamics equation expression is: where, M ii represents the upper left mass sub - matrix, M ij represents the upper right mass sub - matrix, M ji represents the lower left mass sub - matrix, M jj represents the lower right mass sub - matrix, represents the reciprocal of time, C ii represents the upper left damping sub - matrix, C ij represents the upper right damping sub - matrix, C ji represents the lower left damping sub - matrix, C jj represents the lower right damping sub - matrix, represents u i , the derivative of u j with respect to time, K ii represents the upper left stiffness sub - matrix, K ij the upper right stiffness sub - matrix, K ji the lower left stiffness sub - matrix, K jj the lower right stiffness sub - matrix, u i represents the displacement sub - vector corresponding to the internal degrees of freedom, u j represents the displacement sub - vector corresponding to the physical degrees of freedom, F j represents the external force; Adopt fixed interface modal synthesis to reduce the degrees of freedom of the modal substructure. The transformation relationship between the physical coordinates and the modal coordinates is: where, Φ ik ∈R i×k is the mass-normalized main mode matrix; Ψ ij ∈R i×j is the constraint mode matrix, and Ψ ij = -(K ii ) -1 K ij ; 0 jk is the 0 matrix; is the identity matrix; q k ∈R k×1 is the main mode coordinate; Γ represents the reduction matrix; q represents the reduced displacement vector; i, j, k respectively represent the number of internal degrees of freedom, the number of interface degrees of freedom, and the number of retained modes; Premultiply both sides of the kinetic equation by Γ T , and we get: In the formula, represents the second derivative of the reduced displacement vector with respect to time, represents the derivative of the reduced displacement vector with respect to time, and q represents the reduced displacement vector; Pre-multiply both ends by Γ T Expand the kinetic equation after that to obtain the matrix elements with the following corresponding relationships: Wherein, represents the reduced mass matrix; represents the reduced damping matrix; represents the reduced stiffness matrix; I kk represents the identity matrix, represents the reduced upper right mass sub - matrix, represents the reduced lower left mass sub - matrix, represents the reduced lower right mass sub - matrix, represents the second derivative of the principal modal coordinate with respect to time, represents the reduced upper left damping sub - matrix, represents the reduced upper right damping sub - matrix, represents the reduced lower left damping sub - matrix, represents the reduced lower right damping sub - matrix, represents the derivative of the principal modal coordinate with respect to time, represents the reduced upper left stiffness sub - matrix, represents the reduced upper right stiffness sub - matrix, q k represents the principal modal coordinate, represents the reduced external force vector.

3. A rub-impact dynamics simulation method for complex swept and curved blades of 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 movement of the blade during the rubbing process; the coordinate system XYZO is a global fixed coordinate system. In the global fixed coordinate system, the origin O is taken from the front end of the engine axis, and the Z axis coincides with the engine axis; the coordinate system X R -Y R -Z R -O R It is the global rotating coordinate system of the wheel center in the static state. R -Y R -Z R -O R In, Z R The axis coincides with the Z axis of the global fixed coordinate system XYZO, O R The origin is located on the engine axis, and X R -O R -Y R Located at the tip of the leaf and Z R In the plane perpendicular to the axis, the global rotation coordinate system X R -Y R -Z R -O R The speed is the same as the speed ω. In the initial state, X R Axis is parallel to the X axis; coordinate system X r -Y r -Z r -O r It is a local rotating coordinate system fixed to the center of the wheel. r -Y r -Z r -O r middle, origin O r Located at the axis of the wheel, each coordinate axis is connected to the global rotating coordinate system X R -Y R -Z R -O R corresponding parallels; The parameters of the blade and casing structure include the casing radius R c (z), the disk radius R d (z) and the blade length L(z), where z represents the axial coordinate; S21. According to the initial position of the blade, denote the position vector of any microelement on the blade tip in the local rotation coordinate system X r -Y r -Z r -O r as wherein, represents the x-component of the vector, represents the y-component of the vector, represents the z-component of the vector, T represents being in the local rotation coordinate system; the position vector of any microelement in the static state is determined according to the blade disk size parameters R d and L and the relative position of the microelement in the local rotation coordinate system X r -Y r -Z r -O r ; under the condition that the blade vibrates due to unsteady excitation loads, denote the blade vibration displacement in the local rotation coordinate system X r -Y r -Z r -O r as represents the x-component of the vector, represents the y-component of the vector, represents the z-component of the vector; under the condition of blade vibration, obtain the position vector of any microelement on the blade tip in the local rotation coordinate system X r -Y r -Z r -O r as : In the formula, represents the x - component of the vector x, represents the y - component of the vector y, represents the z - component of the vector z; S22. According to the local rotation coordinate system X r -Y r -Z r -O r and the global rotation coordinate system X R -Y R -Z R -O R , the position vector of any microelement at the blade tip in the global rotation coordinate system X R -Y R -Z R -O R is as follows: That is: where, Δp R is the position vector of the origin O r in the local rotating coordinate system X r -Y r -Z r -O r in the global rotating coordinate system X R -Y R -Z R -O R When there is no whirl of the disk - blade rotor, Δp R = 0; when there is whirl of the disk - blade rotor, and are the projection coordinates of the disk - center whirl vector of the disk - blade rotor on the X R axis and Y R axis; S23. Position vector p of the blade tip in the global fixed coordinate system X-Y-Z-O b1 It is expressed as: where Δz is the axial distance between the origin of the global rotating coordinate system X R -Y R -Z R -O R and the origin of the global fixed coordinate system X-Y-Z-O, and T is the rotation transformation matrix between the two coordinate systems. The expression of the rotation transformation matrix is: In the formula, t represents time; S24. Let the unit vectors of the X-axis, Y-axis, and Z-axis of the global fixed coordinate system X-Y-Z-O be i, j, and k respectively. Denote the tip position vector p b1 as: In the formula, represents the x - component of the vector, represents the y - component of the vector, represents the z - component of the vector; Under the conditions of rotor whirl and tip vibration, the intrusion amount δ between the tip microelement at any axial position and the casing is: In the formula, R c represents the radius of the casing.

4. A rubbing dynamics simulation method for complex swept and curved blades of an engine according to claim 3, characterized in that The process of establishing the dynamic equation of blade-casing rubbing in Step 3 is as follows: S31. Treat the blade rub as an equivalent result after the blade and the casing are in surface-to-surface contact, and introduce the rubbing pressure to process the distribution clearance at the blade tip. The normal rubbing force F n The expression is as follows: F n = ∫ A P n dA = ∫ A k c δdA; Wherein, P n = k c δ is the normal contact pressure generated by the intrusion of any micro-element at the blade tip and the casing; k c is the rubbing stiffness per unit area corresponding to the position of the micro-element; δ is the intrusion amount of the blade tip at the position of the micro-element and the casing; here, the direction of the normal rubbing load of any micro-element is regarded as the same; theoretically, there are differences between the k c corresponding to any micro-element, but k c depends on the structural mechanical properties of the casing, and the change in the value of k c is negligible. Regarding k c as a constant, the expression for the normal rubbing friction force F n is: F n = k c ∫ A δdA; Denote the rubbing stiffness between the blade and the casing as K c , and the rubbing stiffness per unit area as k c The expression is as follows: S32: Obtain the tangential contact stress corresponding to any tip microelement according to the Coulomb friction law, and construct the rubbing load model for the rubbing analysis of the blade and the casing in three-dimensional space: where P t represents the tangential frictional load, and μ represents the friction coefficient at the contact point; S33. Calculate the rubbing normal pressure load P on the blade tip microelement n , and the expression is as follows: where \(n\) is the unit normal vector of the blade tip microelement and the casing at the contact point in the global fixed coordinate system \(X - Y - Z - O\), and the direction of the rubbing normal force is parallel to the \(X\) r -O r -Y r plane. The expression for the unit normal vector \(n\) is obtained as follows: S34. Obtain the tangential friction load P of the rub-impact according to Coulomb's friction law t , and the expression is as follows: P t = -μP n t; where μ is the friction coefficient at the contact point, and P n = ||P n || is the magnitude of the normal load, and t is the unit tangential vector in the direction of the tangential relative velocity at the contact point between the blade tip element and the casing; S35: According to the expressions of the rubbing normal pressure load and the tangential friction load, obtain the three components of the rubbing load received by the tip microelement in the global fixed coordinate system X-Y-Z-O as: In the formula, represents the derivative with respect to time; S36. According to each component of the rub-impact load, the rub-impact load of the blade tip microelement is denoted as P = [P x P y P z in the global fixed coordinate system X-Y-Z-O T ; According to the transformation matrix between the local rotating coordinate system X r -Y r -Z r -O r and the global fixed coordinate system X-Y-Z-O, the rubbing load vector P r -Y r -Z r -O r under the local rotating coordinate system X r is as follows: where, T T represents the transpose of T, represents the x - component of the vector P r , represents the y - component of the vector P r , represents the z - component of the vector P r ; The rubbing force on any unit at the blade tip in the local rotating coordinate system X r -Y r -Z r -O r is as follows is: In the formula, represents the x - component of the vector , represents the y - component of the vector , represents the z - component of the vector , A e is the area of the tip element, the rubbing load P in any element is r constant. Regarding the discrete element as a micro - element body, the rubbing force is: Assume that ni represents the i-th node of the discrete element and s represents the number of nodes of the discrete element, and obtain: In the formula, represents the node rubbing force, represents the x-component of the vector , represents the y-component of the vector , represents the z-component of the vector . ni corresponds to the m-th group among the physically reduced degrees of freedom, and the corresponding reduced force vector is expressed as: where k + 3m - 2, k + 3m - 1, and k + 3m represent the number of rows of the vector, and the rubbing forces at the entire blade tip are obtained by assembling Establish the dynamic equation of blade-casing rubbing as:

5. A rubbing dynamics simulation method for complex swept and curved blades of an engine according to claim 4, characterized in that, The process of calculating the blade rubbing stress in Step 4 is as follows: Based on the control equation, use the Newmark numerical integration method for solution, restore the time-series motion law of the entire blade according to the reduction theory, and combine the finite element theory to calculate the blade rubbing force through the physical equation and geometric equation; S41. Construct the element linear displacement function of any blade tetrahedral element in the local rotation coordinate system X r -Y r -Z r -O r and the displacements of each node of the element are as follows: and The displacements of each node of the element are: In the formula, represents the x - component of the vector , represents the y - component of the vector , represents the z - component of the vector . Obtain the coordinates of each node through the dynamic equation as: In the formula, represents the x-component of the vector ; represents the y-component of the vector ; represents the z-component of the vector . The unit linear displacement function is expressed as: where x r is the position variable in the x - direction of the space inside the element, y r is the position variable in the y - direction of the space inside the element, z r represents the position variable in the z - direction of the space inside the element, represents the displacement in the x - direction of Node 1 in the element, represents the displacement in the x - direction of Node 2 in the element, represents the displacement in the x - direction of Node 3 in the element, represents the displacement in the x - direction of Node 4 in the element, represents the displacement in the y - direction of Node 1 in the element, represents the displacement in the y - direction of Node 2 in the element, represents the displacement in the y - direction of Node 3 in the element, represents the displacement in the y - direction of Node 4 in the element, represents the displacement in the z - direction of Node 1 in the element, represents the displacement in the z - direction of Node 2 in the element, represents the displacement in the z - direction of Node 3 in the element, represents the displacement in the z - direction of Node 4 in the element, H is the determinant composed of node coordinates, h ab represents the algebraic cofactor of the element in the a - th row and b - th column in the determinant, and the expression of H is: In the formula, represents the initial x-position of Node No. 1 in the element, represents the initial x-position of Node No. 2 in the element, represents the initial x-position of Node No. 3 in the element, represents the initial x-position of Node No. 4 in the element, represents the initial y-position of Node No. 1 in the element, represents the initial y-position of Node No. 2 in the element, represents the initial y-position of Node No. 3 in the element, represents the initial y-position of Node No. 4 in the element, represents the initial z-position of Node No. 1 in the element, represents the initial z-position of Node No. 2 in the element, represents the initial z-position of Node No. 3 in the element, represents the initial z-position of Node No. 4 in the element; S42. Derive the element strain ε by differentiating the position coordinates [x r y r z r T with respect to the linear displacement function: where ε x represents the strain in the x direction, ε y represents the strain in the y direction, ε z represents the strain in the z direction, γ xy represents the shear strain in the xy direction, γ yz represents the shear strain in the yz direction, γ zx represents the shear strain in the zx direction; S43: Establish the unit stress-strain relationship through the unit constitutive matrix D, and obtain the rubbing stress σ as: σ = Dε; In the formula, the constitutive matrix D is related to the unit elastic modulus E and Poisson's ratio v, and the expression is as follows:

Images