CN113254863A - Self-adaptive orthogonal decomposition method for order reduction of multi-disk rotor system model - Google Patents

Self-adaptive orthogonal decomposition method for order reduction of multi-disk rotor system model Download PDF

Info

Publication number
CN113254863A
CN113254863A CN202110492921.4A CN202110492921A CN113254863A CN 113254863 A CN113254863 A CN 113254863A CN 202110492921 A CN202110492921 A CN 202110492921A CN 113254863 A CN113254863 A CN 113254863A
Authority
CN
China
Prior art keywords
matrix
order
equation
expressed
reduced
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
CN202110492921.4A
Other languages
Chinese (zh)
Other versions
CN113254863B (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.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical 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 Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202110492921.4A priority Critical patent/CN113254863B/en
Publication of CN113254863A publication Critical patent/CN113254863A/en
Application granted granted Critical
Publication of CN113254863B publication Critical patent/CN113254863B/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
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G01M13/045—Acoustic or vibration analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Acoustics & Sound (AREA)
  • Computing Systems (AREA)
  • Complex Calculations (AREA)

Abstract

The invention provides a self-adaptive orthogonal decomposition method for reducing orders of a multi-disc rotor system model, which is characterized in that solutions of system dynamics control equations under different modeling parameters are obtained in advance through a four-order Runge Kutta method, an autocorrelation matrix is obtained according to a snapshot matrix formed by the solutions, reduced orders POM under new parameters are obtained by calculating characteristic values of the autocorrelation matrix, and finally the reduced orders POM simplifies complex original system dynamics control equations under the new model parameters. The invention can obtain the optimal reduced order model of the complex rotor system, can keep the orthogonal normalization characteristic of the reduced order POM, eliminates the local characteristic of the ITGM method, is used in a wide frequency parameter range, and has accurate calculation, higher calculation efficiency and strong robustness.

Description

Self-adaptive orthogonal decomposition method for order reduction of multi-disk rotor system model
Technical Field
The invention relates to the field of dynamics, in particular to an orthogonal decomposition method for model order reduction.
Background
A rotor bearing system is one of the most important components in mechanical systems such as aircraft engines, gas turbines, large generators, and the like. Generally, the rotor bearing system comprises a plurality of rotating discs, which are very complex and have high dimension, and it is necessary to reduce the dimension of the rotor bearing system by using a proper method for reducing the dimension of the system.
The POD method obtains the most important components of an infinite dimensional continuous system or a finite high dimensional dofs (degrees of freedom) system through intrinsic orthogonal modes (POMs). The POD method can greatly reduce the degree of freedom of the complex system, and meanwhile, can improve the calculation efficiency of the complex system on the premise of keeping higher calculation accuracy, so that the POD method is widely applied to the related field of fluid mechanics. However, the Reduced Order Model (ROM) obtained by the POD method generally lacks robustness when system parameters are changed. Compared with the IGM method, the ITGM method needs more reduced-order modes for approaching the original system and has poorer approaching precision.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides an adaptive orthogonal decomposition method for reducing the order of a multi-disk rotor system model, namely an IGM method. The method comprises the steps of obtaining solutions of system dynamics control equations under different modeling parameters in advance through a four-order Runge Kutta method, obtaining an autocorrelation matrix according to a snapshot matrix formed by the solutions, further obtaining reduced POM under a new parameter by calculating a characteristic value of the autocorrelation matrix, and finally simplifying a complex original system dynamics control equation under the new model parameter through the reduced POM.
To solve the above problems in the POD method, the present invention proposes an adaptive POD method (IGM method) based on the grassmann manifold and the lagrange interpolation. Firstly, obtaining solutions of system dynamics control equations under different modeling parameters in advance through a fourth-order Runge Kutta method, then obtaining an autocorrelation matrix according to a snapshot matrix of the solutions, obtaining a reduced POM under a new parameter through the autocorrelation matrix, and finally simplifying a complex original system dynamics control equation under the new model parameter through the reduced POM. By analyzing the response result of the model, the IGM method not only obtains the optimal ROM (Reduced-Order Models) of a complex physical system, maintains the orthogonal normalization characteristic of the Reduced-Order POM, but also eliminates the limitation of the ITGM method. The IGM method is accurate in calculation in a wider parameter domain, high in calculation efficiency and strong in robustness.
The technical scheme adopted by the invention for solving the technical problem comprises the following steps:
step one, obtaining a differential equation of motion of a high-dimensional nonlinear rotor system;
the rotor bearing system is a multi-disc rotor bearing system, is supported by two oil lubrication short shaft journal bearings at two ends, has a base loosening fault at a left support position, and is provided with 14 discs on a shaft;
the shafts are elastic, and the mass of each section of shaft is concentrated on the mass center of the corresponding disc, and the equivalent rigidity of the shaft between the turntables is k without considering the gyro moment of the disciI 1.. 15, the horizontal and vertical displacements of the geometric center of the disk and the center of mass of the support end are denoted Xi,YiI 1.. 16, the nonlinear oil film force in the horizontal and vertical directions is denoted as Fx,FyVibration with loosening failure only in the vertical direction, corresponding to a displacement of Y17Base klLoose stiffness and base clIs a piecewise linear function expressed as:
Figure BDA0003053104370000021
Figure BDA0003053104370000022
the differential equation of motion for a rotor bearing system is expressed as:
Figure BDA0003053104370000023
wherein X ═ X1,...X16,Y1...Y16,Y17]M, C, K are the mass, damping and stiffness matrices, respectively, and F is the resultant of the nonlinear fluid film force, the imbalance force and the gravity of each disk, expressed as:
Figure BDA0003053104370000024
M＝diag(m1,…,m16,m1,…m16,m17)
C＝diag(c1,…,c16,c1,…c16,c17)
Figure BDA0003053104370000031
carrying out dimensionless transformation on a dynamic control equation of the system, wherein the formula (3) is dimensionless by the following formula:
Figure BDA0003053104370000032
wherein f isxAnd fyIs a dimensionless non-linear oil film force,
Figure BDA0003053104370000033
is the Sommerfeld number, μ is the lubricant viscosity, c is the oil film thickness, L is the bearing length, R is the bearing radius, ω is the angular frequency, τ is the dimensionless time, W is half the rotor system weight;
thereby establishing a high-dimensional nonlinear rotor system model;
inputting the response signal into a snapshot matrix and calculating an autocorrelation matrix;
manifold G (M, N)s) Is provided with a plurality of sampling time sequences NsThe M-dimensional column vector coverage of the corresponding solution;
n precomputed snapshot signal sets of dimension M system are set as
Figure BDA0003053104370000034
ωi∈[ωα,ωβ]Wherein each parameter ωiSample length t ofsSame, snapshot set
Figure BDA0003053104370000035
Viewed as manifold G (M, N)s) N points of
Figure BDA0003053104370000036
Selecting a point S of the manifold0As a reference point, obtaining a snapshot matrix y with a corresponding parameter omega on a manifold by a Lagrange interpolation method*(ω,ts) Expressed as:
Figure BDA0003053104370000037
manifold G (M, N)s) The upper interpolation point is represented as:
y(ωi,ts)＝UiΣiVi t,i＝0,...,N-1 (7)
wherein, sigmai∈RM×MIs a matrix of singular values and is,
Figure BDA0003053104370000041
Vi∈RM×Mis an orthogonal matrix, then equation (7) is rewritten as:
Figure BDA0003053104370000042
let y (omega)i,ts) The solution representing the original system with parameter value ω, also belongs to manifold G (M, N)s) By thin singular value decomposition, y (ω)i,ts) Expressed as:
y(ω,ts)＝UΣVT (9)
according to LagrangianInterpolation method, truncation error epsilonN(ω,ts) Is defined as:
Figure BDA0003053104370000043
where the superscript (N +1) denotes the (N +1) th derivative of the parameter ω, expressed as
Figure BDA0003053104370000044
When N → ∞ is reached, the truncation error εN(ω,ts)→0；y(ωi,ts) The rewrite is:
Figure BDA0003053104370000045
in conjunction with equation (9), equation (11) is expressed as:
Figure BDA0003053104370000046
diagonal matrix sigma-diag [ sigma ]1,σ2,...,σM]Is y (ω)i,ts) Singular values of where σ1≥σ2≥,...,≥σM≧ 0, the singular value energy ratio ε is expressed as:
Figure BDA0003053104370000047
UΣVTis divided into main parts
Figure BDA0003053104370000048
And a truncated error part epsilonm(ω,ts) Substituting the main part and the error part into the formula (12) to obtain:
Figure BDA0003053104370000051
wherein epsilonm(ω,ts) And εN(ω,ts) Are all minor parts omitted, and therefore equation (14) is rewritten as:
Figure BDA0003053104370000052
autocorrelation matrix ScIs defined as:
Figure BDA0003053104370000053
solving a dimensionless control equation of the original system by a fourth-order Runge-Kutta (RK4) method;
inputting different modeling parameters to calculate displacement response, and storing response signals in a snapshot matrix
Figure BDA0003053104370000054
The autocorrelation matrix S is calculated according to equation (16)cComprises the following steps:
Figure BDA0003053104370000055
thirdly, by an autocorrelation matrix ScObtaining reduced order POM under new parameters
The reduced order POM under the parameter omega is calculated,
Figure BDA0003053104370000056
expressed as:
Figure BDA0003053104370000057
equation (17) is rewritten as:
Figure BDA0003053104370000058
by calculating a correlation matrix ScObtaining reduced-order POMs from the feature vectors;
four, order-reduced POM for reducing order of kinetic differential equation
The reduced POM is defined by the eigenvector of the autocorrelation matrix, and is sampled at a rotation speed ω' with a sampling length tsLet us order
Figure BDA0003053104370000059
The eigenvectors representing the autocorrelation matrix are arranged in descending order of the corresponding eigenvalues
Figure BDA00030531043700000510
Representing a linear subspace having n reduced-order POMs;
projecting the transient response of the complex high-dimensional space onto a low-dimensional subspace through coordinate transformation:
Figure BDA0003053104370000061
wherein the content of the first and second substances,
Figure BDA0003053104370000062
Figure BDA0003053104370000063
dimensionless displacement response, coordinate transformation matrix and POD modal coordinate of the original system are respectively, and the ordinary differential equation (10) is simplified as follows:
Figure BDA0003053104370000064
the IGM method with k points for the same ROM degrees of freedom only needs to perform singular value decomposition once, so that the computation efficiency of the IGM method is much higher than that of the ITGM method.
The singular value energy ratio epsilon is more than or equal to 95 percent.
When the dimensionless control equation is solved by a fourth-order Runge-Kutta (RK4) method, initial displacement and speedAre each Xi＝Yi＝Y17＝0.1，
Figure BDA0003053104370000065
The time step is 0.001 pi.
The invention has the beneficial effects that:
(1) the invention can not only obtain the optimal reduced order model of the complex rotor system, but also keep the orthogonal normalization characteristic of the reduced order POM and eliminate the local characteristic of the ITGM method.
(2) The method can be used in a wide frequency parameter range, and has the advantages of accurate calculation, high calculation efficiency and strong robustness. The method is not only suitable for the rotor bearing model, but also suitable for other dynamic models.
Drawings
FIG. 1 is a rotor bearing system for a 600MW turbine engine.
Figure 2 is a schematic view of a multi-disc rotor bearing system,
FIG. 3 is a graph of a Grassmann manifold based adaptive interpolation method, FIG. 3(a) is a four-point interpolation of the ITGM method; fig. 3(b) shows four-point interpolation by the IGM method.
FIG. 4 is a comparison of the two-point interpolation IGM and ITGM methods at a rotation speed ω 550 rad/s: FIG. 4(a) is a displacement response, with the horizontal axis representing the number of rotations and the vertical axis representing the longitudinal displacement; fig. 4(b) is an axis trace diagram, and the horizontal and vertical coordinates are the horizontal and vertical displacements of the rotating shaft.
Fig. 5 is a comparison between the IGM and the ITGM method for three-point interpolation when the rotation speed ω is 550rad/s, and fig. 5(a) shows the displacement response, with the horizontal axis representing the number of rotations and the vertical axis representing the longitudinal displacement; FIG. 5(b) is an axial trace diagram, in which the horizontal and vertical coordinates are the horizontal and vertical displacements of the rotating shaft.
FIG. 6 is a comparison of the five-point interpolated IGM and ITGM methods at a speed ω 550 rad/s: FIG. 6(a) is a displacement response, with the horizontal axis representing the number of rotations and the vertical axis representing the longitudinal displacement; fig. 6(b) is an axis trace diagram, and the horizontal and vertical coordinates are the horizontal and vertical displacements of the rotating shaft.
Fig. 7 is a convergence curve of the ROM response of the IGM method when the rotation speed ω is 550rad/s, fig. 7(a) is a convergence curve of the 2-point ROM response, fig. 7(b) is a convergence graph of the 3-point ROM response, fig. 7(c) is a convergence curve of the 5-point ROM response, the horizontal axis is the number of rotations, and the vertical axis is the longitudinal displacement.
Detailed Description
The invention is further illustrated with reference to the following figures and examples.
The method comprises the following specific implementation steps:
firstly, establishing a high-dimensional nonlinear rotor system model, calculating the transient response of a rotor by a four-order Runge Kutta method under different modeling parameter values, and inputting a response signal into a snapshot matrix
Figure BDA0003053104370000071
In (1).
Subsequently, an autocorrelation matrix S is calculatedcAnd by means of an autocorrelation matrix ScTo obtain the reduced-order POM after the model parameters are changed.
And finally, deriving a coordinate transformation matrix by the reduced-order POM, and simplifying a complex high-dimensional nonlinear rotor system differential equation through coordinate conversion to achieve the purpose of reducing the order.
Step one, obtaining a differential equation of motion of a high-dimensional nonlinear rotor system
Fig. 1 is a rotor bearing system for a 600MW turbine engine, which is a multi-disc rotor bearing system as shown in fig. 2, with a base loosening failure at the left support of fig. 2, with 14 discs on the shaft supported by two oil lubricated short journal bearings at both ends. Two oil lubrication short shaft journal bearings at two ends are used for supporting, a base loosening fault occurs at the left supporting position, and 14 discs are arranged on a shaft.
Assuming that the axes are elastic and the mass of each segment axis is concentrated at the centroid of the corresponding disk, the gyro moment has less influence on the system response because the distribution of the disks is symmetrical about the midpoint of the axes, and thus the gyro moment of the disks is not considered.
O1,O16Representing the geometric centers of the left and right bearings, respectively. The geometric center of each disk is Oi(i 2.. 15), the centroid is O (i 2.. 15), and the eccentricity is ei(i ═ 2.. 15). Equivalent mass of each disc and corresponding axial mass in mi(i 1.. 16) representsThe equivalent damping of each disc is ci(i 1.. 16). Further, the equivalent stiffness of the shaft between the individual turntables is ki(i 1.. 15), the horizontal and vertical displacements of the geometric center of the disk and the center of mass of the support end are denoted Xi,Yi(i 1.. 16), the nonlinear oil film forces in the horizontal and vertical directions are denoted as Fx,Fy. Assuming that the vibration of the loosening fault occurs only in the vertical direction, the corresponding displacement is Y17. Base klLoose stiffness and base clThe looseness damping of (a) is generally a piecewise linear function expressed as:
Figure BDA0003053104370000081
Figure BDA0003053104370000082
considering the nonlinear fluid film force, the unbalance force and the gravity of the whole system, the motion differential equation of the rotor bearing system is expressed as:
Figure BDA0003053104370000083
wherein X ═ X1,...X16,Y1...Y16,Y17]M, C, K are the mass, damping and stiffness matrices, respectively, and F is the resultant of the nonlinear fluid film force, the imbalance force and the gravity of each disk, expressed as:
Figure BDA0003053104370000084
M＝diag(m1,…,m16,m1,…m16,m17)
C＝diag(c1,…,c16,c1,…c16,c17)
Figure BDA0003053104370000085
in order to solve the transient response of the multi-degree-of-freedom system conveniently, the dynamic control equation of the system is usually subjected to dimensionless transformation, and the formula (3) is dimensionless by the following formula:
Figure BDA0003053104370000091
wherein f isxAnd fyIs a dimensionless non-linear oil film force,
Figure BDA0003053104370000092
is the Sommerfeld number, μ is the lube viscosity, c is the oil film thickness, L is the bearing length, R is the bearing radius, ω is the angular frequency, τ is the dimensionless time, and W is half the rotor system weight.
Obtaining the model of the step one for establishing the high-dimensional nonlinear rotor system
Secondly, inputting the response signal into the snapshot matrix and calculating the autocorrelation matrix
Based on the kinetic system theory, the entire solution set of the second order ordinary differential equation is a microfluidics, and the order of the microfluidics is greater than 2. Thus, the snapshot matrix of responses under different parameters can be considered as one point on the grassmann manifold. Manifold G (M, N)s) Is provided with a plurality of sampling time sequences NsThe M-dimensional column vector coverage of the corresponding solution.
N precomputed snapshot signal sets of dimension M system are set as
Figure BDA0003053104370000093
ωi∈[ωα,ωβ]Wherein each parameter ωiSample length t ofsThe same is true. Then, the snapshot set
Figure BDA0003053104370000094
Can be viewed as manifold G (M, N)s) N points of
Figure BDA0003053104370000095
Selecting a point S of the manifold0As a reference point, a snapshot matrix y with a corresponding parameter omega on a manifold can be obtained by a Lagrange interpolation method*(ω,ts) As shown in fig. 3. Expressed as:
Figure BDA0003053104370000096
according to Grassmann manifold correlation theory in differential geometry and singular value decomposition basic theory, manifold G (M, N)s) The upper interpolation point is represented as:
y(ωi,ts)＝UiΣiVi t,i＝0,...,N-1 (28)
wherein, sigmai∈RM×MIs a matrix of singular values and is,
Figure BDA0003053104370000101
Vi∈RM×Mis an orthogonal matrix. Then equation (7) can be rewritten as:
Figure BDA0003053104370000102
let y (omega)i,ts) Representing the solution of the original system with parameter values omega, which also belongs to manifold G (M, N)s). Using a thin singular value decomposition method, y (ω)i,ts) Expressed as:
y(ω,ts)＝UΣVT (30)
according to the Lagrange interpolation method, the truncation error epsilonN(ω,ts) Is defined as:
Figure BDA0003053104370000103
in which the superscript (N +1) denotes the parameter ωDerivative of the (N +1) th order, expressed as
Figure BDA0003053104370000104
When N → ∞ is reached, the truncation error εN(ω,ts) → 0, indicates that the more interpolation points, the closer to the solution of the original system at parameter ω. y (omega)i,ts) Is rewritten as
Figure BDA0003053104370000105
In conjunction with equation (9), equation (11) is expressed as:
Figure BDA0003053104370000106
diagonal matrix sigma-diag [ sigma ]1,σ2,...,σM]Is y (ω)i,ts) Singular values of where σ1≥σ2≥,...,≥σMIs more than or equal to 0. The singular value energy ratio epsilon is expressed as:
Figure BDA0003053104370000107
the energy contribution criterion is typically ε ≧ 95%, or even ε ≧ 99%, and the value of m is much smaller than the dimensionality of the original system (this condition is typically satisfied). Therefore, U Σ VTIs divided into main parts
Figure BDA0003053104370000111
And a truncated error part epsilonm(ω,ts) Substituting the main part and the error part into the formula (12) to obtain:
Figure BDA0003053104370000112
wherein epsilonm(ω,ts) And εN(ω,ts) All of the minor portions may be omitted, and therefore, formula (14) is rewrittenComprises the following steps:
Figure BDA0003053104370000113
Figure BDA0003053104370000114
the dimensionality is high, resulting in a snapshot matrix y*(ω,ts) The singular value decomposition of (a) can take a lot of time and memory. Autocorrelation matrix ScIs defined as:
Figure BDA0003053104370000115
the dimensionless control equation of the original system is solved by a fourth-order Runge-Kutta (RK4) method, which is a direct numerical integration method, and the initial displacement and the velocity are X respectivelyi＝Yi＝Y17＝0.1，
Figure BDA0003053104370000116
The time step is 0.001 pi. Inputting different modeling parameters to calculate displacement response, and storing response signals in a snapshot matrix
Figure BDA0003053104370000117
The snapshot matrix has NsThe snapshot matrix with rows and columns 33 sampled from the transient response may better reduce the model order. Because the transient response of the system contains free vibration, forced vibration and more dynamic information in the original system. Therefore, the present invention chooses to obtain a snapshot matrix with different frequency parameter values from the transient response. The autocorrelation matrix S can be calculated according to the formula (16)c：
Figure BDA0003053104370000118
Thirdly, by an autocorrelation matrix ScObtaining reduced order POM under new parameters
The reduced order POM under the parameter omega is calculated,
Figure BDA0003053104370000119
expressed as:
Figure BDA00030531043700001110
equation (17) is rewritten as:
Figure BDA00030531043700001111
as can be seen from the above, by calculating the correlation matrix ScThe feature vectors of (a) obtain reduced order POMs. Despite the snapshot matrix
Figure BDA0003053104370000121
It is not the solution of the original system under parameter ω, but it is located in the vicinity of the solution of the original system. In addition, due to the snapshot matrix
Figure BDA0003053104370000122
Multiple motion states (e.g., quasi-periodic motion, chaotic motion) are contained at corresponding interpolation points in the wide parameter region, and thus, the motion state is derived from the snapshot matrix
Figure BDA0003053104370000123
The obtained reduced-order POM has richer modality information than from the case of the ordinary snapshot matrix. Thus, the adaptive POM method not only can obtain the optimal POM of a complex physical system, but also can keep the orthogonal normalization of the reduced POM. Meanwhile, the method can be used in a wide parameter range because the limitation of the local property of the point tangent space on the manifold is eliminated.
Four, order-reduced POM for reducing order of kinetic differential equation
The reduced POM is defined by the eigenvectors of the autocorrelation matrix. Let us assume that the sampling is done at a rotational speed ω' with a sampling length ts. Order to
Figure BDA0003053104370000124
Feature vectors representing the autocorrelation matrix, the feature vectors being arranged in descending order of the corresponding feature values. Is provided with
Figure BDA0003053104370000125
Representing a linear subspace having n reduced-order POMs.
Through coordinate transformation, the transient response of a complex high-dimensional space can be projected onto a low-dimensional subspace:
Figure BDA0003053104370000126
wherein the content of the first and second substances,
Figure BDA0003053104370000127
Figure BDA0003053104370000128
the dimensionless displacement response, coordinate transformation matrix and POD modal coordinates of the original system are respectively. Ordinary differential equation (10) is simplified as:
Figure BDA0003053104370000129
fifthly, calculating results of IGM method
For simple motion cases, as shown in fig. 4-6, the displacement response and the axis trajectory of the ROM are obtained by interpolation of IGM and ITGM methods at 2, 3, and 5 points, and compared to the results for ω being 500 rad/s. In fig. 4, 17 reduced-order POMs are obtained by the IGM method during two-point interpolation, and then a response is calculated for the POMs, which is consistent with the response of the original system. In fig. 5, the response consistent with the original system can be calculated by using only 12 reduced POMs by using the IGM method during three-point interpolation. In fig. 6, only 4 reduced-order POMs are needed for 5-point interpolation to perfectly match the response of the original system.
The convergent response curves of the ROM with different degrees of freedom at a rotation speed ω 550rad/s are obtained by the IGM methods of 2 point, 3 point and 5 point, respectively, as shown in fig. 7. It can be concluded that the more degrees of freedom of the ROM, the more accurate the response obtained using the IGM method at 3 different interpolation points. In addition, as the number of interpolation points increases, the DOF number of the ROM decreases.
During complex movements, the ITGM and IGM methods can be found to be equivalent in the local parameter region. However, since the ITGM method cannot calculate dense interpolation points in a wider parameter area, it is not suitable for complex motion in the case of a wider parameter area.
For the same ROM degree of freedom, the ITGM method of the k point needs to execute singular value decomposition of k +1, the IGM method of the k point only needs to execute singular value decomposition once, and meanwhile, the calculation time of the IGM five-point method is more than 20 times faster than the numerical integration speed of an original system, so that the calculation efficiency of the IGM method is far higher than that of the ITGM method.

Claims (3)

1. A self-adaptive orthogonal decomposition method for reducing the order of a multi-disk rotor system model is characterized by comprising the following steps:
step one, obtaining a differential equation of motion of a high-dimensional nonlinear rotor system;
the rotor bearing system is a multi-disc rotor bearing system, is supported by two oil lubrication short shaft journal bearings at two ends, has a base loosening fault at a left support position, and is provided with 14 discs on a shaft;
the shafts are elastic, and the mass of each section of shaft is concentrated on the mass center of the corresponding disc, and the equivalent rigidity of the shaft between the turntables is k without considering the gyro moment of the disciI 1.. 15, the horizontal and vertical displacements of the geometric center of the disk and the center of mass of the support end are denoted Xi,YiI 1.. 16, the nonlinear oil film force in the horizontal and vertical directions is denoted as Fx,FyVibration with loosening failure only in the vertical direction, corresponding to a displacement of Y17Base klLoose stiffness and base clIs a piecewise linear function expressed as:
Figure FDA0003053104360000011
Figure FDA0003053104360000012
the differential equation of motion for a rotor bearing system is expressed as:
Figure FDA0003053104360000013
wherein X ═ X1,...X16,Y1...Y16,Y17]M, C, K are the mass, damping and stiffness matrices, respectively, and F is the resultant of the nonlinear fluid film force, the imbalance force and the gravity of each disk, expressed as:
Figure FDA0003053104360000014
M＝diag(m1,…,m16,m1,…m16,m17)
C＝diag(c1,…,c16,c1,…c16,c17)
Figure FDA0003053104360000021
carrying out dimensionless transformation on a dynamic control equation of the system, wherein the formula (3) is dimensionless by the following formula:
Figure FDA0003053104360000022
wherein f isxAnd fyIs dimensionlessThe non-linear oil film force is,
Figure FDA0003053104360000023
is the Sommerfeld number, μ is the lubricant viscosity, c is the oil film thickness, L is the bearing length, R is the bearing radius, ω is the angular frequency, τ is the dimensionless time, W is half the rotor system weight;
thereby establishing a high-dimensional nonlinear rotor system model;
inputting the response signal into a snapshot matrix and calculating an autocorrelation matrix;
manifold G (M, N)s) Is provided with a plurality of sampling time sequences NsThe M-dimensional column vector coverage of the corresponding solution;
n precomputed snapshot signal sets of dimension M system are set as
Figure FDA0003053104360000024
Wherein each parameter ωiSample length t ofsSame, snapshot set
Figure FDA0003053104360000025
Viewed as manifold G (M, N)s) N points of
Figure FDA0003053104360000026
Selecting a point S of the manifold0As a reference point, obtaining a snapshot matrix y with a corresponding parameter omega on a manifold by a Lagrange interpolation method*(ω,ts) Expressed as:
Figure FDA0003053104360000031
manifold G (M, N)s) The upper interpolation point is represented as:
y(ωi,ts)＝UiΣiVi t,i＝0,...,N-1 (7)
wherein, sigmai∈RM×MIs a matrix of singular values and is,
Figure FDA0003053104360000032
Vi∈RM×Mis an orthogonal matrix, then equation (7) is rewritten as:
Figure FDA0003053104360000033
let y (omega)i,ts) The solution representing the original system with parameter value ω, also belongs to manifold G (M, N)s) By thin singular value decomposition, y (ω)i,ts) Expressed as:
y(ω,ts)＝UΣVT (9)
according to the Lagrange interpolation method, the truncation error epsilonN(ω,ts) Is defined as:
Figure FDA0003053104360000034
where the superscript (N +1) denotes the (N +1) th derivative of the parameter ω, expressed as
Figure FDA0003053104360000035
When N → ∞ is reached, the truncation error εN(ω,ts)→0；y(ωi,ts) The rewrite is:
Figure FDA0003053104360000036
in conjunction with equation (9), equation (11) is expressed as:
Figure FDA0003053104360000037
diagonal matrix sigma-diag [ sigma ]1,σ2,...,σM]Is y (ω)i,ts) Singular values of where σ1≥σ2≥,...,≥σM≧ 0, the singular value energy ratio ε is expressed as:
Figure FDA0003053104360000041
UΣVTis divided into main parts
Figure FDA0003053104360000042
And a truncated error part epsilonm(ω,ts) Substituting the main part and the error part into the formula (12) to obtain:
Figure FDA0003053104360000043
wherein epsilonm(ω,ts) And εN(ω,ts) Are all minor parts omitted, and therefore equation (14) is rewritten as:
Figure FDA0003053104360000044
autocorrelation matrix ScIs defined as:
Figure FDA0003053104360000045
solving a dimensionless control equation of the original system by a fourth-order Runge-Kutta (RK4) method;
inputting different modeling parameters to calculate displacement response, and storing response signals in a snapshot matrix
Figure FDA0003053104360000046
The autocorrelation matrix S is calculated according to equation (16)cComprises the following steps:
Figure FDA0003053104360000047
thirdly, by an autocorrelation matrix ScObtaining reduced order POM under new parameters
The reduced order POM under the parameter omega is calculated,
Figure FDA0003053104360000048
expressed as:
Figure FDA0003053104360000049
equation (17) is rewritten as:
Figure FDA00030531043600000410
by calculating a correlation matrix ScObtaining reduced-order POMs from the feature vectors;
four, order-reduced POM for reducing order of kinetic differential equation
The reduced POM is defined by the eigenvector of the autocorrelation matrix, and is sampled at a rotation speed ω' with a sampling length tsLet us order
Figure FDA00030531043600000411
The eigenvectors representing the autocorrelation matrix are arranged in descending order of the corresponding eigenvalues
Figure FDA0003053104360000051
Representing a linear subspace having n reduced-order POMs;
projecting the transient response of the complex high-dimensional space onto a low-dimensional subspace through coordinate transformation:
Figure FDA0003053104360000052
wherein the content of the first and second substances,
Figure FDA0003053104360000053
Figure FDA0003053104360000054
dimensionless displacement response, coordinate transformation matrix and POD modal coordinate of the original system are respectively, and the ordinary differential equation (10) is simplified as follows:
Figure FDA0003053104360000055
the IGM method with k points for the same ROM degrees of freedom only needs to perform singular value decomposition once, so that the computation efficiency of the IGM method is much higher than that of the ITGM method.
2. The method of claim 1, wherein the step-down adaptive orthogonal decomposition comprises:
the singular value energy ratio epsilon is more than or equal to 95 percent.
3. The method of claim 1, wherein the step-down adaptive orthogonal decomposition comprises:
when the dimensionless control equation is solved by a fourth-order Runge-Kutta method, the initial displacement and the speed are respectively Xi＝Yi＝Y17＝0.1，
Figure FDA0003053104360000056
The time step is 0.001 pi.
CN202110492921.4A 2021-05-07 2021-05-07 Adaptive orthogonal decomposition method for model reduction of multi-disk rotor system Active CN113254863B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110492921.4A CN113254863B (en) 2021-05-07 2021-05-07 Adaptive orthogonal decomposition method for model reduction of multi-disk rotor system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110492921.4A CN113254863B (en) 2021-05-07 2021-05-07 Adaptive orthogonal decomposition method for model reduction of multi-disk rotor system

