CN113434983A - Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system - Google Patents

Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system Download PDF

Info

Publication number
CN113434983A
CN113434983A CN202110769503.5A CN202110769503A CN113434983A CN 113434983 A CN113434983 A CN 113434983A CN 202110769503 A CN202110769503 A CN 202110769503A CN 113434983 A CN113434983 A CN 113434983A
Authority
CN
China
Prior art keywords
bearing
oil film
sliding bearing
nonlinear
formula
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110769503.5A
Other languages
Chinese (zh)
Other versions
CN113434983B (en
Inventor
裴世源
张庆山
洪军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN202110769503.5A priority Critical patent/CN113434983B/en
Publication of CN113434983A publication Critical patent/CN113434983A/en
Application granted granted Critical
Publication of CN113434983B publication Critical patent/CN113434983B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Sliding-Contact Bearings (AREA)

Abstract

The invention discloses a method for rapidly calculating the nonlinear dynamic characteristics of a sliding bearing rotor system, which solves the problems of low efficiency and limited applicability range of the traditional method for calculating the nonlinear oil film force of a bearing under the special working conditions of explosion, impact, swing, rapid acceleration/deceleration and the like of the sliding bearing rotor system; according to the method, the rigidity damping coefficient of the node in the feasible region of the axis locus is fitted into a function at any position in the feasible region, and when the nonlinear dynamics of the rotor is subsequently solved, the nonlinear oil film force vector of the system at the t moment under a given working condition can be quickly obtained, so that the complex processes of repeated iteration and calling are avoided, and the efficiency of the nonlinear dynamics calculation process is higher; the calculation method provided by the invention comprises two parts of nonlinear oil film force calculation and rotor nonlinear dynamics calculation, each parameter index of the system is considered in detail, and a comprehensive and standard rapid calculation method is provided for the research of the nonlinear dynamics characteristics of the sliding bearing rotor system.

Description

Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system
Technical Field
The invention belongs to the field of rotor dynamics, and particularly relates to a method for rapidly calculating nonlinear dynamic characteristics of a sliding bearing rotor system.
Background
The sliding bearing rotor system is an important component of rotary mechanical equipment and is widely applied to the fields of aerospace, engineering machinery, national defense and military industry and the like. With continuous progress of science and technology and special requirements of markets, the working environment of rotary mechanical equipment gradually tends to high speed, high temperature and heavy load, so that a sliding bearing rotor system always works under special working conditions such as explosion, impact, swing, rapid acceleration/deceleration and the like, the load acting on a bearing is continuously changed along with time, the balance position of the axis of the rotor and the thickness of a bearing oil film are also changed at each moment, and the bearing can generate nonlinear oil film force, so that the sliding bearing rotor system shows very remarkable nonlinear characteristics. Relevant researches find that the nonlinear characteristics can greatly affect the reliability and safety of rotating mechanical equipment, so that the accurate, effective and rapid analysis of the nonlinear dynamic behavior of the sliding bearing rotor system is more and more important.
In a sliding bearing rotor system, the nonlinear oil film force of a bearing is one of the main factors of the system showing the nonlinear characteristic, so the key for obtaining the nonlinear characteristic of the system is the calculation of the nonlinear oil film force of the bearing. However, the traditional rigidity and damping calculation method is only suitable for the condition that the axle center track of the bearing changes in a small range, the calculation time is too long, and the calculation efficiency is low. The migration rate method and the impedance method are only suitable for the circularly symmetric sliding bearing, and have certain limitations on non-circularly symmetric sliding bearings such as oval, three-oil-leaf and step surfaces. Therefore, in order to solve the nonlinear dynamic characteristics of the sliding bearing rotor system under special working conditions and solve the problems of low calculation efficiency and limited application range of the traditional method, a calculation method for quickly solving the nonlinear dynamic characteristics of the sliding bearing rotor system is urgently needed.
Disclosure of Invention
The invention aims to provide a method for rapidly calculating nonlinear dynamic characteristics of a sliding bearing rotor system aiming at the problem that the nonlinear oil film force calculation time of the sliding bearing rotor system is too long when the sliding bearing rotor system works under special working conditions such as explosion, impact, swing, sudden acceleration/deceleration and the likeVector Fg, unbalance force vector F1And F2Axle center track feasible region S, bearing rigidity damping coefficient and nonlinear oil film force vector FoilAnd (4) calculating. Known conditions include rotor parameters, disc parameters, sliding bearing parameters, material parameters, and lubrication parameters.
The invention is realized by adopting the following technical scheme:
a method for rapidly calculating the nonlinear dynamic characteristics of a sliding bearing rotor system comprises the following steps:
1) establishing a finite element model and a motion differential equation of an Euler-Bernoulli beam of a sliding bearing rotor system according to given rotor parameters, disc parameters, sliding bearing parameters and material parameters, wherein each node of the finite element model of the Euler-Bernoulli beam of the sliding bearing rotor system has four degrees of freedom;
2) calculating a total mass matrix M, a total rigidity matrix K, a total gyro matrix G, a total gravity vector Fg and a total unbalance force vector F according to the finite element model and the motion differential equation in the step 1)1And F2；
3) Calculating the oil film thickness h (theta) of the bearing according to given sliding bearing parameters, and enabling h (theta) to be equal to zero to solve the axis track feasible region S of the bearing;
4) performing mesh division on the feasible region S by adopting a Delaunay triangulation algorithm, and refining the mesh at the position close to the boundary;
5) solving a Reynolds equation by using a finite difference method, and calculating the oil film pressure distribution p of the sliding bearing;
6) calculating the oil film bearing capacity F by using the oil film pressure distribution p obtained in the step 5)x、Fy；
7) Calculating the rigidity coefficient of the sliding bearing at each node in the feasible region S by using the oil film bearing capacity obtained in the step 6) and adopting a small disturbance method
Figure BDA0003152298480000021
Damping coefficient
Figure BDA0003152298480000022
i is a node in the feasible region SThe number of (2);
8) fitting the rigidity and damping coefficient calculated in the step 7) into a function related to the coordinates in the feasible region S by adopting a two-dimensional curved surface random interpolation function grirdata, and further obtaining a rigidity coefficient k of a sliding bearing at any position in the feasible region S of the axis trackxx、kxy、kyx、kyyDamping coefficient cxx、cxy、cyx、cyy；
9) Determining the eccentric position (x) of the bearing at the initial moment based on the eccentricity parametert0,yt0) And velocity
Figure BDA0003152298480000031
Determining iteration times and iteration termination time t according to given working conditionsmaxLet t be t 0;
10) determining a Delaunay triangular area according to the eccentric position of the bearing at the time t, and determining a node Q which is closest to the eccentric position in the Delaunay triangular area by using a Voronoi diagram;
11) invoking the stiffness coefficient k obtained by fitting in the step 8)xx、kxy、kyx、kyyAnd damping coefficient cxx、cxy、cyx、cyyAnd 10) solving the oil film force F of the bearing at the t moment according to the eccentric position and the node Q of the bearing at the t momentoilx、FoilyAnd assembling to obtain a nonlinear oil film force vector Foil；
12) Calling M, K, G, Fg and F obtained by calculation in the step 2)1、F2And the nonlinear oil film force vector F calculated in the step 11)oilSubstituting the system motion differential equation into the system motion differential equation established in the step 1), and solving the transient dynamic response of the system at the time t under the given working condition;
13) assuming that the time interval is delta t, calculating the eccentric position (x) of the bearing at the moment t + delta t according to the initial eccentric position, the external excitation load and the oil film pressure distribution in the step 5)t+Δt,yt+Δt) And velocity
Figure BDA0003152298480000032
14) Repeating the steps 10) to 13) and judging whether t is more than or equal to tmaxIf t < tmaxThen continue iteration; if t ≧ tmaxAnd then, finishing iteration to obtain the axis track of the bearing and the transient dynamic response of the system at each moment, and further analyzing the nonlinear dynamic characteristics of the system.
The invention is further improved in that, in step 1), the motion differential equation is as follows:
Figure BDA0003152298480000033
in the formula: m is a total mass matrix, G is a total gyro matrix, K is a total stiffness matrix, and Fg is a total weight vector; f1And F2The total unbalance force vector of the system in the x and y directions, FoilIs the sliding bearing nonlinear oil film force vector, q (t) is the generalized nodal displacement, and omega is the rotor speed.
The invention further improves that in the step 2), the total mass matrix M, the total rigidity matrix K, the total gyro matrix G, the total gravity vector Fg and the total unbalance force vector F1And F2Respectively assembled by each unit matrix, wherein the unit quality matrix MeCell stiffness matrix KeUnit gyro matrix GeCan be directly obtained according to the rotor dynamics theory, and the unit force vector Fge、F1 eAnd F2 eThe expression is as follows:
unit gravity vector Fge：
Fg e＝[0,Mig,0,0]T (2)
In the formula: miIs the mass of the rotor shaft section at the ith beam unit;
cell imbalance force vector F1 eAnd F2 e：
Figure BDA0003152298480000041
In the formula: m ised is the mass of the unbalance,
Figure BDA0003152298480000042
is the phase angle of the imbalance, omega is the rotor speed, F1 eRepresenting the x-direction cell imbalance force vector, F2 eRepresenting a y-direction cell imbalance force vector;
the total force vector Fg and the total unbalance force vector F obtained by assembly1And F2The expression is as follows:
total weight vector Fg:
Fg＝[0,M1g,0,0,……0,Mig,0,0,……,0,MN-1g,0,0]T (4)
in the formula: miThe mass of the rotor shaft section at the ith beam unit is shown, and N is the total number of nodes in the finite element model;
total unbalance force vector F1And F2：
Figure BDA0003152298480000043
In the formula: m ised is the mass of the unbalance,
Figure BDA0003152298480000044
is the unbalance phase angle, Ω is the rotor speed, j is the node number at the unbalance mass position, F1Representing the x-direction imbalance force vector, F2Representing the y-direction imbalance force vector.
The further improvement of the invention is that in the step 3), the expression of the oil film thickness of the fixed bush bearing is shown as a formula (6); for the tilting pad bearing, the expression of the oil film thickness is shown as the formula (7):
h(θ)＝cp-Xjcos(θ)-Yjsin(θ)-cpmcos(θ-θp) (6)
h(θ)＝cp-Xjcos(θ)-Yjsin(θ)-cpmcos(θ-θp)-(R+t)δsin(θ-θp) (7)
in the formula: theta is the circumferential position angle, cpIs a radial clearance of the plain bearing, XjIs the journal circumferential displacement, YjIs the axial displacement of the shaft neck, R is the radius of the bearing, m is the preload coefficient, and t is the pre-eccentricity.
The invention is further improved in that the Reynolds equation in the step 5) is as follows:
Figure BDA0003152298480000051
in the formula: gxIs a circumferential turbulence correction factor, GyIs an axial turbulence correction factor, eta is the viscosity of the lubricating oil, rho is the density of the lubricating oil, p is the pressure of an oil film, and x and y are circumferential and axial coordinates respectively;
the boundary conditions are as follows:
Figure BDA0003152298480000052
in the formula: psFor supply pressure, B is the bearing width.
The further improvement of the invention is that in the step 6), the calculation formula of the oil film bearing capacity of the sliding bearing in the x and y directions is as follows;
Figure BDA0003152298480000053
in the formula: theta is the circumferential position angle theta1Initial angle of oil film, theta2Oil film break angle.
The invention further improves the method that in the step 7), the journal generates slight disturbance in the x and y directions under the condition that the bearing axis is located at the ith node in a feasible region, and the axis position of the journal is moved from (x, y) to (x, y)0，y0) Then, then
Figure BDA0003152298480000054
Figure BDA0003152298480000055
The expression between is:
Figure BDA0003152298480000056
in the formula:
Figure BDA0003152298480000057
is the oil film bearing capacity of the axle center position of the shaft neck at (x, y),
Figure BDA0003152298480000058
the axle center position of the shaft neck is located at (x)0，y0) Oil film bearing capacity of the oil film;
the stiffness damping coefficient at the ith node in the feasible region is as follows:
Figure BDA0003152298480000061
in the formula:
Figure BDA0003152298480000062
respectively represents the direct rigidity of the ith node in the x direction and the y direction,
Figure BDA0003152298480000063
representing the cross stiffness at the ith node;
Figure BDA0003152298480000064
respectively represents the direct damping of the x and y directions at the ith node,
Figure BDA0003152298480000065
indicating the cross-damping at the ith node.
The further improvement of the invention is that in the step 8), the function obtained by fitting the two-dimensional curved surface random interpolation function grirdata is as follows:
Figure BDA0003152298480000066
in the formula: (x, y) are the coordinates of any point in the feasible region S, respectively, and the fitting functions all satisfy the data points at the nodes in step 7).
The invention is further improved in that in step 11), the oil film force F of the bearing at the time toilx、FoilyComprises the following steps:
Figure BDA0003152298480000067
in the formula: (x)t，yt) Is the eccentric position of the bearing at time t,
Figure BDA0003152298480000068
the velocities in the x and y directions at the eccentric position at time t, kxx、kxy、kyx、kyyIs (x)t，yt) Stiffness coefficient at location, cxx、cxy、cyx、cyyIs (x)t，yt) Damping coefficient at position, FQx、FQyIs the oil film bearing capacity, Δ x, at node Qt、ΔytThe distances between the eccentric position and the node Q in the x and y directions, respectively;
assembled nonlinear oil film force vector FoilThe expression is as follows:
Figure BDA0003152298480000069
in the formula: n is the total number of nodes in the finite element model, and p and q are the node numbers of the positions of the left sliding bearing and the right sliding bearing in the finite element model respectively.
A further improvement of the present invention is that, in step 13), the expressions of the eccentric positions of the bearing at the time t and at the time t + Δ t are as follows:
Figure BDA0003152298480000071
in the formula: (x)t，yt) Is the eccentric position of the bearing at time t, (x)t+Δt，yt+Δt) Is the eccentric position of the bearing at time t + at, at is the time interval,
Figure BDA0003152298480000072
the speed in x and y directions at the eccentric position at the time t respectively,
Figure BDA0003152298480000073
and
Figure BDA0003152298480000074
can be determined from the fact that the pressure profile p is equal to the external excitation load.
The invention has at least the following beneficial technical effects:
1) the method for rapidly calculating the nonlinear dynamic characteristics of the sliding bearing rotor system solves the problems of low calculation efficiency and limited application range of the nonlinear oil film force of the bearing under the special working conditions of explosion, impact, swing, sudden acceleration/deceleration and the like solved by the traditional calculation method.
2) According to the method, the rigidity and the damping coefficient of the node in the feasible region of the axis locus are fitted into the two-dimensional surface function at any position in the feasible region, and when the nonlinear dynamics of the rotor is subsequently solved, the fitting result, the eccentric position, the speed and the Voronoi diagram are combined to quickly obtain the nonlinear oil film force vector of the system at the t moment under the given working condition, so that the complex processes of repeated iteration and calling are avoided, the nonlinear dynamics calculation process is simpler, and the efficiency is higher.
3) The calculation method provided by the invention comprises two parts of nonlinear oil film force calculation and rotor nonlinear dynamics calculation, and specifically comprises a system mass matrix M, a rigidity matrix K, a gyro matrix G, a gravity vector Fg and an unbalanced force vector F1And F2Axle center track feasible region S, bearing rigidity damping coefficient and nonlinear oil film force vectorQuantity FoilThe calculation of the method takes various parameter indexes of the system into consideration in detail, and provides a comprehensive and standard rapid calculation method for the research of the nonlinear dynamic characteristics of the sliding bearing rotor system.
Drawings
Fig. 1 is a flow chart of a method for rapidly calculating the nonlinear dynamic characteristics of a sliding bearing rotor system.
FIG. 2 is a Euler-Bernoulli beam finite element model of a sliding bearing rotor system.
FIG. 3 is a schematic diagram of matrix assembly.
Fig. 4 is a schematic diagram of the feasible region meshing of the axis locus.
FIG. 5 is a schematic diagram of an axial trace.
Detailed Description
The invention is described in further detail below with reference to the attached drawing figures:
as shown in fig. 1, which is a flowchart of a fast calculation method, the fast calculation method for nonlinear dynamics characteristics of a sliding bearing rotor system provided by the present invention includes the following steps:
1) and establishing a finite element model (4 degrees of freedom of each node) and a motion differential equation of the Euler-Bernoulli beam according to the given rotor parameters, disc parameters and material parameters. The finite element model of the Euler-Bernoulli beam is shown in FIG. 2, and the motion differential equation is as follows:
Figure BDA0003152298480000081
in the formula: m is a total mass matrix, G is a total gyro matrix, K is a total stiffness matrix, and Fg is a total weight vector; f1And F2The total unbalance force vector of the system in the x and y directions, FoilIs the nonlinear oil film force vector of the sliding bearing, q (t) is the generalized node displacement, and omega is the rotor speed;
2) calculating a total mass matrix M, a total rigidity matrix K, a total gyro matrix G, a total gravity vector Fg and a total unbalance force vector F according to the finite element model and the motion differential equation in the step 1)1And F2. Each overall matrix is assembled by each unit matrix, wherein the unit quality matrix MeCell stiffness matrix KeUnit gyro matrix GeCan be directly obtained according to the rotor dynamics theory, and the unit force vector Fge、F1 eAnd F2 eThe expression is as follows:
unit gravity vector Fge：
Fg e＝[0,Mig,0,0]T (2)
In the formula: miIs the mass of the rotor shaft section at the ith beam element.
Cell imbalance force vector F1 eAnd F2 e：
Figure BDA0003152298480000091
In the formula: m ised is the mass of the unbalance,
Figure BDA0003152298480000092
is the phase angle of the imbalance, omega is the rotor speed, F1 eRepresenting the x-direction cell imbalance force vector, F2 eRepresenting the y-direction cell imbalance force vector.
The overall mass matrix M, the overall stiffness matrix K and the overall gyro matrix G are assembled in the manner shown in FIG. 3, and the overall weight vector Fg and the overall unbalance vector F1And F2The assembly is as follows:
total weight vector Fg:
Fg＝[0,M1g,0,0,……0,Mig,0,0,……,0,MN-1g,0,0]T (4)
in the formula: miIs the mass of the rotor shaft section at the ith beam element, and N is the total number of nodes in the finite element model.
Total unbalance force vector F1And F2：
Figure BDA0003152298480000093
In the formula: m ised is the mass of the unbalance,
Figure BDA0003152298480000094
is the unbalance phase angle, Ω is the rotor speed, j is the node number at the unbalance mass position, F1Representing the x-direction imbalance force vector, F2Representing the y-direction imbalance force vector.
3) And calculating the oil film thickness h (theta) of the bearing according to given sliding bearing parameters, and enabling h (theta) to be equal to zero to solve the feasible region S of the axis track of the bearing. For the fixed bush bearing, the expression of the oil film thickness is shown as the formula (6); for the tilting pad bearing, the expression of the oil film thickness is shown as formula (7);
h(θ)＝cp-Xjcos(θ)-Yjsin(θ)-cpmcos(θ-θp) (6)
h(θ)＝cp-Xjcos(θ)-Yjsin(θ)-cpmcos(θ-θp)-(R+t)δsin(θ-θp) (7)
in the formula: theta is the circumferential position angle, cpIs a radial clearance of the plain bearing, XjIs the journal circumferential displacement, YjIs the axial displacement of the shaft neck, R is the radius of the bearing, m is the preload coefficient, and t is the pre-eccentricity.
4) And adopting a Delaunay triangulation algorithm to perform meshing on the feasible region. In order to improve the accuracy of the subsequent calculation result, the grid needs to be refined near the boundary, as shown in fig. 4;
5) solving a Reynolds equation by using a finite difference method, and calculating the oil film pressure distribution p of the sliding bearing, wherein the expression of the Reynolds equation is as follows:
Figure BDA0003152298480000101
in the formula: gxIs a circumferential turbulence correction factor, GyIs axial turbulence repairPositive factors, η is the viscosity of the lubricating oil, ρ is the density of the lubricating oil, p is the oil film pressure, and x and y are the circumferential and axial coordinates, respectively.
The boundary conditions are as follows:
Figure BDA0003152298480000102
in the formula: psFor supply pressure, B is the bearing width.
6) Calculating the oil film bearing capacity F by using the oil film pressure distribution p obtained in the step 5)x、FyThe specific calculation formula is as follows:
Figure BDA0003152298480000103
in the formula: theta is the circumferential position angle theta1Initial angle of oil film, theta2Oil film break angle
7) Calculating the rigidity coefficient of the sliding bearing at each node in the feasible region by using the oil film bearing capacity obtained in the step 6) and adopting a small disturbance method
Figure BDA0003152298480000104
Damping coefficient
Figure BDA0003152298480000105
i is the number of nodes within the feasible region S. Assuming that the bearing axis is located at the ith node in the feasible region, the journal generates small disturbance in the x and y directions, and the axis position of the journal moves from (x, y) to (x0，y0) Then, then
Figure BDA0003152298480000106
The expression between is:
Figure BDA0003152298480000111
in the formula:
Figure BDA0003152298480000112
is the oil film bearing capacity of the axle center position of the shaft neck at (x, y),
Figure BDA0003152298480000113
the axle center position of the shaft neck is located at (x)0，y0) The oil film bearing capacity of (c).
The stiffness damping coefficient at the ith node in the feasible region is as follows:
Figure BDA0003152298480000114
in the formula:
Figure BDA0003152298480000115
respectively represents the direct rigidity of the ith node in the x direction and the y direction,
Figure BDA0003152298480000116
representing the cross stiffness at the ith node;
Figure BDA0003152298480000117
respectively represents the direct damping of the x and y directions at the ith node,
Figure BDA0003152298480000118
indicating the cross-damping at the ith node.
8) Fitting the rigidity and damping coefficient of the node calculated in the step 7) into a function related to any coordinate in a feasible domain by using a two-dimensional curved surface random interpolation function griddata, and further obtaining a rigidity coefficient k of a sliding bearing at any position in the feasible domain of the axis trackxx、kxy、kyx、kyyDamping coefficient cxx、cxy、cyx、cyy. The function resulting from the above fitting is as follows:
Figure BDA0003152298480000119
in the formula: (x, y) are the coordinates of any point in the feasible region S, respectively, and the fitting functions all satisfy the data points at the nodes in step 7).
9) Determining the eccentric position (x) of the bearing at the initial moment based on the eccentricity parametert0,yt0) And velocity
Figure BDA00031522984800001110
Determining iteration times and iteration termination time t according to given working conditionsmaxLet t be t 0;
10) and determining a Delaunay triangular area according to the eccentric position of the bearing at the time t, and determining a node Q which is closest to the eccentric position in the Delaunay triangular area by using a Voronoi diagram.
11) Invoking the stiffness coefficient k obtained by fitting in the step 8)xx、kxy、kyx、kyyAnd damping coefficient cxx、cxy、cyx、cyyAnd 10) solving the oil film force F of the bearing at the t moment according to the eccentric position and the node Q of the bearing at the t momentoilx、FoilyAnd assembling to obtain a nonlinear oil film force vector Foil. Oil film force F of bearing at time toilx、FoilyThe calculation formula is as follows:
Figure BDA0003152298480000121
in the formula: (x)t，yt) Is the eccentric position of the bearing at time t,
Figure BDA0003152298480000122
the velocities in the x and y directions at the eccentric position at time t, kxx、kxy、kyx、kyyIs (x)t，yt) Stiffness coefficient at location, cxx、cxy、cyx、cyyIs (x)t，yt) Damping coefficient at position, FQx、FQyOil film bearing capacity, Δ x, at node Qt、ΔytAre respectively inclined toThe distance between the heart position and the node Q in the x and y directions.
Assembled nonlinear oil film force vector FoilThe expression is as follows:
Figure BDA0003152298480000123
in the formula: n is the total number of nodes in the finite element model, and p and q are the node numbers of the positions of the left sliding bearing and the right sliding bearing in the finite element model respectively.
12) Calling M, K, G, Fg and F obtained by calculation in the step 2)1、F2And 11) calculating the obtained nonlinear oil film force vector FoilSubstituting the system motion differential equation into the system motion differential equation established in the step 1), and solving the transient dynamic response of the system at the time t under the given working condition;
13) assuming that the time interval is delta t, calculating the eccentric position (x) of the bearing at the moment t + delta t according to the initial eccentric position, the external excitation load and the oil film pressure distribution in the step 3)t+Δt,yt+Δt) And velocity
Figure BDA0003152298480000124
the expression of the eccentric position of the bearing at the time t and the time t + deltat is as follows:
Figure BDA0003152298480000125
in the formula: (x)t，yt) Is the eccentric position of the bearing at time t, (x)t+Δt，yt+Δt) Is the eccentric position of the bearing at time t + at, at is the time interval,
Figure BDA0003152298480000126
the speed in x and y directions at the eccentric position at the time t respectively,
Figure BDA0003152298480000127
and
Figure BDA0003152298480000128
can be determined from the fact that the pressure profile p is equal to the external excitation load.
14) Repeating the steps 10) to 13) and judging whether t is more than or equal to tmax. If t < tmaxThen continue iteration; if t ≧ tmaxThen the iteration ends. Thus, the axial locus and the transient dynamic response of the system at each moment as shown in fig. 5 can be obtained, and the nonlinear dynamic characteristics of the system can be further analyzed.
In conclusion, the invention provides a method for rapidly calculating the nonlinear dynamic characteristic of a sliding bearing rotor system, the feasible region of the axle center track of the bearing is subdivided by a Delaunay triangulation algorithm, the rigidity and the damping coefficient of the nodes in the feasible region of the axle center track are fitted into a function of any coordinate in the feasible region by adopting a two-dimensional curved surface random interpolation function griddata, when the nonlinear dynamics of the rotor is subsequently solved, the fitting result, the eccentric position, the speed and the Voronoi diagram are combined to quickly obtain the nonlinear oil film force vector of the system at the t moment under a given working condition, so that the complex processes of repeated iteration and calling are avoided, the defects of low calculation efficiency, poor applicability and the like of the traditional method are overcome, and an effective method is provided for quickly calculating the nonlinear dynamics characteristics of the sliding bearing rotor system under special working conditions such as explosion, impact, swing, sudden acceleration/deceleration and the like.

