RESEARCH ARTICLE

Method for solving the nonlinear inverse problem in gas face seal diagnosis based on surrogate models

  • Yuan YIN ,
  • Weifeng HUANG ,
  • Decai LI ,
  • Qiang HE ,
  • Xiangfeng LIU ,
  • Ying LIU
Expand
  • State Key Laboratory of Tribology in Advanced Equipment, Tsinghua University, Beijing 100084, China

Received date: 13 Sep 2021

Accepted date: 07 Apr 2022

Published date: 15 Sep 2022

Copyright

2022 Higher Education Press 2022

Abstract

Physical models carry quantitative and explainable expert knowledge. However, they have not been introduced into gas face seal diagnosis tasks because of the unacceptable computational cost of inferring the input fault parameters for the observed output or solving the inverse problem of the physical model. The presented work develops a surrogate-model-assisted method for solving the nonlinear inverse problem in limited physical model evaluations. The method prepares a small initial database on sites generated with a Latin hypercube design and then performs an iterative routine that benefits from the rapidity of the surrogate models and the reliability of the physical model. The method is validated on simulated and experimental cases. Results demonstrate that the method can effectively identify the parameters that induce the abnormal signal output with limited physical model evaluations. The presented work provides a quantitative, explainable, and feasible approach for identifying the cause of gas face seal contact. It is also applicable to mechanical devices that face similar difficulties.

Cite this article

Yuan YIN , Weifeng HUANG , Decai LI , Qiang HE , Xiangfeng LIU , Ying LIU . Method for solving the nonlinear inverse problem in gas face seal diagnosis based on surrogate models[J]. Frontiers of Mechanical Engineering, 2022 , 17(3) : 33 . DOI: 10.1007/s11465-022-0689-z

1 Introduction

A gas face seal is a type of noncontacting mechanical face seal, such as the type designed to work with a pair of end faces separated by a gas film. However, various faults can occur inside the seal and result in failures [1]. Many studies [27] have reported the use of online monitoring techniques to obtain real-time signals from mechanical face seals (including gas face seals). However, quantitatively determining the cause of an abnormal observation found by the sensor(s) remains challenging.
Physical models for gas face seals have long been applied in design tasks [810]; however, these models are absent in diagnosis tasks, which should have been able to provide quantitative and explainable expert knowledge (unlike machine learning approaches [1115] suitable for systems with no or incomplete physical models). This scenario is mainly attributed to computational reasons. The physics-based diagnosis task can be abstracted as a process of determining the input that results in an output consistent with the signals observed through the physical model or a process of solving an inverse problem. When the model is highly nonlinear and has no explicit inverse expressions, a large number of evaluations are needed to find the correct input (particularly when multiple parameters are undetermined). However, the evaluation of the physical model is costly. Such difficulties are common in systems subject to complex numerical models and do not have explicit equations for the inverse problems.
The gas face seal dynamics model, which directly governs the leakage and friction of a seal, is taken as a typical example. The model discretizes the transient Reynolds equation in space and time and couples it with dynamic equations [16,17]. Solving the model demands considerable time and storage costs. Some existing dynamic models that linearize the transient Reynolds equation for efficiency improvement, namely, semianalytical models [18,19], are effective when applied to conditions close to a particular presupposed state. However, they are not applicable in diagnosis tasks, where the seal state may deviate over an extensive range.
An alternative approach is using surrogate models (metamodels or response surface methodology) [20,21] to replace the physical model for massive iterations by fitting a limited count of existing results. The surrogate model technique has already been widely recognized as an optimal design tool in various fields [2125]. In terms of applications in the field of mechanical face seals (including gas face seals), Liu et al. [26] established a quadratic polynomial response surface model fitting a seal ring vibration model for structure optimization. Additional studies are needed to demonstrate the applicability of surrogate models in mechanical face seals.
In the present work, targeting the task of determining the fault parameters of a gas face seal using the monitored acoustic emission signals is established to solve the inverse problem of the seal dynamic model. This method combines the advantage of the physical models and surrogate models. The method starts with a small initial database and uses rapid surrogate models instead of the original physical model to speculate the input with massively iterative optimizers. It also uses the physical model results as a feedback mechanism to overcome the limited accuracy of the surrogate models with small databases. The method is demonstrated in both simulated and experimental cases.

2 Model of gas face seal dynamics

A typical gas face seal structure is shown in Fig.1. A stator with a planar end face is flexibly mounted on a fixed seat through springs and a secondary seal, and the rotor is rigidly mounted on the rotating shaft. The rotor has spiral grooves on its face (Fig.2). The pressure is high at the outer radius and low at the inner radius. It is sealed by a pair of relatively rotating faces. The following are brief descriptions of the model.
Fig.1 Typical structure of a spiral groove gas face seal.

Full size|PPT slide

Fig.2 Spiral grooves on the rotor face.

Full size|PPT slide

