Complete machine vibration analysis method for aero-engine under time-varying overload

Abstract

The invention discloses a whole machine vibration analysis method for an aero-engine under time-varying overload. Through a finite element method and a Lagrange equation, a disc unit and beam unit motion differential equation is deduced under the aircraft maneuvering flight time-varying overload condition, a dynamic finite element model of a rotor system under the time-varying overload is obtained through system final assembly, an additional damping matrix, an additional stiffness matrix and an additional exciting force vector caused by the time-varying overload are obtained, and the time-varying overload is calculated. Introducing the nonlinear force of the bearing and a finite element model of the casing, finally establishing a finite element model of the dynamics of the complete machine system of the aero-engine under the time-varying overload, and performing simulation analysis on the vibration response of the complete machine system under the time-varying overload by using a numerical integration method.

Description

Technical Field

[0001] The present invention belongs to the field of aero-engine overall vibration analysis, and specifically provides a method for analyzing the overall vibration characteristics of an aero-engine under time-varying overloads. Background Art

[0002] An aero-engine needs to bear not only the high temperature, high pressure, and high-speed loads of itself, but also the time-varying overloads during the maneuvering flights of the aircraft, such as accelerating and decelerating flights, rapid dives and pulls, small-radius turns, and turns. Time-varying overload refers to the ratio of the loads received by each direction of the aircraft's center of gravity position to the aircraft's mass that continuously changes with time during the maneuvering flight of the aircraft, including time-varying normal overload, lateral overload, and longitudinal overload. Under time-varying overloads, the additional loads acting on the engine will cause certain changes in the vibration characteristics of the aero-engine system. How to accurately predict the vibration characteristics of the overall system under time-varying overloads is the key to ensuring the structural safety and reliability of the engine.

[0003] At present, the research on the vibration characteristics of aero-engines mainly focuses on the characteristics of the engine rotor system and the characteristics of the support bearings, etc. There is less research on the characteristics of the overall engine and it is concentrated on the characteristics under general load conditions. There is less research on the vibration characteristics of the overall engine under time-varying overloads. And as the power core of the aircraft, the engine is also the prerequisite for ensuring the maneuverability of the aircraft. To avoid flight failures caused by excessive vibration of the aero-engine under time-varying overloads, it is necessary to design a method for analyzing the overall vibration of the engine under time-varying overloads. Summary of the Invention

[0004] To solve the problems existing in the prior art, the present invention discloses a method for analyzing the overall vibration under time-varying overloads, which is used for the dynamic modeling of the overall aero-engine and the analysis of the overall dynamic characteristics.

[0005] Technical Solution:

[0006] The present invention discloses a method for analyzing the overall vibration of an aero-engine under time-varying overloads, including:

[0007] Establish a coordinate system, respectively establish corresponding model units for the aero-engine casing, the main shaft and the disk of the overall aero-engine rotor system, and obtain the motion differential equations of each model unit;

[0008] According to the motion differential equations, considering the bearing nonlinear forces of the rotor system, establish a dynamic model of the overall system with bearing nonlinear forces under time-varying overloads;

[0009] Determine the rotation speed, operation period, and number of sampling points within the period of the rotor system, and initialize the displacement, velocity, and acceleration of the rotor system;

[0010] Solve the dynamic model of the whole engine system to obtain the unbalance excitation, time-varying overload excitation, and gravity excitation received by the rotor system, and then iteratively solve the displacements, velocities, and accelerations of the whole engine system at subsequent moments to predict the vibration characteristics of the whole aeroengine.

[0011] Further, the coordinate system includes:

[0012] The ground coordinate system, an inertial coordinate system fixed on the earth's surface, is used to describe the spatial position, flight speed, and acceleration of the aircraft;

[0013] The rotor coordinate system is used to describe the whirling position, speed, and acceleration of the rotor around its axis;

[0014] The disk coordinate system is fixed on the disk and rotates with the disk around the axis, and is used to describe the rotation of the disk.

[0015] Further, the model unit includes:

[0016] The disk element refers to the finite element modeling of the disk in the rotor system;

[0017] The beam element refers to the finite element modeling of the engine casing and the main shaft in the rotor system using the Timoshenko beam model.

[0018] Further, the motion differential equation of the disk element is:

[0019]

