CN112632849A - Mechanical seal state analysis method based on numerical model - Google Patents
Mechanical seal state analysis method based on numerical model Download PDFInfo
- Publication number
- CN112632849A CN112632849A CN202011453641.4A CN202011453641A CN112632849A CN 112632849 A CN112632849 A CN 112632849A CN 202011453641 A CN202011453641 A CN 202011453641A CN 112632849 A CN112632849 A CN 112632849A
- Authority
- CN
- China
- Prior art keywords
- parameter
- model
- mechanical seal
- state
- optimization
- 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.)
- Pending
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 23
- 238000005457 optimization Methods 0.000 claims abstract description 74
- 238000000034 method Methods 0.000 claims abstract description 43
- 238000004088 simulation Methods 0.000 claims abstract description 8
- 238000004422 calculation algorithm Methods 0.000 claims description 25
- 230000003068 static effect Effects 0.000 claims description 16
- 239000012530 fluid Substances 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 10
- 238000002922 simulated annealing Methods 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000011144 upstream manufacturing Methods 0.000 claims description 6
- 230000002068 genetic effect Effects 0.000 claims description 3
- 239000002245 particle Substances 0.000 claims description 3
- 238000012706 support-vector machine Methods 0.000 claims description 3
- 238000012423 maintenance Methods 0.000 abstract description 9
- 238000007789 sealing Methods 0.000 description 11
- 238000009826 distribution Methods 0.000 description 10
- 239000011159 matrix material Substances 0.000 description 8
- 238000013461 design Methods 0.000 description 5
- 230000009191 jumping Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000013016 damping Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000006073 displacement reaction Methods 0.000 description 3
- 238000005461 lubrication Methods 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- 230000007704 transition Effects 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000007667 floating Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 238000013519 translation Methods 0.000 description 2
- 230000014616 translation Effects 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 238000000137 annealing Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000005314 correlation function Methods 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002706 hydrostatic effect Effects 0.000 description 1
- 238000009776 industrial production Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000012933 kinetic analysis Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 230000001953 sensory effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003746 surface roughness Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 238000011179 visual inspection Methods 0.000 description 1
Images
Classifications
-
- 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/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- 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
Abstract
The invention relates to a mechanical seal state analysis method based on a numerical model, which comprises the following steps: s1, confirming at least one undetermined state parameter to be analyzed in the mechanical seal, and giving a value range of the undetermined state parameter; s2 obtaining at least one online observation parameter in the mechanical seal; s3, establishing a numerical model for simulating the working state of the mechanical seal, wherein the numerical model takes the parameter to be determined as an input parameter and outputs a simulation observation parameter; s4, adjusting the value of the parameter to be determined within the value range by using an optimization method until the difference between the simulated observation parameter and the online observation parameter is within the error allowable range, and recording the current parameter to be determined as the target state parameter. The mechanical seal state analysis method based on the numerical model can conjecture the state parameters which cannot be directly monitored on line, further know the real-time state of the friction pair, and is convenient for implementing accurate and efficient maintenance.
Description
Technical Field
The invention relates to the technical field of sealing state simulation, in particular to a mechanical sealing state analysis method based on a numerical model.
Background
In some instances, it is desirable to prevent or reduce the migration of species between two adjacent spaces, such migration being referred to as "leakage", and "sealing" refers to an effective means of preventing or reducing leakage. Unlike the common simple static seal in life, the situation that sealing needs to be implemented between parts moving relatively is encountered in industrial production, and the mechanical seal has good performance under high parameters, wherein the dry gas seal is particularly suitable for the working condition with high parameters. A mechanical seal is a face-to-face dynamic seal that requires reduced or eliminated frictional wear of the friction pair (formed by the two faces in relative motion and the fluid medium) to extend life while maintaining low or no leakage. Due to the fact that the principle is complex, the structure is compact, the working state of the mechanical seal in the operation process is difficult to know, a user is forced to adopt maintenance strategies with high cost such as regular maintenance, and accurate and efficient maintenance cannot be implemented based on the real-time state of the friction pair.
Disclosure of Invention
Therefore, it is necessary to provide a mechanical seal state analysis method based on a numerical model for simulating a real-time friction state of a friction pair, aiming at the problem that accurate and efficient maintenance cannot be implemented based on the real-time state of the friction pair at present.
A mechanical seal state analysis method based on a numerical model comprises the following steps:
s1, confirming at least one undetermined state parameter to be analyzed in the mechanical seal, and giving a value range of the undetermined state parameter;
s2 obtaining at least one online observation parameter in the mechanical seal;
s3, establishing a numerical model for simulating the working state of the mechanical seal, wherein the numerical model takes the parameter to be determined as an input parameter and outputs a simulation observation parameter;
s4, adjusting the value of the parameter to be determined within the value range by using an optimization method until the difference between the simulated observation parameter and the online observation parameter is within the error allowable range, and recording the current parameter to be determined as the target state parameter.
In one embodiment, the mechanical seal to be analyzed comprises a static ring and a dynamic ring, wherein the static ring is a compensation ring, and the dynamic ring is a non-compensation ring; the undetermined state parameters comprise dynamic ring end jump and static ring bearing abnormity, and the static ring bearing abnormity is expressed by two components of axial force and deflection moment.
In one embodiment, the online observation parameter comprises an average fluid leakage rate, and the online observation parameter further comprises an end face contact state between the static ring and the dynamic ring.
In one embodiment, a rotating speed parameter of the rotating shaft and a pressure parameter on and downstream of the mechanical seal end surface are measured, and the rotating speed parameter and the pressure parameter are used as determination parameters and input into the numerical model.
In one embodiment, a rotating speed sensor is used for measuring the rotating speed of the rotating shaft, a pressure sensor is used for measuring pressure parameters on the upstream and downstream of the end face of the mechanical seal, a flow sensor is used for measuring the flow of fluid and calculating the average leakage rate of the fluid, and an acoustic emission sensor is used for detecting the end face contact state between the movable ring and the static ring.
In one embodiment, a global optimization method is adopted to adjust the value of the parameter to be determined until the difference value between the simulated observation parameter and the online observation parameter is within an error allowable range; the global optimization method adopts a simulated annealing algorithm, a genetic algorithm and/or a particle swarm algorithm and the like.
In one embodiment, the global optimization method comprises a first precision optimization model and a second precision optimization model, wherein the second precision optimization model has higher precision than the first precision optimization model; firstly, a first precision optimization model is adopted for global optimization to obtain a preliminary optimization result, then the preliminary optimization result is used as an initial value, a second precision optimization model is adopted for local optimization to obtain a final optimization result, and a target state parameter is determined.
In one embodiment, the numerical model for simulating the working state of the mechanical seal comprises a physical model and a proxy model, wherein the physical model and the proxy model respectively take the parameter to be determined as an input parameter, and the physical model and the proxy model respectively output a simulation observation parameter; the physical model is directly established based on the working state of the mechanical seal, the proxy model is established based on the calculation result set of the physical model, and the calculation speed of the proxy model is higher than that of the physical model.
In one embodiment, the proxy model is built based on the physical model by a polynomial fitting, support vector machine regression, and/or gaussian random process.
In one embodiment, the target state parameter is found by a Bayesian optimization method.
The mechanical seal state analysis method based on the numerical model can conjecture the state parameters which cannot be directly monitored on line, further know the real-time state of the friction pair, and is convenient for implementing accurate and efficient maintenance. The mechanical seal state analysis method based on the numerical model can be used for estimating future performance change and sudden failure risk and providing basis for making a maintenance plan and automatically responding to an emergency. The mechanical sealing state analysis method based on the numerical model is based on the numerical model based on the sealing working principle, and the obtained result has good interpretability.
Drawings
Fig. 1 is a flowchart illustrating a method for analyzing a mechanical seal state based on a numerical model according to an embodiment of the present invention;
FIG. 2 is a schematic structural diagram of a mechanical seal device according to an embodiment of the present invention;
FIG. 3 is a schematic structural diagram of a moving ring and a stationary ring according to an embodiment of the present invention;
FIG. 4 is a flowchart of a simulated annealing algorithm provided in accordance with an embodiment of the present invention;
fig. 5 is a flowchart illustrating a method for hierarchical optimization according to an embodiment of the present invention to estimate a parameter of an undetermined state;
fig. 6 is a flowchart of estimating a parameter of an undetermined state based on a proxy model and bayesian optimization according to an embodiment of the present invention;
fig. 7 is a diagram of an idea for building a physical model according to an embodiment of the present invention.
Wherein: 100. a moving ring; 200. a stationary ring; 400. an acoustic emission sensor.
Detailed Description
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in detail below. In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present invention. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein.
In the description of the present invention, it is to be understood that the terms "central," "longitudinal," "lateral," "length," "width," "thickness," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," "outer," "clockwise," "counterclockwise," "axial," "radial," "circumferential," and the like are used in the orientations and positional relationships indicated in the drawings for convenience in describing the invention and to simplify the description, and are not intended to indicate or imply that the referenced device or element must have a particular orientation, be constructed and operated in a particular orientation, and are not to be considered limiting of the invention.
Furthermore, the terms "first", "second" and "first" are used for descriptive purposes only and are not to be construed as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one such feature. In the description of the present invention, "a plurality" means at least two, e.g., two, three, etc., unless specifically limited otherwise.
In the present invention, unless otherwise expressly stated or limited, the terms "mounted," "connected," "secured," and the like are to be construed broadly and can, for example, be fixedly connected, detachably connected, or integrally formed; can be mechanically or electrically connected; they may be directly connected or indirectly connected through intervening media, or they may be connected internally or in any other suitable relationship, unless expressly stated otherwise. The specific meanings of the above terms in the present invention can be understood by those skilled in the art according to specific situations.
In the present invention, unless otherwise expressly stated or limited, the first feature "on" or "under" the second feature may be directly contacting the first and second features or indirectly contacting the first and second features through an intermediate. Also, a first feature "on," "over," and "above" a second feature may be directly or diagonally above the second feature, or may simply indicate that the first feature is at a higher level than the second feature. A first feature being "under," "below," and "beneath" a second feature may be directly under or obliquely under the first feature, or may simply mean that the first feature is at a lesser elevation than the second feature.
It will be understood that when an element is referred to as being "secured to" or "disposed on" another element, it can be directly on the other element or intervening elements may also be present. When an element is referred to as being "connected" to another element, it can be directly connected to the other element or intervening elements may also be present. The terms "vertical," "horizontal," "upper," "lower," "left," "right," and the like as used herein are for illustrative purposes only and do not denote a unique embodiment.
The invention provides an analysis method for predicting a real-time friction state of a friction pair in a mechanical sealing device. Mechanical seals include contact mechanical seals, hydrodynamic/hydrostatic mechanical seals, and dry gas seals. The following embodiments are described with respect to dry gas seals, it being understood that the principles and steps may be applied to other types of mechanical seals with reasonable modifications.
In an embodiment of the present invention, a method for analyzing a mechanical seal state based on a numerical model includes the following steps: s1, confirming at least one undetermined state parameter to be analyzed in the mechanical seal, and giving a value range of the undetermined state parameter; s2 obtaining at least one online observation parameter in the mechanical seal; s3, establishing a numerical model for simulating the working state of the mechanical seal, wherein the numerical model takes the parameter to be determined as an input parameter and outputs a simulation observation parameter; s4, adjusting the value of the parameter to be determined within the value range by using an optimization method until the difference between the simulated observation parameter and the online observation parameter is within the error allowable range, and recording the current parameter to be determined as the target state parameter. The general analysis idea of the method for analyzing a mechanical seal state based on a numerical model provided in this embodiment is shown in fig. 1.
It should be noted that the undetermined state parameter refers to a plurality of state parameters which cannot be observed online, may be inconsistent with an offline measurement state, and may significantly affect the operation of the mechanical seal, and an undetermined state parameter list and possible value ranges thereof need to be given according to experience; the on-line observation parameters may be obtained by a sensory inspection (such as visual inspection) without using a tool, or by a sensor mounted on the mechanical seal device, and only the on-line observation parameters that may be changed due to a change in the parameter to be determined are considered in this embodiment.
The mechanical seal state analysis method based on the numerical model can conjecture the state parameters which cannot be directly monitored on line, further know the real-time state of the friction pair, and is convenient for implementing accurate and efficient maintenance. The mechanical seal state analysis method based on the numerical model can be used for estimating future performance change and sudden failure risk and providing basis for making a maintenance plan and automatically responding to an emergency. The mechanical sealing state analysis method based on the numerical model is based on the numerical model based on the sealing working principle, and the obtained result has good interpretability.
In dry gas seals (and more generally mechanical seals), the rotating one of the two seal rings is referred to as the rotating ring 100 and the non-rotating one is referred to as the stationary ring 200. One of the two seals is mounted floating, called the compensator ring, and the other is mounted fixed, called the non-compensator ring (in special cases it may happen that both seal rings are compensator rings, this embodiment is not discussed).
In an embodiment of the present invention, as shown in fig. 2-3, when the stationary ring 200 is a compensation ring and the moving ring 100 is an uncompensated ring, the stationary ring 200 is constrained not to rotate but can still move along another 5 degrees of freedom (in general, the mechanical field considers that a rigid body has 6 degrees of freedom, including 3 translations and 3 rotations), but in these 5 degrees of freedom, only 3 are considered in the classical mechanical sealing dynamics: axial (Z-axis) translation, rotation about the X-axis, and rotation about the Y-axis, as is the case in this embodiment. In addition, the stationary ring 200 may be a non-compensation ring, and the moving ring 100 may be a compensation ring, and the stationary ring 200 may be completely fixed. The special design of the dry gas seal allows the compensating ring to follow the movement of the non-compensating ring in an attempt to keep the shape of the gas film stable (too thick a film would result in excessive leakage, and a film less than a certain value would mean solid contact).
As an implementation manner, on a set of spiral groove dry gas seal (a specific form of mechanical seal) equipment, a rotating speed sensor is used for measuring the rotating speed of a rotating shaft, a pressure sensor (one of an upstream and a downstream) is used for measuring pressure parameters of the upstream and the downstream of the end surface of the mechanical seal, a flow sensor is used for measuring the flow of fluid and calculating the average leakage rate of the fluid, and an acoustic emission sensor 400 is used for detecting an end surface contact dynamic signal between the moving ring 100 and the static ring 200 by means of filtering.
Firstly, confirming undetermined state parameters and an online observation mode of a mechanical seal to be analyzed. In this example, the rotation speed and the upstream and downstream pressures are considered as the determined parameters in the model. The measurements of the flow sensor and acoustic emission sensor 400 are considered to be on-line observed parameters. The undetermined state parameters comprise the end jump of the dynamic ring 100 and the bearing abnormity of the static ring 200, and the support abnormity of the static ring 200 is expressed by two components of axial force and yawing moment (without loss of generality, and the influence of the yawing moment phase is not considered).
In a first specific embodiment of the invention, the numerical model is, for example, y (y) and (x), where x [ γ, F, M ═ y (x) ]]To be determined, y (x) ═ q (t; x), σ (t; x)]Are observations predicted by the model, where q (t), σ (t) are the leakage rate and the root mean square of the short-term acoustic emission after band-pass filtering (170kHz + -40 kHz), respectively. Based on the dry gas seal end face lubrication mechanism and the kinetic analysis of the compensating ring, a physical equation (given below) can be established. These equations can be solved numerically either by writing a program or using commercially available finite element software. Monitoring the leak rate does not yield its immediate value, but rather an average. The online observation y' includes
And σ' (t).
And establishing a difference measurement rule for the monitoring result and the simulation result. The difference is defined as
KqAnd KσAre weight coefficients. In an algorithmic implementation, the integrals in the above equations are replaced with numerical integrals.
By adopting a global optimization method suitable for multi-modal and derivative-free problems, the embodiment takes a simulated annealing algorithm as an example (a genetic algorithm, a particle swarm algorithm, etc. can also be used), as shown in fig. 4:
(1) initialization
Let the iteration number i equal to 0. Determining an initial set of x, denoted x0. Calculating an objective function L at an initial value0＝L(x0(ii) a y'). Determination of the initial temperature T0(the "temperature" in the simulated annealing algorithm does not represent a physical temperature, and accordingly, "simulated annealing" does not represent a physical annealing process, and only means that the algorithm was designed based on the physical phenomenon);
(2) generating a new state
According to the state transition rule, from xiRandomly generating xi+1. The present example employs the following state transition rules:
first of all, calculate
Wherein Δγ，ΔF，ΔMIs in the range of [ -0.5, 0.5 [)]Upper uniformly distributed independent random variable, constant Rγ，RF，RMRepresenting the transfer range scale. For gamma*，F*，M*Until it enters the feasible domain, thereby obtaining x*＝[γ*，F*，M*]For determining whether to accept the new state.
Easy to find, TiThe smaller the range of state transitions.
(3) Determining whether to accept the new state
Calculating L*＝L(x*(ii) a y'), calculating Δ L ═ L*-LiAnd determining according to the value of the delta L:
if Δ L is less than or equal to 0, accepting a new solution;
if Δ L > 0, a new solution is accepted with probability p ═ exp (- Δ L/T).
Ordering x upon acceptance of new solutioni+1＝x*，Li+1＝L*Otherwise, let xi+1＝xi，Li+1＝Li. In the algorithm implementation, x, which historically achieved the smallest L, is always retained.
(4) Let Ti+1＝βTi. Where β is a constant, this example may take 0.995.
(5) Let i equal i + 1. And (4) judging whether the termination condition is met, and if not, jumping to (2).
Solving for x, which minimizes L (x; y') by the above method, is recorded as
This is the estimation given by the method to the state parameter to be determined.
In a second specific embodiment of the present invention, a hierarchical optimization strategy is employed based on the first embodiment. Two optimization models are respectively established, as shown in fig. 5, including a low-precision optimization model and a high-precision optimization model. In other embodiments described in this embodiment, if not specifically stated, a high-precision optimization model is used.
When the optimization is implemented, the result obtained by the low-precision optimization model is optimized by using a global optimization method, and then the optimization result is used as an initial value, and the high-precision optimization model is used for implementing the optimization by using a local optimization method.
It should be noted that, a local optimization algorithm usually searches based on a local gradient, and can converge with a small number of iteration steps, but the converged value is usually only guaranteed to be a local extremum; the global optimization algorithm requires more iteration steps, but it has the ability to search for the minimum (or maximum) value over the entire domain of definition.
Here, a low-precision optimization model with higher efficiency is used in a global optimization algorithm that requires a larger number of iteration steps; finally, to obtain an accurate result, a high-precision optimization model is needed, which is inefficient, but since the previous optimization step finds several possible approximate positions of the global optimum, local optimization is only needed, and thus a local optimization algorithm with a small number of iterations can be used.
Specifically, the method comprises the following steps:
a calculation program with normal algorithm parameters (parameters for controlling the accuracy of a numerical algorithm under a certain physical model, including grid fineness, time step length and tolerance) is used as a high-accuracy optimization model, and a calculation program with more loose algorithm parameters (namely, thicker grid, larger time step length and larger tolerance) is used as a low-accuracy optimization model.
The output of the low-precision optimization model is first optimized using a simulated annealing algorithm (as in the first embodiment) (i.e., L (x;y') is adopted, and in the process, not only x corresponding to the history minimum value of L is returned, but also iteration history is returned to be a set H { (x)i，Li) I ═ 0, 1, 2,. and n (where n is the number of iterations) and the optimization continues with a high accuracy optimization model as follows:
(1) initializing selected low-precision optimization model iteration history set table
Local optimization result set of high-precision optimization model
(2) (x, L) with the smallest L is selected from H and removed from H;
(3) and judging whether the (x, L) selected in the last step is close to at least one element in the C (wherein the 'close' is defined as being smaller than a preset threshold value in at least one dimension of gamma, F and M). If not, adding the (x, L) into C;
(4) using the (x, L) as an initial value, and adopting a local optimization algorithm to find the (x) which obtains the minimum value under a high-precision optimization model*，L*) And is recorded in R;
(5) and judging whether the termination condition is met. And if not, jumping to (2).
As an implementation of the local optimization algorithm in the step (4), the newton method based on finite difference gradient estimation includes the following steps:
(1) giving an initial value (x, L);
(2) performing finite difference estimation on the Hessian matrix H and the gradient g of the function l (x) around x ═ γ, F, M ];
(3) x-H-1g is taken as the new value of x;
(4) and judging whether the termination condition is met. And if not, jumping to (2).
Finally finding the x which obtains the minimum L under the high-precision optimization model as the x estimated by the method
In a third specific embodiment of the present invention, based on the first embodiment, as shown in fig. 6, the number of times of calculation of the physical model is reduced as much as possible by using a proxy model technique and a bayesian optimization technique. Since the physical model is typically slower to compute, it is desirable to analyze sufficiently on the faster-computing proxy model to determine a most valuable point at which to compute the physical model, and then to add the physical model computation results for this point to the database of proxy models, and so on iteratively to obtain results with fewer physical model computations. In other embodiments of the present invention, the "models" are all physical models, unless otherwise specified.
The model for calculating L (x; y') on the basis of the physical mechanism as in the first embodiment and the second embodiment is referred to as a "physical model". The proxy model refers to a model in which: set S { (x) based on some results that have been calculated by a physical modeli，Li) I ═ 0, 1, 2,.., m }, where x for each input can be given the value of L (x; y') of the value (noted
) Some forms of the model can additionally give a visual indication of the distance between L (x; y') (note cumulative probability function as Φ (L; x; S)). The computation time is much lower than that of the physical model. Such a model is called a proxy model. Proxy models can be constructed by a variety of methods, such as polynomial fitting, support vector machine regression, gaussian random processes, and the like.
This embodiment takes a proxy model based on the gaussian random process as an example, which gives a gaussian distribution estimate for L with an average of
Standard deviation of
Wherein
LSIs a column vector composed of L in S, F is a polynomial characteristic of x in S, and R is a correlation matrix of x in S and the matrix. f is the polynomial characteristic of x to be estimated, and r is the correlation matrix of x to be estimated and x in S.
Where the polynomial characteristic of x refers to a row vector consisting of several multiplicative forms of the individual components of x. Typically, one expression form of the second-order polynomial characteristic of x ═ γ, F, M ] is
f2nd＝[1 γ F M γ2 F2 M2 FM γM γF]
The correlation matrix is calculated as follows:
a correlation function is determined. The following form of function may be generally chosen:
wherein θ ═ θ1 θ2 … θP]And p ═ p1 p2 … pP]Referred to as a hyper-parameter, may be predetermined or determined with an additional optimization step.
r＝[R(x，x1) R(x，x2) … R(x，xm)]
In the Bayes optimization method, an acquisition function alpha (x; S) is calculated and optimized according to the output result of the proxy model. Commonly used acquisition functions are: estimated values, confidence limits, boosting probabilities, boosting expectations, mean square error, Thompson sampling, and the like. This example is slightly improved based on the probability of boosting, setting the acquisition function as:
α(x；S)＝1-Φ(min(LS)-ε；x；S)
i.e., α (x) is the probability that the proxy model gives L (x; y') greater than the historical minimum of L, or that the absolute value of the difference is less than ε, although it is reduced (to similarly apply the simulated annealing algorithm given in example 1, it is defined that there is no lift or insufficient lift to translate into a minimization problem). Since the term objective function herein means the characterization of the difference between the simulated result and the measured result, and the minimum value should be not less than 0 and close to 0, the following dynamic epsilon can be used:
ε＝λmin(LS)
λ is a small number, and can be 0.1.
In addition, it is also necessary to establish an initial database by calculating a test design, which can be generally used as a Latin hypercube test design, and the method is as follows:
the test protocol was recorded as:
V＝{vij}N×P
each row of the matrix V corresponds to one trial and each column corresponds to one factor. The number of trials is denoted as N, the number of factors (or parameter dimensions) is denoted as P (in this example P ═ 3), and each row of V is a random full permutation of 1, 2. The N discretization levels of each factor have a linear correspondence with the physical parameter values:
where l and u represent the lower and upper limits of the feasible region, respectively, by the amounts indicated in their subscripts (assuming that the feasible region constrains each dimension independently).
Based on the technology, the Bayesian optimization process for estimating the parameters to be determined comprises the following steps:
(1) generating an initial set of x's using Latin hypercube design techniques1，x2，...，xNCalculate their corresponding Li＝L(xi(ii) a y'), thereby forming S { (x)i，Li) I ═ 0, 1, 2, ·, N }. Initializing i-N;
(2) establishing a proxy model based on S, searching x for making the acquisition function alpha (x; S) obtain the minimum value by adopting a global optimization method (such as a simulated annealing algorithm given in embodiment 1), and marking as xi+1；
(3) Calculating Li+1＝L(xi+1(ii) a y'), mixing (x)i+1，Li+1) Adding into S;
(4) let i equal i + 1. And (4) judging whether the termination condition is met, and if not, jumping to (2).
After termination, find x with the smallest L in S, i.e. the estimated value
As shown in fig. 7, the following describes the process of establishing the physical equation of the dry gas seal in examples 1 to 3, and it can be understood that the numerical model in examples 1 to 2 and the physical model in example 3 are physical models described in detail herein. It is not necessarily fully applicable to more general mechanical seals, but there are similarities. It will be appreciated that other reasonable physical models may be created based on actual conditions, such as for purposes of improving computational efficiency and/or computational accuracy.
First, second
The stationary ring 200, which is a compensation ring, is considered to have 3 degrees of freedom: axial displacement zsAngle of inclination gamma about a coordinate axisxsAnd gammays. As shown in FIG. 3, zsWith the direction of the uncompensated ring (in this example, the moving ring 100) pointing towards the compensated ring (in this example, the stationary ring 200) being positive, γxsAnd gammaysPositive in the right-hand helical direction about the X-axis and Y-axis, respectively. Between the stationary ring 200 and the stationary ring 200 seat, which is considered to be a fixed part, there are floating bearings consisting of springs and secondary seals, which are considered to be linear stiffness and damping as a whole, including axial stiffness and damping kzs，czsSum angular stiffness and damping kγs，cγs. Further, the stationary ring 200 has a mass and moment of inertia m, J. Unlike the floating-mounted stationary ring 200, the axial displacement z of the stationary-mounted, axially rotating ring 100rAngle of inclination gamma about a coordinate axisxrAnd gammayrIt is considered to be determined. The relative displacement of the stationary ring 200 and the moving ring 100 and the end face height pattern determine the size of the gap, i.e., the distribution of the film thickness.
Where (x, y) represents the location of a point in the seal zone. The sealing ring band is shown as
The film thickness distribution h is considered as a field defined on the plane. By two specialtiesField magnitude h ofdesAnd hdefIndicating an end face height pattern, the former indicating a pre-designed, normal end face height pattern, such as grooves, pre-designed waviness and conicity, texture, etc., and the latter indicating a non-pre-designed end face height pattern, i.e., waviness and conicity due to deformation. The remaining scalar h0I.e. represents Us-Ur0 and hdes，hdefFilm thickness at zero everywhere. The surface height fluctuations of the roughness scale are not taken into account in the film thickness distribution, but are taken into account in the calculation of lubrication and contact.
With respect to formula (2), U is as previously describedrIs considered to be a determined value. It is generally the case that the end face of the rotating ring 100 is not perpendicular to the axis, but has an inclination angle γrThus having
Secondly, lubrication
The film thickness distribution determines the distribution of the medium pressure on the end face. The medium is considered to be an ideal gas, ignoring its inertia, and is considered to have a constant temperature and viscosity, and the boundary is considered to have no slip. The absolute pressure p of the air film at this time satisfies:
where μ is the dynamic viscosity of the medium, t is the time, and r and θ are the polar diameter and polar angle, respectively, in polar coordinates. pvsCorresponding to the source term of the incoming fluid, typically medium from an orifice (in this case vsThe rate at which it flows into the end face), and v is taken for a seal without this structures＝0。
Is the flow factor. Among the problems that the film thickness is large compared to the surface roughness,
1 can be taken directly. If the problem may include a case where the film thickness is small, the correction should be made, and here:
wherein sigmaaIs the standard deviation of the surface height distribution.
As a performance output, the leak rate q can be calculated by:
three, contact
The film thickness distribution determines the distribution of the contact pressure on the end face. It is important to note that dry gas seals are not in contact during normal, steady operation, but may be in contact under some abnormal conditions, and it is precisely in this case that there is a need to identify the working conditions of the seal, so the model involves considerations of contact (and most of the prior art documents do not).
The contact force is expressed as a pressure pc determined by the thickness h of the local film, using the GW model
The RMS value σ of the acoustic emission signal, as output from the monitorable element, is also given by the analysis of the contact:
this model suggests that the contact action between rough surfaces occurs between the rough peaks, and that each rough peak occurs only at the rough peaksAnd (4) elastically deforming. The problem of contact between two surfaces is equivalently treated as an ideal rigid smooth surface and an equivalent rough surface contact, with the following surface parameters: y is equivalent modulus of elasticity, v is linear velocity of relative motion, ηsR is the roughness peak density (i.e., the number of roughness peaks per unit area)sMean radius of curvature of roughness peak, σsIn order to obtain the standard deviation of the height of the rough peak,
the average height of the roughness peaks (referenced to the average height of the surface), K the transfer coefficient characterizing the propagation and monitoring process, and Φ the probability density function of a standard gaussian distribution.
Dynamics of
The kinetic equation for the stationary ring 200 is set forth below:
here, the inertia of the stationary ring 200 is integrated in the matrix MDamping and stiffness of the bearings are determined by the matrix CsAnd KsA description is given. QfDescribing the force from the gas film, QcThe force from the contact is described.
Equilibrium generalized force QbWith only axial component FbFrom the pressure of the medium surrounding the stationary ring 200 and the pre-tension F in the design statePre：
Various abnormal factors directly acting on the static ring 200, including the error of pretightening force, uneven support, or other various ways of bearing the static ring 200, are simulated by additionally bearing a generalized force Q by the static ring 200. In various embodiments of the present disclosure, the phase of the static ring 200 subjected to the torque is not considered, so that the process is further simplified to
The technical features of the embodiments described above may be arbitrarily combined, and for the sake of brevity, all possible combinations of the technical features in the embodiments described above are not described, but should be considered as being within the scope of the present specification as long as there is no contradiction between the combinations of the technical features.
The above-mentioned embodiments only express several embodiments of the present invention, and the description thereof is more specific and detailed, but not construed as limiting the scope of the invention. It should be noted that, for a person skilled in the art, several variations and modifications can be made without departing from the inventive concept, which falls within the scope of the present invention. Therefore, the protection scope of the present patent shall be subject to the appended claims.
Claims (10)
1. A mechanical seal state analysis method based on a numerical model is characterized by comprising the following steps:
s1, confirming at least one undetermined state parameter to be analyzed in the mechanical seal, and giving a value range of the undetermined state parameter;
s2 obtaining at least one online observation parameter in the mechanical seal;
s3, establishing a numerical model for simulating the working state of the mechanical seal, wherein the numerical model takes the parameter to be determined as an input parameter and outputs a simulation observation parameter;
s4, adjusting the value of the parameter to be determined within the value range by using an optimization method until the difference between the simulated observation parameter and the online observation parameter is within the error allowable range, and recording the current parameter to be determined as the target state parameter.
2. The numerical model-based mechanical seal state analysis method according to claim 1, wherein the mechanical seal device to be analyzed includes a stationary ring and a moving ring, the stationary ring being a compensation ring, the moving ring being a non-compensation ring; the undetermined state parameters comprise dynamic ring end jump and static ring bearing abnormity, and the static ring bearing abnormity is expressed by two components of axial force and deflection moment.
3. The numerical model-based machine seal state analysis method according to claim 2, wherein the online observed parameter includes a fluid average leakage rate, and the online observed parameter further includes a face contact state between the stationary ring and the moving ring.
4. A mechanical seal state analysis method based on a numerical model according to claim 3, characterized in that a rotation speed parameter of a rotating shaft and a pressure parameter upstream and downstream of a mechanical seal end face are measured, and the rotation speed parameter and the pressure parameter are input to the numerical model as determination parameters.
5. The method of analyzing a state of a mechanical seal according to claim 3, wherein a rotational speed sensor is used to measure a rotational speed of the rotating shaft, a pressure sensor is used to measure a pressure parameter upstream and downstream of the end face of the mechanical seal, a flow sensor is used to measure a flow rate of the fluid and calculate an average leakage rate of the fluid, and an acoustic emission sensor is used to detect an end face contact state between the moving ring and the stationary ring.
6. The method for analyzing the state of the mechanical seal based on the numerical model according to any one of claims 1 to 5, wherein a global optimization method is adopted to adjust the value of the parameter to be determined until the difference between the simulated observation parameter and the on-line observation parameter is within an error allowable range; the global optimization method adopts a simulated annealing algorithm, a genetic algorithm and/or a particle swarm algorithm.
7. The numerical model-based mechanical seal state analysis method according to claim 6, wherein the global optimization method includes a first precision optimization model and a second precision optimization model, the second precision optimization model having a higher precision than the first precision optimization model; firstly, a first precision optimization model is adopted for global optimization to obtain a preliminary optimization result, then the preliminary optimization result is used as an initial value, a second precision optimization model is adopted for local optimization to obtain a final optimization result, and a target state parameter is determined.
8. The method for analyzing the state of the mechanical seal based on the numerical model according to any one of claims 1 to 5, wherein the numerical model for simulating the working state of the mechanical seal comprises a physical model and a proxy model, the physical model and the proxy model respectively take the parameter to be determined as an input parameter, and the physical model and the proxy model respectively output a simulation observation parameter; the physical model is built based on the working state of the mechanical seal, the proxy model is built based on the calculation result set of the physical model, and the calculation speed of the proxy model is higher than that of the physical model.
9. A numerical model based mechanical seal state analysis method according to claim 8, wherein the proxy model is established based on the physical model by a polynomial fitting, support vector machine regression, and/or gaussian random process.
10. The numerical model-based mechanical seal state analysis method according to claim 9, wherein the target state parameter is found using a bayesian optimization method.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011453641.4A CN112632849A (en) | 2020-12-12 | 2020-12-12 | Mechanical seal state analysis method based on numerical model |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011453641.4A CN112632849A (en) | 2020-12-12 | 2020-12-12 | Mechanical seal state analysis method based on numerical model |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112632849A true CN112632849A (en) | 2021-04-09 |
Family
ID=75310175
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011453641.4A Pending CN112632849A (en) | 2020-12-12 | 2020-12-12 | Mechanical seal state analysis method based on numerical model |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112632849A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113536637A (en) * | 2021-07-20 | 2021-10-22 | 北京交通大学 | Magnetic liquid sealing pressure resistance analysis method based on MATLAB and COMSOL combined simulation |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120027521A1 (en) * | 2010-07-30 | 2012-02-02 | Fci Holdings Delaware, Inc. | Engineered Mine Seal |
| CN108256160A (en) * | 2017-12-21 | 2018-07-06 | 中国石油天然气集团公司 | A kind of Forecasting Methodology of thermal recovery operating mode Special threading connector compression in sealing contact |
| CN108364350A (en) * | 2018-01-22 | 2018-08-03 | 青岛理工大学 | A kind of Meso-level Structure of Concrete model three-dimensional rebuilding method |
| CN108591451A (en) * | 2018-07-02 | 2018-09-28 | 清华大学 | A kind of mechanical seal intelligent regulating system and method |
| CN109577966A (en) * | 2018-11-29 | 2019-04-05 | 四川富利斯达石油科技发展有限公司 | Using the method for tracer monitoring individual well residual oil saturation |
| CN110083968A (en) * | 2019-05-08 | 2019-08-02 | 中国船舶重工集团公司第七0三研究所 | The compressor characteristics prediction technique of numerical model is influenced based on amendment sealing gland amount of leakage |
| CN110159764A (en) * | 2019-05-31 | 2019-08-23 | 清华大学 | Intelligent mechanical sealing system and its implementation |
| CN111827941A (en) * | 2020-07-07 | 2020-10-27 | 中国石油大学(华东) | Intelligent oil field injection-production real-time optimization and regulation simulation experiment system and method |
-
2020
- 2020-12-12 CN CN202011453641.4A patent/CN112632849A/en active Pending
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120027521A1 (en) * | 2010-07-30 | 2012-02-02 | Fci Holdings Delaware, Inc. | Engineered Mine Seal |
| CN103069110A (en) * | 2010-07-30 | 2013-04-24 | Fci特拉华控股有限公司 | Engineered mine seal |
| CN108256160A (en) * | 2017-12-21 | 2018-07-06 | 中国石油天然气集团公司 | A kind of Forecasting Methodology of thermal recovery operating mode Special threading connector compression in sealing contact |
| CN108364350A (en) * | 2018-01-22 | 2018-08-03 | 青岛理工大学 | A kind of Meso-level Structure of Concrete model three-dimensional rebuilding method |
| CN108591451A (en) * | 2018-07-02 | 2018-09-28 | 清华大学 | A kind of mechanical seal intelligent regulating system and method |
| CN109577966A (en) * | 2018-11-29 | 2019-04-05 | 四川富利斯达石油科技发展有限公司 | Using the method for tracer monitoring individual well residual oil saturation |
| CN110083968A (en) * | 2019-05-08 | 2019-08-02 | 中国船舶重工集团公司第七0三研究所 | The compressor characteristics prediction technique of numerical model is influenced based on amendment sealing gland amount of leakage |
| CN110159764A (en) * | 2019-05-31 | 2019-08-23 | 清华大学 | Intelligent mechanical sealing system and its implementation |
| CN111827941A (en) * | 2020-07-07 | 2020-10-27 | 中国石油大学(华东) | Intelligent oil field injection-production real-time optimization and regulation simulation experiment system and method |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113536637A (en) * | 2021-07-20 | 2021-10-22 | 北京交通大学 | Magnetic liquid sealing pressure resistance analysis method based on MATLAB and COMSOL combined simulation |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Aerens et al. | Force prediction for single point incremental forming deduced from experimental and FEM observations | |
| Liu et al. | Analytical modeling for thermal errors of motorized spindle unit | |
| Lin et al. | Effect of large-area texture/slip surface on journal bearing considering cavitation | |
| Pilch et al. | Measurement of whole-field surface displacements and strain using a genetic algorithm based intelligent image correlation method | |
| Dutta et al. | Tool condition monitoring in turning by applying machine vision | |
| Li et al. | Prediction of surface wear of involute gears based on a modified fractal method | |
| Zhao et al. | Dynamic analysis of spiral-groove rotary seal ring for wet clutches | |
| Choo et al. | The effects of three-dimensional model surface roughness features on lubricant film thickness in EHL contacts | |
| Prajapati et al. | Topography analysis of random anisotropic Gaussian rough surfaces | |
| Liu et al. | Investigation on the influence of interference fit on the static and dynamic characteristics of spindle system | |
| Yang et al. | Wear analysis of angular contact ball bearing in multiple-bearing spindle system subjected to uncertain initial angular misalignment | |
| Chudzik et al. | Effect of radial internal clearance on the fatigue life of the radial cylindrical roller bearing | |
| CN112632849A (en) | Mechanical seal state analysis method based on numerical model | |
| Houpert | An engineering approach to non-hertzian contact elasticity—Part II | |
| Liu et al. | Transmissibility enhanced inverse chatter stability solution | |
| Yao et al. | Position-dependent milling process monitoring and surface roughness prediction for complex thin-walled blade component | |
| Hua et al. | Modeling of lubrication in micro contact | |
| Wang et al. | Orthogonal analysis of multisensor data fusion for improved quality control | |
| Brecher et al. | Method for Calculating Normal Pressure Distribution of High Resolution and Large Contact Area | |
| Sunil Kumar et al. | Influence of flatness and waviness of rough surfaces on surface contact conductance | |
| Zhang et al. | Robust, fractal theory, and FEM-based temperature field analysis for machine tool spindle | |
| Ren et al. | Surface variation modeling by fusing multiresolution spatially nonstationary data under a transfer learning framework | |
| Zhou et al. | A fully coupled 3D elastohydrodynamic model built with MITC element for static performance analysis of gas foil bearings | |
| Wang et al. | Lubrication characteristics of the worn slipper in the slipper-swashplate pair | |
| Yang et al. | Modeling the line contact on an elastic half-space with the statistical approach: Self-affine fractal roughness and numerical framework |
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 |