Claims (10)

1. A method for rapidly calculating the nonlinear dynamic characteristics of a sliding bearing rotor system is characterized by comprising the following steps:
1) establishing a finite element model and a motion differential equation of an Euler-Bernoulli beam of a sliding bearing rotor system according to given rotor parameters, disc parameters, sliding bearing parameters and material parameters, wherein each node of the finite element model of the Euler-Bernoulli beam of the sliding bearing rotor system has four degrees of freedom;
2) calculating a total mass matrix M, a total rigidity matrix K, a total gyro matrix G, a total gravity vector Fg and a total unbalance force vector F according to the finite element model and the motion differential equation in the step 1)1And F2；
3) Calculating the oil film thickness h (theta) of the bearing according to given sliding bearing parameters, and enabling h (theta) to be equal to zero to solve the axis track feasible region S of the bearing;
4) performing mesh division on the feasible region S by adopting a Delaunay triangulation algorithm, and refining the mesh at the position close to the boundary;
5) solving a Reynolds equation by using a finite difference method, and calculating the oil film pressure distribution p of the sliding bearing;
6) calculating the oil film bearing capacity F by using the oil film pressure distribution p obtained in the step 5)x、Fy；
7) Calculating the rigidity coefficient of the sliding bearing at each node in the feasible region S by using the oil film bearing capacity obtained in the step 6) and adopting a small disturbance method
Figure FDA0003152298470000011
Damping coefficient
Figure FDA0003152298470000012
i is the serial number of the nodes in the feasible region S;
8) fitting the rigidity and damping coefficient calculated in the step 7) into a function related to the coordinates in the feasible region S by adopting a two-dimensional curved surface random interpolation function grirdata, and further obtaining a rigidity coefficient k of a sliding bearing at any position in the feasible region S of the axis trackxx、kxy、kyx、kyyDamping coefficient cxx、cxy、cyx、cyy；
9) Determining the eccentric position (x) of the bearing at the initial moment based on the eccentricity parametert0,yt0) And velocity
Figure FDA0003152298470000013
Determining iteration times and iteration termination time t according to given working conditionsmaxLet t be t 0;
10) determining a Delaunay triangular area according to the eccentric position of the bearing at the time t, and determining a node Q which is closest to the eccentric position in the Delaunay triangular area by using a Voronoi diagram;
11) invoking the stiffness coefficient k obtained by fitting in the step 8)xx、kxy、kyx、kyyAnd damping coefficient cxx、cxy、cyx、cyyAnd 10) solving the oil film force F of the bearing at the t moment according to the eccentric position and the node Q of the bearing at the t momentoilx、FoilyAnd assembling to obtain a nonlinear oil film force vector Foil；
12) Calling M, K, G, Fg and F obtained by calculation in the step 2)1、F2And the nonlinear oil film force vector F calculated in the step 11)oilSubstituting the system motion differential equation into the system motion differential equation established in the step 1), and solving the transient dynamic response of the system at the time t under the given working condition;
13) assuming that the time interval is delta t, calculating the eccentric position (x) of the bearing at the moment t + delta t according to the initial eccentric position, the external excitation load and the oil film pressure distribution in the step 5)t+Δt,yt+Δt) And velocity
Figure FDA0003152298470000021
14) Repeating the steps 10) to 13) and judging whether t is more than or equal to tmaxIf t < tmaxThen continue iteration; if t ≧ tmaxAnd then, finishing iteration to obtain the axis track of the bearing and the transient dynamic response of the system at each moment, and further analyzing the nonlinear dynamic characteristics of the system.
2. The sliding bearing rotor system nonlinear dynamics fast calculation method according to claim 1, characterized in that in step 1), the motion differential equation is:
Figure FDA0003152298470000022
in the formula: m is a total mass matrix, G is a total gyro matrix, K is a total stiffness matrix, and Fg is a total weight vector; f1And F2Are respectively asOverall unbalance force vector of system in x, y direction, FoilIs the sliding bearing nonlinear oil film force vector, q (t) is the generalized nodal displacement, and omega is the rotor speed.
3. The method for rapidly calculating the nonlinear dynamics of the sliding bearing rotor system according to claim 2, wherein in the step 2), the total mass matrix M, the total stiffness matrix K, the total gyro matrix G, the total weight vector Fg and the total unbalance force vector F1And F2Respectively assembled by each unit matrix, wherein the unit quality matrix MeCell stiffness matrix KeUnit gyro matrix GeCan be directly obtained according to the rotor dynamics theory, and the unit force vector Fge、F1 eAnd F2 eThe expression is as follows:
unit gravity vector Fge：
Figure FDA0003152298470000023
In the formula: miIs the mass of the rotor shaft section at the ith beam unit;
cell imbalance force vector F1 eAnd F2 e：
Figure FDA0003152298470000031
In the formula: m ised is the mass of the unbalance,
Figure FDA0003152298470000032
is the phase angle of the imbalance, omega is the rotor speed, F1 eRepresenting the x-direction cell imbalance force vector, F2 eRepresenting a y-direction cell imbalance force vector;
the total force vector Fg and the total unbalance force vector F obtained by assembly1And F2The expression is as followsShown in the figure:
total weight vector Fg:
Fg＝[0,M1g,0,0,……0,Mig,0,0,……,0,MN-1g,0,0]T (4)
in the formula: miThe mass of the rotor shaft section at the ith beam unit is shown, and N is the total number of nodes in the finite element model;
total unbalance force vector F1And F2：
Figure FDA0003152298470000033
In the formula: m ised is the mass of the unbalance,
Figure FDA0003152298470000034
is the unbalance phase angle, Ω is the rotor speed, j is the node number at the unbalance mass position, F1Representing the x-direction imbalance force vector, F2Representing the y-direction imbalance force vector.
4. The method for rapidly calculating the nonlinear dynamics characteristic of the sliding bearing rotor system according to the claim 3, characterized in that in the step 3), the expression of the oil film thickness of the fixed pad bearing is shown as the formula (6); for the tilting pad bearing, the expression of the oil film thickness is shown as the formula (7):
h(θ)＝cp-Xjcos(θ)-Yjsin(θ)-cpmcos(θ-θp) (6)
h(θ)＝cp-Xjcos(θ)-Yjsin(θ)-cpmcos(θ-θp)-(R+t)δsin(θ-θp) (7)
in the formula: theta is the circumferential position angle, cpIs a radial clearance of the plain bearing, XjIs the journal circumferential displacement, YjIs the axial displacement of the shaft neck, R is the radius of the bearing, m is the preload coefficient, and t is the pre-eccentricity.
5. The method for rapidly calculating the nonlinear dynamics characteristic of the sliding bearing rotor system according to claim 4, wherein the Reynolds equation in the step 5) is as follows:
Figure FDA0003152298470000041
in the formula: gxIs a circumferential turbulence correction factor, GyIs an axial turbulence correction factor, eta is the viscosity of the lubricating oil, rho is the density of the lubricating oil, p is the pressure of an oil film, and x and y are circumferential and axial coordinates respectively;
the boundary conditions are as follows:
Figure FDA0003152298470000042
in the formula: psFor supply pressure, B is the bearing width.
6. The method for rapidly calculating the nonlinear dynamics characteristic of the sliding bearing rotor system according to claim 5, characterized in that in the step 6), the calculation formula of the oil film bearing capacity of the sliding bearing in the x and y directions is as follows;
Figure FDA0003152298470000043
in the formula: theta is the circumferential position angle theta1Initial angle of oil film, theta2Oil film break angle.
7. The method for rapidly calculating the nonlinear dynamics characteristics of the sliding bearing rotor system according to claim 6, wherein in step 7), assuming that the bearing axis is located at the ith node in the feasible region, the journal generates small disturbances in the x and y directions, and the axis position of the journal moves from (x, y) to (x, y)0，y0) Then F isx i、Fy i、Fx0 i、Fy0 iThe expression between is:
Figure FDA0003152298470000044
in the formula: fx i、Fy iIs the oil film bearing capacity of the journal axis position at (x, y), Fx0 i、Fy0 iThe axle center position of the shaft neck is located at (x)0，y0) Oil film bearing capacity of the oil film;
the stiffness damping coefficient at the ith node in the feasible region is as follows:
Figure FDA0003152298470000051
in the formula:
Figure FDA0003152298470000052
respectively represents the direct rigidity of the ith node in the x direction and the y direction,
Figure FDA0003152298470000053
representing the cross stiffness at the ith node;
Figure FDA0003152298470000054
respectively represents the direct damping of the x and y directions at the ith node,
Figure FDA0003152298470000055
indicating the cross-damping at the ith node.
8. The method for rapidly calculating the nonlinear dynamics characteristic of the sliding bearing rotor system according to claim 7, wherein in the step 8), the function obtained by fitting the two-dimensional curved surface random interpolation function griddata is as follows:
Figure FDA0003152298470000056
in the formula: (x, y) are the coordinates of any point in the feasible region S, respectively, and the fitting functions all satisfy the data points at the nodes in step 7).
9. The method for rapidly calculating the nonlinear dynamics characteristic of the sliding bearing rotor system according to claim 8, wherein in the step 11), the oil film force F of the bearing at the time toilx、FoilyComprises the following steps:
Figure FDA0003152298470000057
in the formula: (x)t，yt) Is the eccentric position of the bearing at time t,
Figure FDA0003152298470000058
the velocities in the x and y directions at the eccentric position at time t, kxx、kxy、kyx、kyyIs (x)t，yt) Stiffness coefficient at location, cxx、cxy、cyx、cyyIs (x)t，yt) Damping coefficient at position, FQx、FQyIs the oil film bearing capacity, Δ x, at node Qt、ΔytThe distances between the eccentric position and the node Q in the x and y directions, respectively;
assembled nonlinear oil film force vector FoilThe expression is as follows:
Figure FDA0003152298470000059
in the formula: n is the total number of nodes in the finite element model, and p and q are the node numbers of the positions of the left sliding bearing and the right sliding bearing in the finite element model respectively.
10. The method for rapidly calculating the nonlinear dynamics characteristic of the sliding bearing rotor system according to claim 9, wherein in the step 13), the expressions of the eccentric positions of the bearings at the time t and at the time t + Δ t are as follows:
Figure FDA0003152298470000061
in the formula: (x)t，yt) Is the eccentric position of the bearing at time t, (x)t+Δt，yt+Δt) Is the eccentric position of the bearing at time t + at, at is the time interval,
Figure FDA0003152298470000062
the speed in x and y directions at the eccentric position at the time t respectively,
Figure FDA0003152298470000063
and
Figure FDA0003152298470000064
can be determined from the fact that the pressure profile p is equal to the external excitation load.
CN202110769503.5A 2021-07-07 2021-07-07 Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system Active CN113434983B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110769503.5A CN113434983B (en) 2021-07-07 2021-07-07 Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110769503.5A CN113434983B (en) 2021-07-07 2021-07-07 Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system