[0020] Where, q d = [x y θ x θ y T , is the generalized coordinate of the disk element, M d is the inertia matrix of the disk element without considering the maneuvering flight condition, C d is the damping matrix, Ω represents the rotational speed of the rotor, G d is the gyro matrix of the disk element without considering the maneuvering flight condition, C db represents the additional damping matrix caused by time-varying overload under maneuvering flight, K d is the stiffness matrix of the shaft section where the disk is located, K db represents the additional stiffness matrix caused by time-varying overload under maneuvering flight, F db represents the additional excitation force vector caused by time-varying overload under maneuvering flight, F du represents the unbalanced force received by the disk, Q d represents the generalized external force acting on the disk element.

[0021] Further, the motion differential equation of the beam element is: ​

[0022]

[0023] Among them, q s is the nodal displacement vector of the beam element, M s is the inertia matrix of the beam element, G s is the gyro matrix of the beam element, C sb is the additional damping matrix caused by time-varying overload under maneuvering flight, K sb is the additional stiffness matrix caused by time-varying overload under maneuvering flight, F sb is the additional excitation force caused by time-varying overload under maneuvering flight, C s is the damping matrix, Q s represents the generalized external force acting on the beam element.

[0024] Furthermore, the component of the bearing nonlinear force F B in the x and y directions in the rotor coordinate system is:

[0025]

[0026] Among them, k n is the comprehensive equivalent stiffness between the roller and the raceway, k i is the contact stiffness between the roller and the inner ring, k o is the contact stiffness between the roller and the outer ring, θ j is the rolling angle of the j-th roller, u θj is the radial elastic contact deformation of this roller, ξ is any position on the beam element, N b is the number of rollers.

[0027] Furthermore, the dynamic model of the whole machine system with bearing nonlinearity under time-varying overload is expressed as:

[0028]

[0029] Among them, M represents the inertia matrix of the whole machine system without considering the maneuvering flight condition, C represents the damping matrix of the whole machine system without considering the maneuvering flight condition, G represents the gyro matrix of the whole machine system without considering the maneuvering flight condition, K represents the stiffness matrix of the whole machine system without considering the maneuvering flight condition, F u is the unbalanced force of the rotor system, F b is the additional force caused by maneuvering flight and the vector of bearing nonlinear force, Q represents the generalized external force received by the whole machine system, including the supporting force at the supporting position of the rotor system, C b represents the additional damping matrix caused by time-varying overload under the maneuvering flight condition, K b represents the additional stiffness matrix caused by time-varying overload under the maneuvering flight condition, F b respectively represent the additional excitation force vectors caused by time-varying overload under the maneuvering flight condition.

[0030] Further, the solution process is as follows:

[0031] Determine the rotational speed Ω of the rotor system, assemble each unit matrix, and obtain the M, C, G, and K matrices of the dynamic model of the whole machine system;

[0032] Initialize the displacement u, velocity and acceleration responses of the rotor system, set the operating period T of the rotor system and the number of sampling points n within the simulation period, and obtain the time step Δt between sampling points.

[0033] Determine the calculation parameters α and β in the Newmark-β method, and the values of Δt, β, and α satisfy the following conditions:

[0034] β ≥ 0.5, α ≥ 0.25(0.5 + β) 2

[0035] β = 0.6, α = 0.3(0.5 + β) 2

[0036] Calculate the unbalanced excitation F u received by the rotor system, the time-varying overload excitation F b under maneuvering flight, and the gravity excitation G r ;

[0037] Use the Newmark-β and Newton-Raphson methods to iteratively solve for the displacement u i+1 and velocity and acceleration

[0038] Further, the iterative solution using the Newmark-β and Newton-Raphson methods includes the following steps:

[0039] Step (1) Initialize the displacement u i+1 and set the error threshold eps, the iteration increment h, and the maximum number of loop iterations cou;

[0040] Step (2) Calculate the total force F received by the whole machine system;

[0041] Step (3) Use the Newmark-β method to calculate the effective stiffness and effective load The expressions are as follows:

[0042]

[0043] Step (4) Solve for the displacement response And define F1, the expression is as follows:

[0044]

[0045] Step (5) Increment each term in it by the iteration increment h, calculate the new non-linear force, and then calculate the motion differential equation of the whole machine system to obtain Define F i , the i-th column of the Jacobian matrix is YKB(i), and the expression is as follows:

[0046]

[0047] YKB(i) = (F i - F1) / h

[0048] In step (6), use the Newton-Raphson method to calculate the new response numerical solution of the whole machine system, and the expression is as follows:

[0049]

[0050] Compare |YKB -1 × F1| with the set error threshold eps: If the former is greater than the latter, continue the iterative calculation until it is less than the latter; when the former is less than the latter, the calculation result is used as the n+1 response numerical solution at time t, and enter step (7);

[0051] If, when iterating to the maximum number of loop iterations cou, |YKB -1 × F1| is still greater than the set error threshold eps, it is determined that the iteration does not converge, mark the calculation as incomplete, and enter step (8);

[0052] Step (7) Use the response numerical solution at time t n+1 as u i+1 , and calculate the velocity and acceleration response, and the expression is as follows:

[0053]

[0054] At the same time, mark the calculation as completed and enter step (8);

[0055] Step (8) If the calculation is not completed, save the result file at this rotor rotation speed and return to step (1); if the calculation is completed, end the loop.

[0056] On the other hand, the present invention also discloses a computer device, including a memory, a processor, and a computer program stored on the memory, and the processor executes the computer program to implement the steps of the foregoing method.

[0057] Beneficial effects:

[0058] Compared with the prior art that only focuses on the characteristic research of local parts such as the engine rotor system and the pivot bearing, the present invention conducts dynamic modeling on the entire engine, comprehensively considering the time-varying mass matrix, stiffness matrix, gyro matrix, and external excitation vector under the influence of maneuvering actions. For the entire aero-engine system under time-varying overload coefficients and time-varying maneuvering parameters, by deriving the motion differential equations of the disk element and the beam element under the time-varying overload conditions of maneuvering flight, a finite element dynamic model of the entire engine under time-varying overload is established.

[0059] An aero-engine vibration analysis method for the entire engine disclosed by the present invention not only establishes a dynamic model of the rotor under time-varying overload, but also establishes dynamic models of the casing, bearings, elastic supports, etc. connected to the rotor, considers the bearing nonlinear force, and gives a vibration response solution algorithm for the entire nonlinear system and a key parameter value-taking method, realizing accurate prediction of the vibration characteristics of the entire aero-engine, and can effectively improve the safety and reliability of the engine structure. Brief description of the drawings

[0060] Figure 1 A flowchart of an aero-engine vibration analysis method for the entire engine under time-varying overload according to the present invention;

[0061] Figure 2 A schematic diagram of the spatial position relationship of three coordinate systems used in the present invention;

[0062] Figure 3 A schematic diagram of the entire aero-engine model in the embodiment of the present invention;

[0063] Figure 4 The time-varying overload coefficient of the aircraft in the embodiment of the present invention;

[0064] Figure 5 The time-varying flight angle of the aircraft in the embodiment of the present invention;

[0065] Figure 6 The time-varying flight angular velocity of the aircraft in the embodiment of the present invention;

[0066] Figure 7 The time-varying flight speed of the aircraft in the embodiment of the present invention;

[0067] Figure 8 The time-varying flight acceleration of the aircraft in the embodiment of the present invention;

[0068] Figure 9 The time-varying rotor speed of the aircraft in the embodiment of the present invention;

[0069] Figure 10 The vibration velocity in the X direction at the main bearing predicted by the method of the present invention;

[0070] Figure 11 The vibration velocity in the Y direction at the main bearing predicted by the method of the present invention. Detailed implementation manners

[0071] The present invention will be further clarified below in conjunction with the accompanying drawings and specific implementation manners. It should be understood that the following specific implementation manners are only used to illustrate the present invention and not to limit the scope of the present invention. After reading the present invention, various equivalent forms of modification of the present invention by those skilled in the art fall within the scope defined by the appended claims of this application.

[0072] The present invention discloses a method for analyzing the vibration of an entire aero-engine under time-varying overload. Specifically, through the finite element method and Lagrange equations, the differential equations of motion of the disk element and the beam element are derived for the time-varying overload conditions of maneuvering flight, a dynamic finite element model of the rotor system under time-varying overload is obtained, as well as the additional damping matrix, additional stiffness matrix, and additional excitation force vector caused by time-varying overload. The nonlinear force of the bearing and the finite element model of the casing are introduced, and finally, a dynamic finite element model of the entire rotor-bearing-casing system under time-varying overload is established, and the vibration response of the entire system under large maneuvering overload conditions is simulated and analyzed using the numerical integration method. The flow of the method of the present invention is as Figure 1 shown, and the technical solution of the present invention will be described below in conjunction with specific steps.

[0073] (1) Establish a coordinate system.

[0074] To describe the motion of the disk during maneuvering flight, a ground coordinate system, a rotor coordinate system, and a disk coordinate system are established respectively, as Figure 2 shown. Among them, the ground coordinate system OXYZ is an inertial coordinate system fixed on the earth's surface. The coordinate axis OY (normal) is upward along the direction perpendicular to the ground, OX (lateral) and OZ (longitudinal) are in the horizontal direction parallel to the ground. OXYZ is used to describe the spatial position, flight speed, and acceleration of the aircraft. The angles of the aircraft's deflection around the OX, OY, and OZ axes are called the pitch angle, track angle, and roll angle respectively; oxyz is the rotor coordinate system, and the coordinate system oxyz is used to describe the whirling position, speed, and acceleration of the rotor around its axis. It is assumed that one end of the rotor coincides with the center of mass of the aircraft, that is, the origin o is also the center of gravity of the aircraft. At this time, oxyz is also the aircraft body coordinate system; the disk coordinate system o′ξηζ is fixed on the disk and can rotate around the o′ζ axis with the disk. At the initial moment, o′ξηζ is parallel to the corresponding coordinate axes of the coordinate system oxyz. The coordinate system o′ξηζ is used to describe the rotation of the disk.

[0075] (2) Analyze the differential equation of motion of the disk element.

[0076] The translational kinetic energy of the disk is expressed as:

[0077]

[0078] wherein, m is the mass of the disk, r is the vibration displacement of the disk, and v B is the absolute velocity of the aircraft in the ground north-east coordinate system, and ω B,X , ω B,Y , ω B,Z are the angular velocity components of the aircraft (coordinate system oxyz) about the respective coordinate axes of the OXYZ coordinate system, denotes the synthesis angular velocity of ω B,X , ω B,Y , ω B,Z these three angular velocities.

[0079] The rotational kinetic energy of the disk is expressed as:

[0080]

[0081] wherein, I d is the diametral moment of inertia of the disk, I p is the polar moment of inertia of the disk, and ω ξ , ω η , ω ζ are the angular velocity components of the disk about the respective coordinate axes of the disk coordinate system o′ξηζ, and are represented by the Euler angles (φ, β, γ).

[0082] Neglecting the influence of gravity, assuming that the stiffness matrix of the shaft segment where the disk is located in the rotor system is K d , then the elastic deformation potential energy of the shaft segment at the disk position is expressed as:

[0083]

[0084] wherein, K d is a 4×4 matrix, and q d = [x y θ x θ y T , which is the generalized coordinate of the disk.

[0085] If considering that there is linear damping in the rotor system and the damping matrix is C d (4×4), the dissipated energy is expressed by using the Rayleigh function as:

[0086]

[0087] The Lagrange equation of the non-conservative system is:

[0088]

[0089] wherein, L = T t + T​r -V is the Lagrange function of the system, and q dj is the generalized coordinate q d of the j-th degree of freedom, and Q dj is the generalized force acting on this degree of freedom, which here are the components of the unbalanced force acting on the disk in each degree of freedom.

[0090] Substituting the kinetic energy, dissipated energy, and elastic potential energy into the Lagrange equation, the differential equation of motion of the disk element under maneuvering flight can be obtained, and it is written in matrix form as follows:

[0091]

[0092] where M d and G d are respectively the inertia matrix and gyro matrix of the disk without considering the maneuvering flight condition; C db and K db and F db respectively represent the additional damping matrix, additional stiffness matrix, and additional excitation force vector caused by the time-varying overload under maneuvering flight; F du is the unbalanced force acting on the disk; Q d represents the generalized external force acting on the disk.

[0093] Specifically, each matrix is expressed as follows:

[0094]

[0095]

[0096] In the formula, Ω represents the rotational speed of the rotor; e is the eccentricity of the disk; is the initial phase angle; are respectively the velocity components of the aircraft along each coordinate axis in the OXYZ coordinate system.

[0097] (3) Analyze the differential equation of motion of the beam element.

[0098] Using the Timoshenko beam considering shear deformation and rotational inertia to perform finite element modeling on the engine casing and the main shaft of the rotor system. The Timoshenko beam element contains 2 nodes, and each node has 4 degrees of freedom, which are the translational degrees of freedom along the x-axis and y-axis and the rotational degrees of freedom around the x-axis and y-axis respectively.

[0099] The nodal displacement vector q s of the beam element is expressed as:

[0100] q s = [x1 y1 θ x1 θy1 x2 y2θ x2 θ y2 T

[0101] The displacement at any ξ of the beam element is expressed in terms of the nodal displacements as follows:

[0102]

[0103]

[0104] where l s is the axial length, I s and A s are the cross-sectional moment of inertia and the cross-sectional area of the beam element respectively, κ is the shear deformation coefficient, E s is the Young's modulus of the beam element, and G sm is the shear modulus of the beam element.

[0105] The strain energy expression of the element is:

[0106]

[0107] Rewritten as:

[0108]

[0109] where q sxoz =[x s1 θ y1 x s2 θ y2 T , K s is the stiffness matrix of the element, and the expression of the element stiffness matrix is as follows:

[0110]

[0111] The stiffness matrix in the space coordinates is obtained as:

[0112]

[0113] The translational kinetic energy T st and the rotational kinetic energy T sr of the beam element are respectively:

[0114]

[0115] where θ ys +θ B,Y ≈θ ys , and ρ s is the material density of the beam element, r s ​​Indicates the relative position of the infinitesimal segment in the coordinate system oxyz, dI sd and dI sp are the diametral moment of inertia and the polar moment of inertia of the infinitesimal segment, respectively.

[0116] The elastic potential energy of the beam element is expressed as:

[0117]

[0118] Assume that the beam element is subject to external linear damping, and the damping matrix is C s (8×8), then the dissipated energy can be expressed as:

[0119]

[0120] Substituting into the Lagrange equation, the motion differential equation of the Timoshenko beam element under time-varying overload conditions can be obtained:

[0121]

[0122] In the formula, M s is the inertia matrix of the beam element, G s is the gyro matrix of the beam element, Q s represents the generalized external force acting on the beam element, C sb , K sb and F sb are the additional damping matrix, additional stiffness matrix and additional excitation force caused by time-varying overload during maneuvering flight, respectively. Their expressions are as follows:

[0123]

[0124] C sb = 2ρ s A s ω B,Z H1

[0125]

[0126]

[0127] Among them,

[0128]

[0129]

[0130] Among them,

[0131]

[0132] (4) Consider the bearing nonlinear force.

[0133] Assume that there is pure rolling between the roller and the inner and outer rings of the bearing, and the rotational speed of the cage is equal to the revolution speed of the roller. Then, at time t, the rolling angle θ of the j-th rolling element j and the radial elastic contact deformation u θj are expressed as:

[0134]

[0135] where ω bm =(ω in r1 + ω out r2) / (r1 + r2) is the revolution speed of the roller, γ represents the initial radial clearance of the bearing, and N b represents the number of rollers.

[0136] Using Hertz contact theory, the components of the bearing nonlinear force in the x and y directions are obtained as:

[0137]

[0138] where k n is the comprehensive equivalent stiffness between the roller and the raceway, and k i and k o are the contact stiffnesses between the roller and the inner and outer rings respectively.

[0139] When u θj ≤ 0, When u θj > 0, For deep groove ball bearings, ξ is generally taken as 3 / 2.

[0140] (5) Assemble the whole machine model to obtain the dynamic model of the whole machine system.

[0141] The dynamic model of the whole machine system with bearing nonlinearity under time-varying overload is expressed as follows:

[0142]

[0143] where M, C, G, and K represent the inertia matrix, damping matrix, gyro matrix, and stiffness matrix of the whole machine system without considering the maneuvering flight conditions respectively, F u is the unbalanced force of the rotor system, F b is the additional force caused by maneuvering flight and the vector of bearing nonlinear force, Q represents the generalized external force received by the whole machine system, including the supporting force at the supporting position of the rotor system. C b , K b and F b represent the additional damping matrix, additional stiffness matrix, and additional excitation force vector caused by time-varying overload under maneuvering flight conditions respectively.

[0144] (6) Solve the differential equations of motion in the dynamic model of the whole machine system to obtain the vibration prediction of the whole machine.

[0145] The present invention uses the combined method of Newmark-β and Newton-Raphson to numerically integrate and solve the differential equations of motion of the whole machine system. The solution process is as follows:

[0146] 6.1) Determine the self-rotation speed Ω of the rotor system;

[0147] 6.2) Assemble the element matrices to obtain the M, C, G, and K matrices of the dynamic model of the whole machine system;

[0148] 6.3) Initialize the displacement u, velocity and acceleration responses of the rotor system;

[0149] 6.4) Determine the operating period T of the rotor system and the number of sampling points n within the simulation period to obtain the time step Δt between sampling points. Determine the calculation parameters α and β in the Newmark-β method. When the calculation parameters satisfy β≥0.5 and α≥0.25(0.5 + β) 2 , the Newmark-β method can ensure the convergence during nonlinear calculation. Therefore, the values of Δt, β, and α are taken as follows:

[0150] β = 0.6, α = 0.3(0.5 + β) 2

[0151] 6.5) Calculate the unbalanced excitation F u acting on the rotor system, the time-varying overload excitation F b under maneuvering flight, and the gravity excitation G r ;

[0152] 6.6) Use the Newton-Raphson method to iteratively solve for the displacement u i+1 , velocity and acceleration The specific iterative calculation steps are as follows:

[0153] Step (1) Initialize the displacement u i+1 , set the error threshold eps, the iterative increment h, and the maximum number of loop iterations cou;

[0154] Step (2) Calculate the total force F acting on the whole machine system;

[0155] Step (3) Use the Newmark-β method to calculate the effective stiffness and effective load The expressions are as follows:

[0156]

[0157] Step (4) Solve the displacement response And define F1, with the expression as follows:

[0158]

[0159] Step (5) Increment each term in it by the iteration increment h, calculate the new non - linear acting force, and then calculate the motion differential equation of the whole machine system to obtain Define F i , and the i - th column of the Jacobian matrix is YKB(i), with the expression as follows:

[0160]

[0161] YKB(i) = (F i - F1) / h

[0162] Step (6) Use the Newton - Raphson method to calculate the new response numerical solution of the whole machine system, with the expression as follows:

[0163]

[0164] Compare |YKB -1 ×F1| with the set error threshold eps: If the former is greater than the latter, continue the iterative calculation until it is less than the latter; when the former is less than the latter, the calculation result is used as the response numerical solution at t n+1 , and enter Step (7);