(1) Geometry
As shown in Fig.3, the stator moves with three degrees of freedom, denoted as zs (axial displacement, which increases with separation), γsx (tilt angle around X axis), and γsy (tilt angle around Y axis). Besides the rotation speed ω, the motion of the rotor, denoted as zr, γrx, and γry, is assumed to be unaffected by the forces or moments. Here, zr is regarded as 0, and γrx and γry fluctuate periodically because of the perpendicularity error between the rotor face and the rotation axis. Their values are given by the angles around axes X and Y at t=0 (t refers to the time), denoted as γrx0 and γ ry0, respectively.
Fig.3 Motion of the seal rings.

Full size|PPT slide

Us=[ z s γsxγsy],
Ur=[ z r γrxγry]= [0γrx0cos(ω t)+γry0sin(ω t) γrx0sin(ω t)+γry0cos(ω t)],
h= h0+hgrv+[ 1y x](Us Ur),
where Us and U r denote the motion of the stator and the rotor, respectively, h is the gas film thickness defined on the annulus area Ω= {(x, y)|ri x2+y2 ro}, between the inner radius r i and the outer radius ro, and hgrv denotes the contribution of the grooves (see Fig.2), which values the depth h spr in the groove and 0 otherwise, as
hgrv={ h spr, in the spiral grooves,0,on the land.
(2) Lubrication
We consider the gas to be ideal and inertia free with constant temperature and viscosity and the flow to be laminar and nonslipping at the boundary. Under these conditions, the absolute gas pressure in the film, p, is governed by the transient average flow Reynolds equation [27] with Dirichlet boundary conditions, given by
(φph312μ p)= ω(ph)2 θ+(ph)t ,
{p| r=ri= pi, p |r=ro=po ,
where a point is denoted as ( r,θ) in the polar coordinate system, and is adopted for convenience, μ denotes the gas dynamic viscosity, p i and p o are the boundary pressures at the inner radius ri and outer radius r o, respectively, and φ is a flow factor given by [27]
φ=10.9exp(0.56 hσa),
where σa is the standard deviation of surface height.
(3) Contact
Contact does not occur when the gas face seal operates properly. However, a wide range of faults or low-speed operation may result in contact. Fan et al. [28] stated that light face contact can be considered elastic. Given the summit heights obeying a Gaussian distribution, the contact pressure pc is simulated as [29]
pc= pc(h)= 43Eηs Rs 12 h (wh) 32σ s ϕ(w w¯ s σ s)dw,
where ϕ() denotes the probability density function of the standard normal distribution, and w denotes local surface height. The properties of the roughness of the contacted faces are described with E, ηs, R s, w¯s, and σ s, referring to elastic modulus, asperity density, average asperity radius, mean of asperity height, and standard deviation of asperity height of an equivalent contact pair.
Face contact generates acoustic emissions, which have been considered an effective monitoring approach by researchers [1,46]. The short-time root-mean-square (RMS) is modeled as [28,30]
V=K V0= KΩ 25ωr Eη s h (w h)2σsϕ (w w¯ s σs) dw,
where K is the coefficient describing the transmission from the acoustic emission source (with RMS V0) to the signal sent out by the sensor, assumed to be constant, and V is the RMS value of the signal.
(4) Dynamics
The motion of the flexibly mounted stator is governed by
Ms U¨s+ Cs U˙s +KsUs= Pf+ Pc+ Pb+P,
where the inertia of the stator and the support provided by the springs and secondary seal (modeled as linear damping and stiffness) are considered as
{ Ms= diag(m,J,J), Cs =diag(cs z,csγ, csγ), Ks= diag(k sz,k sγ,k sγ).
The inertia properties are described by the inertia matrix Ms (consisting of mass m, moment of inertia J); and the damping and stiffness are described by the damping matrix Cs (consisting of axial and angular damping c sz and c sγ) and the stiffness matrix Ks (consisting of axial and angular stiffness k sz and k sγ). The damping and stiffness are considered to be potentially abnormal, modeled by being scaled by a scaling coefficient Qs, as
[ks z,ksγ, cs z,csγ]= Qs[k sz0, ksγ 0,cs z0, csγ0].
The generalized forces in Eq. (10) are defined as
Pf =[ ΩpdΩ ΩpydΩ ΩpxdΩ]T ,
Pc =[ Ωpc dΩΩpc ydΩ Ωpc xdΩ]T ,
Pb =[Fb00]T ,
P= [FMx My]T.
Pb has only one component, F b, which comes from the pressure at the back of the stator and the preload. Pb is equal to Pf when UsUr =0 ( Pc=0, as in a noncontacting seal at the designed state), which is how h0 is determined. P contains the faults modeled by the forces or moments exerted on the stator.
The six parameters listed in Tab.1 correspond to the common faults in gas face seal applications. These parameters are considered for identification. The axial compactness change can be represented as an axial force F exerted on the stator. The asymmetry of the stator support is abstracted into the moments Mx and My. The stiffness and damping of the stator support may differ from their design performance because of assembly faults or aging. Thus, they are modeled using a scaling coefficient Qs. The rotor tilt components γ rx0 and γ ry0 commonly occur because of assembly errors.
Tab.1 Suspended causes and their respective ranges
Dimension index j Parameter Lower bound lj Upper bound uj
1 Axial force exerted, F/N −100 100
2 Moment exerted around X, Mx/( Nm) −5 5
3 Moment exerted around Y, My/( Nm) −5 5
4 Support stiffness and damping scale Q s 0 100
5 Rotor tilt around X at t=0, γ rx0/mrad −0.5 0.5
6 Rotor tilt around Y at t=0, γ ry0/mrad −0.5 0.5
The fixed parameters are listed in the Appendix.

3 Establishment of surrogate models

3.1 Test design and initial database

A surrogate model must be established based on an existing database consisting of physical model results. The evaluation of a physical model is expensive. Thus, the evaluation sites must be carefully decided, a process known as test design. Popular test design methods for computer experiments include full factorial design, fractional factorial design, central composite design, Box–Behnken design, Latin hypercube design (LHD), maximum entropy design, minimum integrated mean square error design, maximin/minimax design, and ϕ p-design [20,3133]. The LHD approach has a good space-filling ability and can be implemented conveniently. Thus, it is used in the present work.
For M0 sites with P variables (P=6 here), the LHD scheme can be performed as follows [31].
A matrix D is generated. Each column of which is a random permutation of 0 to ( M0 1).
D= {dij} M0×P.
Given that the parameters are uniformly distributed, the parameters are then randomly generated according to the level indices dij as
vij=lj+dij+ δijM0 (uj lj),i =1,2, , M0, j=1,2, ,P,
where vij contains F,Mx,My, Qs, γrx0, and γry0 for (j=1, 2,...,6) at site i, and uj and l j are the upper and lower bounds, respectively. δ ij is uniformly distributed from 0 to 1.
M0=100 sites are designed using LHD, as summarized in Tab.2. The initial database is generated by evaluating the physical model at these sites.
Tab.2 Sites designed by LHD for the initial dataset
No. F /N Mx/ (Nm) My/ (Nm) Qs γrx0/mrad γry0/mrad
1 −92.0 −2.02 −3.09 82.9 0.147 −0.101
2 −60.8 −4.16 2.95 10.7 0.052 0.498
3 63.5 −3.35 −1.86 12.6 −0.075 −0.231
··· ··· ··· ··· ··· ··· ···
99 17.3 −1.61 3.94 17.7 −0.144 −0.433
100 82.5 4.66 0.90 18.3 −0.290 −0.378

3.2 Kriging model

As a preprocessing step, standardization is first performed to provide the data with a mean of 0 and a variance of 1. This step is applied to each input dimension to give x=[x1 x2 xP]. For the output, V0 is sampled at Q=50 uniformly distributed time steps in a period, and the standardization is applied to all time steps, where V0 is sampled, instead of on each dimension, yielding y=[ y1, y2,, yQ]. The training dataset is denoted as W={( x( k),y(k) )|k=1,2,, M} or in matrix form as X = [ x(1)Tx(2)T x( M) T]T and Y=[ y(1)Ty(2)T y(M)T]T=[ Y1Y2YQ].
Multiple methods for establishing surrogate models have been developed, including polynomial response surfaces, support vector machine regression, and the Kriging model [20]. The Kriging model has been frequently used in recent surrogate model studies because of its good performance for problems with strong linearity and small sample. Our preliminary studies demonstrated its advantages. The Kriging model has several variants; in particular, the universal Kriging model is employed here [34].
The output dimension y j is assumed to be a random variable Y j(x) depending on x, composed of a deterministic polynomial model fq(x) βj and a random variable Zj(x) submitting to a Gaussian process of x as
Yj( x) =fq(x) βj+Zj(x),
where Zj(x) has a mean of 0, a variance of σj2, and a covariance of
Cov (Zj(x),Zj(x))=σj2R(x, x) .
The estimated output y^j and its estimated variance s^j2 are
y^j=fβ ^j+rR1( YjFβ ^j),
s^j2=σj 2(1[fr] [0FT FR]1[ fT rT]),
where R denotes the correlation function, which is used in the correlation matrix of the sample R and the correlation vector for a general input r, F denotes the polynomial model output of the sample, f denotes the polynomial model output for a general input, and βj is the Kriging model coefficent found.
R=[ R (x(1),x(1))R(x(1),x(2))R(x(1), x(M))R(x(2), x(1))R(x(2),x(2))R(x(2), x(M)) R( x(M), x(1))R(x(M),x(2))R(x(M), x(M))],
r=[ R (x,x(1) )R(x,x(2))R(x,x(M))],
F=[ fq(x(1))fq(x(2) )fq(x(M) ) ],
f=fq(x),
β^j= (FT R1F)1 FT R1 Yj,
s^j2=1MYjFβ ^jTR1 Yj Fβ ^j.
A Nugget effect [35] is appended to ensure stability: R1 is calculated as (R+N) 1, where a tiny diagonal matrix N=103I is added before calculating the inverse.
The probability density of y j at a targeted value yj can be estimated as
g( x;X,Yj,yj)=1 2π s^jexp ( (y y^j)22s^j2).
A commonly used correlation function family is the anisotropic Gaussian exponential function, with scaling parameters θ= [ θ1 θ2 θP ] and exponential parameters p=[p1 p2 pP], given by
R( x,x)=exp(1P i=1P( θk || xi xi ||pi)).
Note that the parameters θ and p are consistently applied to all dimensions of y. In the present study, they are restricted to the ranges 1/64θi64 (scaled logarithmically during optimization) and 1pi2, optimized by minimizing [36]
L(θ,p)=M2 j=1Qln s^j 2+Q2ln| R|,
which is realized with an interior point method [37].
A static surrogate model can be established from the initial database. Eight test sites are randomly set for testing according to a uniform distribution, as summarized in Tab.3.
Tab.3 Sites randomly set for testing
No. F /N Mx/ (Nm) My/ (Nm) Qs γrx0/mrad γry0/mrad
1 −56.2 −4.53 1.79 67.9 0.435 −0.116
2 3.9 3.31 −4.65 5.3 0.030 0.171
3 −98.5 −1.17 −4.33 41.7 0.187 0.089
4 86.1 3.46 0.27 9.2 0.154 −0.084
5 40.2 4.10 2.62 26.2 −0.453 0.236
6 −34.4 1.33 2.56 99.1 −0.135 −0.253
7 96.5 2.23 2.53 65.2 −0.427 0.132
8 76.9 −2.27 −0.64 76.6 −0.022 −0.262
The discrepancy measurement d is defined as
d( y,y)=j= 1Q(yjyj)2j= 1Q(yj)2+ j=1Q (yj)2.
A comparison of the true output at the test sites and the output estimated by the surrogate model is presented in Fig.4(a), along with the discrepancy defined above. The surrogate model reflects the overall trend for several cases but does not perform well in others. Although one can eventually find a database large enough to generate a surrogate model with satisfactory accuracy, this approach is against the initial purpose of decreasing the number of physical model evaluations. Given the small sample size (100 observations for six dimensions) and the strong nonlinearity of the problem, the result can be regarded as successful in terms of the sparse training set.
Fig.4 Estimation of the surrogate model established on the initial database: (a) comparison of the true output and the estimated output; (b) density of the scaled error.

Full size|PPT slide

Although the accuracy is limited by the database size, the estimation for the distribution provides information about the uncertainty, which can be used to support decision-making. The 95% confidence range is plotted in Fig.4(a). The ratio of the error to the estimated standard deviation (( y^j yj)/(y^ jyj)s^ js^j) is plotted in Fig.4(b), along with the probability density function of the standard normal distribution, which it should obey. The consistency between the two demonstrates that the probability estimation is reasonable.
We establish a probability-based criterion for determining the sites for physical model evaluations by taking advantage of this feature. Thus, an efficient method is presented by dynamically updating the surrogate model while searching for the desired input using the feedback acquired from the physical model.

4 Probability-based iterative method for solving the inverse problem

The present study aims to determine the causes of the abnormal signals. This inverse problem can be abstracted as a process of searching for an identified x^ of the input x for an observed output y, or
y=y( x) y( x^),
where y(x) denotes the physical model of the gas face seal, which fetches the input parameter x and returns the standardized fluctuation of V 0 as output.
The genetic algorithm is considered because the mathematical nature of y(x) (i.e., multimodality, strong nonlinearity) is unclear, and its derivative cannot be analytically deduced [38,39]. Genetic algorithms can handle problems with weak mathematical conditions and undefined derivatives. However, they require a large number of iterations. Thus, they cannot be directly applied to the physical model. Enlightened by the idea of Bayesian optimization [40], we employed the surrogate models to construct the objective function for the genetic algorithm, which can be evaluated rapidly.
Therefore, a probability-based iterative method is established. Based on all available results from the physical model, the next point for evaluation is determined by
x=arg maxxSG( x)=arg maxxSj =1Qg(x;X,Yj,yj) ,
where S denotes the feasible region when searching for x, as listed in Tab.1.
The function G (x) 1/(PQ) is used as the fitness function for the genetic algorithm, where the power 1/ 1(PQ) (PQ) is applied to prevent an extreme product of probability densities. Compared with using the discrepancy definition in Eq. (32) for the optimization, using this probabilistic expression encourages the algorithm to be sensitive to the error on sites with low randomness, allowing early escape from the local minimum.
Then, the physical model is evaluated at x. The result acts as validation because the surrogate model may not always be sufficiently accurate, particularly when a limited database is used. If the iteration continues, the result is added to the database, thereby updating the surrogate model near the site that it previously identified.
The whole method proceeds in the following cycle (illustrated in Fig.5):
Fig.5 Schematic flowchart of the inverse problem solution procedure.

Full size|PPT slide

(1) Prepare the initial database W 0, as described in Section 3.1. The database is denoted by W, which is dynamically updated.
(2) Terminate and return the ( x(i),y(i)) in W with minimum d (y(i), y) as the identification if (x(i),y(i))W:d (y(i), y) d T or the iteration limit kT has been exceeded; otherwise, proceed. Here, discrepancy termination threshold dT and maximum iteration count kT are set to be 0.02 and 100, respectively.
(3) Establish a surrogate model on the database W (as described in Section 3.2).
(4) Find the x by optimizing Eq. (34) using the genetic algorithm.
(5) Evaluate y=y(x) at the maximum point found in step 4 and add the newly obtained (x,y) into the database W. Then, go back to step (2).
The infeasible cost of evaluating the physical model numerous times is avoided using the above method, with the surrogate models used instead. Additionally, the relatively low accuracy of the surrogate models with the limited database is overcome using the reliable physical model results as feedback. Therefore, the method combines the rapid evaluation speed of the surrogate models with the accuracy of the physical model.

5 Results and discussion

5.1 Validation of simulated cases

The simulated outputs of the sites listed in Tab.3 were provided to test the proposed method by identifying the axial force F, moment M=Mx2+My2, support stiffness and damping scale Q s, and rotor tilt γ r=γrx02 +γry02. The moments around the two axes were combined and the tilts around the two axes were also merged because identifying their phases does not make sense in practice.
The identification results of all cases are presented in Tab.4 and Fig.6. For cases 1, 2, 3, 5, and 6, early terminations occurred because highly similar outputs (with discrepancy less than d T=0.02) were found. For the other cases, the iterations were terminated because the maximum iteration count (k T=100) was reached. Moreover, the x corresponding to the minimum discrepancy found was returned. Fig.7 exhibits the iteration process toward the target y in each case. Each iteration costs approximately 5.5 min on average, within which approximately 4.5 min is the cost for solving the physical model with an 8 × 3 GHz computer.
Tab.4 Validation of the simulated cases by comparing the true and identified parameters
No. Parameter type F /N M /( Nm) Qs γr /mrad Qs γr /mrad Iteration count Discrepancy
1 Truth −56.2 4.87 67.9 0.450 30.6 34 0.019
Identified −71.9 4.36 69.0 0.462 31.9
Error −15.7 −0.51 1.1 0.012 1.4
2 Truth 3.9 5.71 5.3 0.174 0.9 8 0.009
Identified −92.8 4.20 39.2 0.034 1.3
Error −96.7 −1.51 33.8 −0.139 0.4
3 Truth −98.5 4.49 41.7 0.207 8.6 9 0.006
Identified −92.9 4.68 85.5 0.105 9.0
Error 5.5 0.20 43.7 −0.102 0.3
4 Truth 86.1 3.47 9.2 0.175 1.6 43/100 0.045
Identified −3.8 2.72 15.1 0.158 2.4
Error −89.9 −0.75 5.9 −0.017 0.8
5 Truth 40.2 4.87 26.2 0.510 13.4 7 0.017
Identified 60.4 4.91 80.0 0.160 12.8
Error 20.2 0.04 53.8 −0.350 −0.6
6 Truth −34.4 2.89 99.1 0.287 28.4 35 0.019
Identified −31.0 2.58 68.5 0.451 35.9
Error 3.3 −0.31 −30.6 0.165 2.5
7 Truth 96.5 3.37 65.2 0.447 29.1 14/100 0.040
Identified 96.9 2.82 81.8 0.439 35.9
Error 0.4 −0.55 16.6 −0.008 6.7
8 Truth 76.9 2.36 76.6 0.263 20.2 56/100 0.048
Identified 87.6 2.49 68.1 0.278 19.0
Error 10.7 0.13 −8.5 0.015 −1.2
Fig.6 Actual values of the queried parameters and their respective identified values. M and Qsγr are identified well. Qs and γr are difficult to be discriminated because of physical reasons.

Full size|PPT slide

Fig.7 Iteration process toward the targeted y′ for the simulated cases. The dashed lines denote the best in the initial database and the improving steps. The black and red solid lines denote the searched target and the result found, respectively. (a)–(h) correspond to cases 1–8.

Full size|PPT slide

The method shows the following performance:
(1) The identifications for M are generally close to the respective ground truths. Except for case 2, where the absolute error reaches 1.51 N∙m, all cases exhibit an absolute error of M less than 0.8 N∙m .
(2) The identifications for F fall into great errors in cases 2 and 4. This scenario is mainly attributed to a weakness in following the high-order oscillation features in contact signals.
(3) The identifications for Q s and γ r show considerable errors. However, the existing studies on seal dynamics concluded [16,41] that the support stiffness and damping affect the relative motion of seal faces by impeding the tracing ability, thereby scaling the effect of rotor tilt. This finding means that the effects of Qs and γr on the contact are highly similar. The analysis of Q s γr is appended in Tab.4 and Fig.6, showing relatively small errors. Therefore, the weakness in recognizing Q s and γ r is due to their physical confounding.
(4) For five of the eight cases, the iteration is terminated to find a highly compatible match (three of which have less than 10 iterations). For the other cases, the requirement of early termination is not met; however, outputs with discrepancies of less than 0.05 are found within the limit of 100 iterations. This finding shows a balance between opportunism and robustness.
In summary, the validation using simulated cases shows that the method has good performance in identifying M and Qs γr in the task of rooting the causes of face contact in limited iterations.

5.2 Validation of experimental cases

Experimental data were employed to validate the performance of the proposed method. In our previous study [4], acoustic emission signals were obtained from a gas face seal running at 600 r/min and 0.4 MPa (see test rig in Fig.8(a)). The signals were processed with a band-pass filter (centered at 160 kHz with 40 kHz width, which successfully filtered out the noises), and the RMS of each short-time wave was calculated.
Fig.8 (a) Test rig for exerting abnormities and online monitoring and (b) filtered acoustic emission RMS during startup.

Full size|PPT slide

Unlike validation on simulated cases, one cannot rule out unintentional errors in experimental processes. The strongest starting process data (found at 60 r/min, see Fig.8(b)) were compared with the simulated low-speed operation results to determine K in Eq. (9) because this moment is when the contact is the strongest (the RMS is not that high because of the low speed). Thus, it is the point that is least affected by unintentional errors. The simulated V0 here is 1.28 W0.5. Thus, K is estimated to be 21.1 mV W0.5.
A loading mechanism appended to the test rig exerted a force at an eccentric point on the stator, as shown in Fig.8(a). In cases 1, 2, and 3, the axial forces of 29.7, 54.6, and 79.5 N were exerted, respectively, at a point 51 mm eccentric from the axis, thereby inducing a moment. These are the only known components. The unintentional components of F and M, as well as the value of Q s and γ r, are inaccessible. However, they are assumed to be minor factors here.
The results are presented in Tab.5 and Fig.6 and Fig.9. From these results, we can state the following:
Tab.5 Validation of the experimental cases by comparing the exerted and identified parameters and checking the consistency of parameters not changed through the experiments
No. Parameter type F /N M /( Nm) Qs γr /mrad Iteration count Discrepancy
1 Exerted −29.7 1.51 88/100 0.195
Identified −37.5 2.16 6.0
Error −7.8 0.65
2 Exerted −54.6 2.78 7/100 0.088
Identified −85.7 3.72 12.1
Error −31.1 0.94
3 Exerted −79.5 4.05 65/100 0.047
Identified −87.5 5.01 5.3
Error −8.0 0.96
Fig.9 Experimental data and results found by the method for each case. The dashed lines denote the initial and improving steps; the black and red solid lines denote the searched target and the result found, respectively. (a)–(c) correspond to cases 1–3.

Full size|PPT slide

(1) The trend of the identified F and M is consistent with the exerted values. The systematic bias in M may be caused by unintentional factors in the experiments that aggravate the contact.
(2) The spire in the experimental observation of case 1 fails to patch with physical model results, represented by high discrepancy (0.195). This phenomenon probably comes from the face deformation and no longer appears explicitly in cases 2 and 3, where the intentionally controlled factors have great dominance.
(3) The values of Q sγ r are approximate across the three cases (compared with its range of 0 to 70.7 mrad). This finding demonstrates the effectiveness of the identification by obtaining consistent results.
The proposed method was proved effective in real-life experimental cases on the basis of the above results. Given the advantage of using the physical model as a basis, one can reconstruct the detailed transient motion and lubrication states of the seal, as shown in Fig.10. This approach is useful for comprehensively explaining the result or establishing a digital twin system.
Fig.10 Reconstructing the detailed relative motion and the separation and pressure distributions by taking experimental case 2 as an example.

Full size|PPT slide

6 Conclusions

(1) A surrogate-model-assisted method is developed to solve the computational difficulty of the inverse problem abstracted from physics-based diagnosis tasks of gas face seals by efficiently inferring the input for a targeted output. The method starts with a small initial database. It performs iterations that take advantage of the rapidity of the surrogate models and the reliability of the physical model.
(2) For the establishment of the surrogate models of the signal output, the Kriging model is validated as capable of using small multivariable databases to provide rapid estimations of the output distribution instead of merely an expectation. Thus, it provides good support in determining the optimal evaluation sites of the physical model while an extensive database is inaccessible.
(3) The rapid surrogate models are employed for the construction of objective functions for large-scale iterations to find optimized sites for physical model evaluations, thus overcoming the infeasibility of applying massively iterative search algorithms to the physical model. The physical model is evaluated at the selected sites to provide detailed feedback to the surrogate models and overcome the limited accuracy of the surrogate models on small databases. The physical model is evaluated only a limited number of times using such an iteration routine. However, the expert knowledge in the physical model is delivered by the surrogate models, and the reliability of the final result is ensured.
(4) The proposed method is validated on simulated and experimental cases of gas face seal dynamics. The results demonstrate that the method can effectively identify the fault parameters inducing the signal output in limited physical model evaluations.
The work presented in this paper provides a quantitative, explainable, and feasible approach for identifying the cause of gas face seal contact. The method is also suitable for diagnosis tasks in mechanical devices that face similar difficulties.

Acknowledgements

This work was supported by the National Key R&D Program of China (Grant No. 2020YFB2010000) and the National Natural Science Foundation of China (Grant No. U1737209). None of the funding bodies influenced the study at any stage.
1
Fan Y E, Gu F S, Ball A. A review of the condition monitoring of mechanical seals. In: Proceedings of ASME the 7th Biennial Conference on Engineering Systems Design & Analysis. Manchester, 2004, 3: 179–184

DOI

2
Phillips R L, Jacobs L E, Merati P. Experimental determination of the thermal characteristics of a mechanical seal and its operating environment. Tribology Transactions, 1997, 40(4): 559–568

DOI

3
DiRussoE. Film Thickness Measurement for Spiral Groove and Rayleigh Step Lift Pad Self-Acting Face Seals. NASA TP-2058, 1982

4
Huang W F, Lin Y B, Liu Y, Liu X F, Gao Z, Wang Y M. Face rub-impact monitoring of a dry gas seal using acoustic emission. Tribology Letters, 2013, 52(2): 253–259

DOI

5
Towsyfyan H, Gu F S, Ball A D, Liang B. Tribological behaviour diagnostic and fault detection of mechanical seals based on acoustic emission measurements. Friction, 2019, 7(6): 572–586

DOI

6
Daraz A, Alabied S, Zhen D, Gu F S, Ball A D. Detection and diagnosis of mechanical seal faults in centrifugal pumps based on acoustic measurement. In: Ball A, Gelman L, Rao B K N, eds. Advances in Asset Management and Condition Monitoring. Cham: Springer, 2020, 963–975

DOI

7
Medina-Arenas M, Sopp F, Stolle J, Schley M, Kamieth R, Wassermann F. Measurement and analysis of inadequate friction mechanisms in liquid-buffered mechanical seals utilizing acoustic emission technique. Vibration, 2021, 4(1): 263–283

DOI

8
Feldman Y, Kligerman Y, Etsion I. Stiffness and efficiency optimization of a hydrostatic laser surface textured gas seal. Journal of Tribology, 2007, 129(2): 407–410

DOI

9
Blasiak S, Zahorulko A V. A parametric and dynamic analysis of non-contacting gas face seals with modified surfaces. Tribology International, 2016, 94: 126–137

DOI

10
Jiang J B, Zhao W J, Peng X D, Li J Y. A novel design for discrete surface texture on gas face seals based on a superposed groove model. Tribology International, 2020, 147: 106269

DOI

11
Zhang Z, Li X H. Acoustic emission monitoring for film thickness of mechanical seals based on feature dimension reduction and cascaded decision. In: Proceedings of 2014 the 6th International Conference on Measuring Technology and Mechatronics Automation. Zhangjiajie: IEEE, 2014, 64–70

DOI

12
YinY, LiuX F, HuangW F, Liu Y, HuS T. Gas face seal status estimation based on acoustic emission monitoring and support vector machine regression. Advances in Mechanical Engineering, 2020, 12(5): 168781402092132

13
ZhangZ H, Min F, ChenG S, ShenS P, WenZ C, ZhouX B. Tri-partition state alphabet-based sequential pattern for multivariate time series. Cognitive Computation, 2021 (in press)

14
Li T Y, Qian Z J, Deng W, Zhang D Z, Lu H H, Wang S H. Forecasting crude oil prices based on variational mode decomposition and random sparse Bayesian learning. Applied Soft Computing, 2021, 113: 108032

DOI

15
Ran X J, Zhou X B, Lei M, Tepsan W, Deng W. A novel k-means clustering algorithm with a noise algorithm for capturing urban hotspots. Applied Sciences, 2021, 11(23): 11202

DOI

16
Miller B A, Green I. Numerical formulation for the dynamic analysis of spiral-grooved gas face seal. Journal of Tribology, 2001, 123(2): 395–403

DOI

17
Green I. A transient dynamic analysis of mechanical seals including asperity contact and face deformation. Tribology Transactions, 2002, 45(3): 284–293

DOI

18
Zirkelback N, San Andre’s L. Effect of frequency excitation on force coefficients of spiral groove gas seals. Journal of Tribology, 1999, 121(4): 853–861

DOI

19
Miller B, Green I. Constitutive equations and the correspondence principle for the dynamics of gas lubricated triboelements. Journal of Tribology, 1998, 120(2): 345–352

DOI

20
Simpson T W, Peplinski J D, Koch P H, Allen J K. Metamodels for computer-based engineering design: survey and recommendations. Engineering with Computers, 2001, 17(2): 129–150

DOI

21
Bhosekar A, Ierapetritou M. Advances in surrogate based modeling, feasibility analysis, and optimization: a review. Computers & Chemical Engineering, 2018, 108: 250–267

DOI

22
Han Z H. Kriging surrogate model and its application to design optimization: a review of recent progress. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3197–3225

DOI

23
Kontogiannis S G, Savill M A. A generalized methodology for multidisciplinary design optimization using surrogate modelling and multifidelity analysis. Optimization and Engineering, 2020, 21(3): 723–759

DOI

24
Garg A, Liu C, Jishnu A K, Gao L, Le Phung M L, Tran V M. A Thompson sampling efficient multi-objective optimization algorithm (TSEMO) for lithium-ion battery liquid-cooled thermal management system: study of hydrodynamic, thermodynamic, and structural performance. Journal of Electrochemical Energy Conversion and Storage, 2021, 18(2): 021009

DOI

25
Tang H S, Ren Y, Kumar A. Optimization tool based on multi-objective adaptive surrogate modeling for surface texture design of slipper bearing in axial piston pump. Alexandria Engineering Journal, 2021, 60(5): 4483–4503

DOI

26
Liu Y, Liu Q X, Yin M, Yin G F. Dynamic analysis and structure optimization of a floating ring system in dry gas seal. Journal of Advanced Mechanical Design, Systems and Manufacturing, 2018, 12(7): JAMDSM0128

DOI

27
Patir N, Cheng H S. An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication. Journal of Tribology, 1978, 100(1): 12–17

DOI

28
FanY B, Gu F S, BallA. Modeling acoustic emissions generated by sliding friction. Wear, 2010, 268(5–6): 811–815

29
Greenwood J A, Williamson J B P. Contact of nominally flat surfaces. Proceedings of the Royal Society of London, 1966, 295(1442): 300–319

DOI

30
Yin Y, Hu S T, Huang W F, Liu X F, Liu Y, Wang Z X. A bi-Gaussian acoustic emission model for sliding friction. IOP Conference Series: Material Science and Engineering, 2019, 686(1): 012026

DOI

31
Mckay M D, Beckman R J, Conover W J. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 1979, 21: 239–245

DOI

32
Morris M D, Mitchell T J. Exploratory designs for computational experiments. Journal of Statistical Planning and Inference, 1995, 43(3): 381–402

DOI

33
Lee B W, Peterson J J, Yin K H, Stockdale G S, Liu Y C, O’Brien A. System model development and computer experiments for continuous API manufacturing. Chemical Engineering Research & Design, 2020, 156: 495–506

DOI

34
Sacks J, Welch W J, Mitchell T J, Wynn H P. Design and analysis of computer experiments. Statistical Science, 1989, 4(4): 409–423

DOI

35
Yin J, Ng S H, Ng K M. Kriging metamodel with modified nugget-effect: the heteroscedastic variance case. Computers & Industrial Engineering, 2011, 61(3): 760–777

DOI

36
Toal D J J, Bressloff N W, Keane A J. Kriging hyperparameter tuning strategies. AIAA Journal, 2008, 46(5): 1240–1252

DOI

37
Byrd R H, Gilbert J C, Nocedal J. A trust region method based on interior point techniques for nonlinear programming. Mathematical Programming, 2000, 89(1): 149–185

DOI

38
Herrera F, Lozano M, Verdegay J L. Tackling real-coded genetic algorithms: operators and tools for behavioural analysis. Artificial Intelligence Review, 1998, 12(4): 265–319

DOI

39
Deng W, Zhang X X, Zhou Y Q, Liu Y, Zhou X B, Chen H L, Zhao H M. An enhanced fast non-dominated solution sorting genetic algorithm for multi-objective problems. Information Sciences, 2022, 585: 441–453

DOI

40
Jones D R, Schonlau M, Welch W J. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 1998, 13(4): 455–492

DOI

41
Hu S T, Huang W F, Liu X F, Wang Y M. Influence analysis of secondary O-ring seals in dynamic behavior of spiral groove gas face seals. Chinese Journal of Mechanical Engineering, 2016, 29(3): 507–514

DOI

Outlines

/