Publications (2)

Publication Number Publication Date
CN113254863A true CN113254863A (en) 2021-08-13
CN113254863B CN113254863B (en) 2023-11-07

Family

ID=77223829

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110492921.4A Active CN113254863B (en) 2021-05-07 2021-05-07 Adaptive orthogonal decomposition method for model reduction of multi-disk rotor system

Country Status (1)

Country Link
CN (1) CN113254863B (en)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE10137909A1 (en) * 2001-01-18 2002-08-08 Siemens Ag Process for the simulation of mechatronic systems
US20060149525A1 (en) * 2005-01-05 2006-07-06 Chang Gung University Multi-point model reductions of VLSI interconnects using the rational arnoldi method with adaptive orders
DE102013007621A1 (en) * 2013-05-02 2014-11-06 André Schleicher Bearing reluctance motor with identical slot pitch in stand and rotor
CN111339706A (en) * 2020-03-09 2020-06-26 西北工业大学 POD-based rotor-bearing system model secondary order reduction method
CN111368447A (en) * 2020-03-09 2020-07-03 西北工业大学 Improved transient time sequence-based nonlinear POD dimension reduction method
CN111563323A (en) * 2020-04-27 2020-08-21 西北工业大学 Transient POD method based on minimum error of bifurcation parameter

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE10137909A1 (en) * 2001-01-18 2002-08-08 Siemens Ag Process for the simulation of mechatronic systems
US20060149525A1 (en) * 2005-01-05 2006-07-06 Chang Gung University Multi-point model reductions of VLSI interconnects using the rational arnoldi method with adaptive orders
DE102013007621A1 (en) * 2013-05-02 2014-11-06 André Schleicher Bearing reluctance motor with identical slot pitch in stand and rotor
CN111339706A (en) * 2020-03-09 2020-06-26 西北工业大学 POD-based rotor-bearing system model secondary order reduction method
CN111368447A (en) * 2020-03-09 2020-07-03 西北工业大学 Improved transient time sequence-based nonlinear POD dimension reduction method
CN111563323A (en) * 2020-04-27 2020-08-21 西北工业大学 Transient POD method based on minimum error of bifurcation parameter

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
姚伟刚;徐敏;叶茂;: "基于特征正交分解的非定常气动力建模技术", 力学学报, no. 04 *
徐敏;安效民;曾宪昂;陈士橹;: "基于BPOD的气动弹性降阶及其在主动控制中的应用", 中国科技论文在线, no. 10 *

Also Published As

Publication number Publication date
CN113254863B (en) 2023-11-07

Similar Documents

Publication Publication Date Title
Jin et al. Nonlinear dynamic analysis of a complex dual rotor-bearing system based on a novel model reduction method
Harvey III et al. Steady and unsteady computation of impeller‐stirred reactors
Gupta Advanced dynamics of rolling elements
Hashish et al. Finite element and modal analyses of rotor-bearing systems under stochastic loading conditions
Erkaya Prediction of vibration characteristics of a planar mechanism having imperfect joints using neural network
CN110532625A (en) Aero-engine turbine disk-twin the modeling method of rotor-bearing system number
Zhang et al. Elastic ring deformation and pedestal contact status analysis of elastic ring squeeze film damper
Xu et al. Vibration characteristics of bearing-rotor systems with inner ring dynamic misalignment
Rashidi et al. Bifurcation and nonlinear dynamic analysis of a rigid rotor supported by two-lobe noncircular gas-lubricated journal bearing system
Rahnejat Computational modelling of problems in contact dynamics
CN111339706A (en) POD-based rotor-bearing system model secondary order reduction method
Meybodi et al. Numerical analysis of a rigid rotor supported by aerodynamic four-lobe journal bearing system with mass unbalance
Ding et al. Dynamic modeling of ultra-precision fly cutting machine tool and the effect of ambient vibration on its tool tip response
Pham et al. Efficient techniques for the computation of the nonlinear dynamics of a foil-air bearing rotor system
CN113434983B (en) Rapid calculation method for nonlinear dynamic characteristics of sliding bearing rotor system
Chen et al. Investigations on the dynamic characteristics of a planar slider-crank mechanism for a high-speed press system that considers joint clearance
Knoll et al. Full dynamic analysis of crankshaft and engine block with special respect to elastohydrodynamic bearing coupling
CN113254863A (en) Self-adaptive orthogonal decomposition method for order reduction of multi-disk rotor system model
CN111368447A (en) Improved transient time sequence-based nonlinear POD dimension reduction method
HASHIMOTO et al. Performance characteristics of elliptical journal bearings in turbulent flow regime
Chen et al. Design and dynamics modeling of a novel 2R1T 3-DOF parallel motion simulator
Taylor et al. Closed-form, steady-state solution for the unbalance response of a rigid rotor in squeeze film damper
Zhang et al. Investigating the dynamic characteristics of the hydrostatic bearing spindle system with active piezoelectric restrictors
Pourashraf et al. A new Galerkin Reduction approach for the analysis of a fully coupled foil air bearing rotor system with bilinear foil model
CN115935687A (en) Method for calculating coupling lubrication and dynamic characteristic parameters of flanged bearing

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