[0165] If, when iterating to the maximum number of loop iterations cou, |YKB -1 ×F1| is still greater than the set error threshold eps, it is determined that the calculation does not converge, mark the calculation as unfinished, and enter Step (8);

[0166] Step (7) Use the response numerical solution at t n+1 as u i+1 , and calculate the velocity and acceleration responses, with the expression as follows:

[0167]

[0168] At the same time, mark the calculation as completed and enter Step (8);

[0169] Step (8) If the calculation is not completed, save the result file at this rotor self - rotation speed and return to Step (1); if the calculation is completed, end the loop.

[0170] The present invention calculates an example of an aero-engine whole machine system with a dual rotor, such as Figure 3 shown. This system consists of an inner rotor and an outer rotor, which respectively simulate the low-pressure and high-pressure rotors in the engine. The support schemes for the inner and outer rotors are 0-1-1 and 1-0-1 respectively. In addition, there are 4 disks and 4 main bearings in the system. The casing is connected to the main bearings. Disk 1 and Disk 4 are distributed on the inner rotor, and Disk 2 and Disk 3 are located on the outer rotor. Among them, Disk 3 and Disk 4 respectively simulate the high-pressure turbine and the low-pressure turbine. The casing-support stiffness is shown in Table 2, the parameters of the intermediate bearing (Support 4) are shown in Table 2, and the parameters of the elastic support are shown in Table 3.