Publications (2)

Publication Number Publication Date
CN113434983A true CN113434983A (en) 2021-09-24
CN113434983B CN113434983B (en) 2022-12-09

Family

ID=77759537

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110769503.5A Active CN113434983B (en) 2021-07-07 2021-07-07 Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system

Country Status (1)

Country Link
CN (1) CN113434983B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114282315A (en) * 2021-11-29 2022-04-05 中国船舶工业集团公司第七0八研究所 Method for calculating critical rotating speed of water-lubricated bearing-multi-annular seal-rotor system
CN114580121A (en) * 2022-05-05 2022-06-03 西安航天动力研究所 Method, device and medium for calculating dynamic characteristics of rotor system based on finite element method

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0969396A2 (en) * 1998-07-01 2000-01-05 Tadahiko Kawai Numerical analyzing method, element analyzing blocks, and processing apparatus for numerical analysis
WO2011017071A2 (en) * 2009-07-28 2011-02-10 University Of Kansas Method and apparatus for pressure adaptive morphing structure
CN108956143A (en) * 2018-06-25 2018-12-07 西安理工大学 A kind of transversal crack fault characteristic value extracting method of rotor-bearing system
CN109829262A (en) * 2019-04-04 2019-05-31 哈尔滨工程大学 A kind of rotor-bearing system nonlinear dynamic analysis method
CN110096784A (en) * 2019-04-25 2019-08-06 西安交通大学 A kind of quick calculating and design method of the bush(ing) bearing with axial pressure difference
US20200401745A1 (en) * 2019-06-19 2020-12-24 National Central University Structure analyzing method, device, and non-transitory computer-readable medium

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0969396A2 (en) * 1998-07-01 2000-01-05 Tadahiko Kawai Numerical analyzing method, element analyzing blocks, and processing apparatus for numerical analysis
WO2011017071A2 (en) * 2009-07-28 2011-02-10 University Of Kansas Method and apparatus for pressure adaptive morphing structure
CN108956143A (en) * 2018-06-25 2018-12-07 西安理工大学 A kind of transversal crack fault characteristic value extracting method of rotor-bearing system
CN109829262A (en) * 2019-04-04 2019-05-31 哈尔滨工程大学 A kind of rotor-bearing system nonlinear dynamic analysis method
CN110096784A (en) * 2019-04-25 2019-08-06 西安交通大学 A kind of quick calculating and design method of the bush(ing) bearing with axial pressure difference
US20200401745A1 (en) * 2019-06-19 2020-12-24 National Central University Structure analyzing method, device, and non-transitory computer-readable medium

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
李强等: "非线性转子-轴承耦合系统润滑及稳定性分析", 《浙江大学学报(工学版)》 *
毛文贵等: "基于流固耦合的滑动轴承非线性油膜动特性研究", 《中国机械工程》 *
郑铁生，许庆余: "滑动轴承动特性系数新算法", 《机械工程学报》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114282315A (en) * 2021-11-29 2022-04-05 中国船舶工业集团公司第七0八研究所 Method for calculating critical rotating speed of water-lubricated bearing-multi-annular seal-rotor system
CN114580121A (en) * 2022-05-05 2022-06-03 西安航天动力研究所 Method, device and medium for calculating dynamic characteristics of rotor system based on finite element method
CN114580121B (en) * 2022-05-05 2022-08-16 西安航天动力研究所 Method, device and medium for calculating dynamic characteristics of rotor system based on finite element method
WO2023213017A1 (en) * 2022-05-05 2023-11-09 西安航天动力研究所 Rotor system dynamic characteristic calculation method and device based on finite element method, and medium