[0171] Table 1 Casing-Support Radial Stiffness

[0172]

[0173] Table 2 Intermediate Bearing Parameters

[0174]

[0175] Table 3 SFD-Related Parameters

[0176]

[0177] Collect the time-varying overload parameters and time-varying maneuvering actions of the aircraft, including the angles, angular velocities, velocities, accelerations, and rotational speeds at each moment in different directions. The obtained data is as Figures 4 to 9 shown. According to the whole machine vibration analysis method provided by the present invention, the Figures 4 - 9 above-mentioned time-varying overload parameters and time-varying maneuvering actions are input into the whole machine system dynamics model of the present invention, and the time-varying vibration response of the aero-engine whole machine system is predicted. The results are as Figure 10 、 Figure 11 shown.

[0178] From Figure 10 、 Figure 11 it can be seen that, different from the steady-state vibration, the radial reaction forces received by the 4 main bearings in the X direction and the Y direction are time-varying reaction forces. Whether it is the amplitude or the low-frequency envelope response, the responses of each support to the time-varying parameters are quite different. The reaction forces in the X direction and the Y direction are asymmetric. In the X direction, the reaction force at Bearing 2 is the largest, while in the Y direction, the reaction force at Bearing 3 is the largest, indicating that the time-varying overload coefficient has a great influence on the reaction forces of the main bearings. If the steady-state non-time-varying method in the prior art is used, the obtained reaction force results cannot characterize the bearing loads caused by the time-varying maneuvering actions.

[0179] The design of the engine main bearing needs to consider not only the steady-state working conditions, but also the time-varying non-steady-state working conditions. A method for analyzing the overall vibration of an aeroengine under time-varying overload disclosed by the present invention performs dynamic modeling on the overall engine, comprehensively considers the time-varying mass matrix, stiffness matrix, gyro matrix, and external excitation vector under the influence of maneuvering actions. For the overall aeroengine system under time-varying overload coefficients and time-varying maneuvering parameters, the nonlinear bearing force is considered, and the time-varying vibration response of the aeroengine can be predicted more accurately. Further, statistical analysis and load spectrum compilation are carried out on these time-varying reaction forces, which can provide support for the fatigue life design of the main bearing.

Claims

1. A method for analyzing the overall vibration of an aero-engine under time-varying overload, characterized in that Including: Establish a coordinate system, respectively establish corresponding model units for the aero-engine brake, the main shaft and the disk of the aero-engine whole rotor system, and obtain the motion differential equations of each model unit; According to the motion differential equations, considering the bearing nonlinear force of the rotor system, establish a dynamic model of the whole machine system with bearing nonlinear force under time-varying overload; Determine the rotation speed, operating period, and number of sampling points within the period of the rotor system, and initialize the displacement, velocity, and acceleration of the rotor system; Solve the dynamic model of the whole machine system to obtain the unbalance excitation, time-varying overload excitation, and gravity excitation received by the rotor system, and then iteratively solve the displacement, velocity, and acceleration of the whole machine system at subsequent moments to predict the vibration characteristics of the aero-engine whole machine.

2. The aero-engine overall vibration analysis method according to claim 1, characterized in that The coordinate system includes: The ground coordinate system, an inertial coordinate system fixed on the earth's surface, used to describe the spatial position, flight speed, and acceleration of the aircraft; The rotor coordinate system, used to describe the whirling position, speed, and acceleration of the rotor around its axis; The disk coordinate system, fixed on the disk and rotating around the axis together with the disk, used to describe the rotation of the disk.

3. The aero-engine overall vibration analysis method according to claim 2, characterized in that The model unit includes: The disk unit refers to the finite element modeling of the disk in the rotor system; The beam unit refers to the finite element modeling of the engine casing and the main shaft in the rotor system using the Timoshenko beam model.

4. The method for analyzing the overall vibration of an aero-engine according to claim 3, characterized in that, The motion differential equation of the disk unit is: where q d = [x y θ x θ y T , is the generalized coordinate of the disk element, M d is the inertia matrix of the disk element without considering the maneuvering flight condition, C d is the damping matrix, Ω represents the rotational speed of the rotor, G d is the gyro matrix of the disk element without considering the maneuvering flight condition, C db represents the additional damping matrix caused by the time-varying overload under maneuvering flight, K d is the stiffness matrix of the shaft segment where the disk is located, K db represents the additional stiffness matrix caused by the time-varying overload under maneuvering flight, F db represents the additional excitation force vector caused by the time-varying overload under maneuvering flight, F du represents the unbalanced force acting on the disk, Q d represents the generalized external force acting on the disk element.​ 5. The method for analyzing the overall vibration of an aeroengine according to claim 4, wherein The motion differential equation of the beam unit is: where q s is the nodal displacement vector of the beam element, M s is the inertia matrix of the beam element, G s is the gyro matrix of the beam element, C sb is the additional damping matrix caused by time-varying overload under maneuvering flight, K sb is the additional stiffness matrix caused by time-varying overload under maneuvering flight, F sb is the additional excitation force caused by time-varying overload under maneuvering flight, C s is the damping matrix, Q s represents the generalized external force acting on the beam element.

6. The method for analyzing the overall vibration of an aeroengine according to claim 5, characterized in that, The bearing nonlinear force F B The components in the x and y directions in the rotor coordinate system are as follows: Among them, k n is the comprehensive equivalent stiffness between the roller and the raceway, k i is the contact stiffness between the roller and the inner ring, k o is the contact stiffness between the roller and the outer ring, θ j is the rolling angle of the j-th roller, u θj is the radial elastic contact deformation of this roller, ξ is any position on the beam element, N b is the number of rollers.

7. The method for analyzing the overall vibration of an aeroengine according to claim 6, characterized in that, The dynamic model of the whole machine system with bearing nonlinearity under time-varying overload is expressed as: Among them, M represents the inertia matrix of the system without considering the maneuvering flight conditions, C represents the damping matrix of the system without considering the maneuvering flight conditions, G represents the gyro matrix of the system without considering the maneuvering flight conditions, K represents the stiffness matrix of the system without considering the maneuvering flight conditions, and F u is the unbalanced force of the rotor system, and F b is the additional force caused by maneuvering flight and the nonlinear force vector of the bearing. Q represents the generalized external force received by the system, including the bearing force at the support position of the rotor system, and C b represents the additional damping matrix caused by the time-varying overload under maneuvering flight conditions, and K b represents the additional stiffness matrix caused by the time-varying overload under maneuvering flight conditions, and F b respectively represent the additional excitation force vectors caused by the time-varying overload under maneuvering flight conditions.

8. The method for analyzing the overall vibration of an aeroengine according to claim 7, wherein, The solution process is: Determine the rotation speed Ω of the rotor system, assemble the matrixes of each unit, and obtain the M, C, G, and K matrixes of the dynamic model of the whole machine system; Initialize the displacement u, velocity and acceleration responses of the rotor system. Set the operating period T of the rotor system and the number of sampling points n within the simulation period to obtain the time step Δt between sampling points. Determine the calculation parameters α and β in the Newmark-β method, and the values of △t, β, and α satisfy the following conditions: β≥0.5,α≥0.25(0.5+β) 2 β＝0.6，α＝0.3(0.5+β) 2 Calculate the unbalanced excitation force F acting on the rotor system u , the time-varying overload excitation force F during maneuvering flight b and the gravity excitation force G r ; The displacement u is iteratively solved using the Newmark-β and Newton-Raphson methods i+1 , velocity , and acceleration 9. The method for analyzing the overall vibration of an aeroengine according to claim 8, wherein The iterative solution using the Newmark-β and Newton-Raphson methods includes the following steps: Step (1) Initialize the displacement u i+1 and set the error threshold eps, the iteration increment h, and the maximum number of loop iterations cou; Step (2) Calculate the total force F received by the whole machine system; Step (3) calculates the effective stiffness of the system using the Newmark-β method and the effective load The expressions are as follows: Step (4) Solve the displacement response And define F1, the expression is as follows: Step (5) Increment each item in it by the iteration increment h, calculate the new non-linear acting force, and then calculate the differential equation of motion of the system to obtain Define F i , the i-th column of the Jacobian matrix is YKB(i), and the expression is as follows: YKB(i) = (F i - F1) / h Step (6) Use the Newton-Raphson method to calculate the new response numerical solution of the system, and the expression is as follows: Comparison | YKB -1 ×F1| Compare with the set error threshold eps: If the former is greater than the latter, continue the iterative calculation until it is less than the latter; when the former is less than the latter, the calculation result is taken as t n+1 The numerical solution of the response at this time enters step (7); If, when iterating to the maximum number of loop iterations cou, |YKB -1 ×F1| is still greater than the set error threshold eps, it is determined that the iteration does not converge, mark the calculation as incomplete, and proceed to step (8); Step (7) takes the numerical solution of the response at the time of the t n+1 as u i+1 , and calculates the velocity and acceleration responses. The expressions are as follows: Mark that the calculation is completed at the same time and enter step (8); Step (8) If the calculation is not completed, save the result file at this rotation speed of the rotor and return to step (1); if the calculation is completed, end the loop.

10. A computer device, comprising a memory, a processor, and a computer program stored on the memory, characterized in that, The processor executes the computer program to implement the steps of the method described in claims 1 to 9.

Images