Also Published As

Publication number Publication date
CN113434983B (en) 2022-12-09

Similar Documents

Publication Publication Date Title
CN113434983B (en) Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system
CN109829262B (en) Nonlinear dynamics analysis method for rotor-bearing system
Bonneau et al. Finite element analysis of grooved gas thrust bearings and grooved gas face seals
Rashidi et al. Bifurcation and nonlinear dynamic analysis of a rigid rotor supported by two-lobe noncircular gas-lubricated journal bearing system
CN110096784B (en) Rapid calculation and design method of radial sliding bearing with axial pressure difference
Gwynllyw et al. On the effects of a piezoviscous lubricant on the dynamics of a journal bearing
Meybodi et al. Numerical analysis of a rigid rotor supported by aerodynamic four-lobe journal bearing system with mass unbalance
Rao et al. An analytical approach to evaluate dynamic coefficients and nonlinear transient analysis of a hydrodynamic journal bearing
Yu et al. Inclination angle effect of tribological performance for hydrostatic bearing having tilting oil pad under variable viscosity conditions
CN107526914B (en) Variable-watershed flow field calculation method of tilting-pad sliding bearing based on structured dynamic grid
Wang Bifurcation and nonlinear analysis of a flexible rotor supported by a relative short spherical gas bearing system
CN115935687B (en) Method for calculating coupling lubrication and dynamic characteristic parameters of flanging bearing
CN110378018B (en) Method for calculating steady-state bearing capacity of hydrodynamic and hydrostatic ball bearing
Fang et al. A comprehensive analysis of factors affecting the accuracy of the precision hydrostatic spindle with mid-thrust bearing layout
Yang et al. Calculation of the dynamic characteristics of ship’s aft stern tube bearing considering journal deflection
Martinez Esparza et al. Design of hybrid hydrostatic/hydrodynamic journal bearings for optimum self-compensation under misaligning external loads
Singh et al. Modeling and simulation of bearing clearance effects on journal center motion trajectories
CN114091314A (en) Vibration prediction method of rotor system model based on magneto-rheological damper
Holmes The effect of sleeve bearings on the vibration of rotating shafts
Yadav et al. Transient response of two lobe aerodynamic journal bearing
Li et al. Dynamic characteristics of opposed-conical gas-dynamic bearings
CN116090135B (en) Damper-rotor system response analysis method
Li et al. On the lubrication performance of journal bearing coupled with the axial movement of journal
CN116976018B (en) Method for predicting relative position deviation of milling tool system
Zadorozhnaya et al. Modelling the thermal state of a turbocharger bearing housing when calculating the rotor dynamics at transient modes

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant