Clamp Nonlinear Modeling and Hysteresis Model Parameter Identification

Based on the Bouc-Wen model, a nonlinear hysteretic restoring force model is established with dynamic equations. The hardening and softening of the material after reaching the yield limit are described by the nonlinear restoring force-displacement hysteretic curve. The Gamultiobj algorithm and the group search optimization (GSO) algorithm are used to fit and identify the material softening and hardening hysteresis curves, respectively. The effectiveness of the optimization algorithm for identifying the parameters of the hysteretic model is verified, and the optimization effects of the two algorithms are compared. The results show that the Gamultiobj algorithm has better parameter identification ability and curve fitting ability for the hysteretic model. Then, the hysteresis curve describing the nonlinear characteristics of a clamp is obtained through the static stiffness experiment of the clamp. The experimental curve is fitted by the Gamultiobj algorithm. As a result, the nonlinear Bouc-Wen model parameters of the clamp are obtained, and the stiffness and damping of the clamp are recognized. The distribution statistics of the obtained parameters are performed, and it is found that each parameter satisfies a certain probability distribution, which indicates that the parameter identification result is reasonable.


I. INTRODUCTION
The external piping system is an important part of the aero engine, and it plays an important role in transporting energy. The clamps support and fix the pipeline at a specific position, which plays a key role of frequency modulation and shock absorption [1]. The material of clamp is metal rubber, which is a nonlinear material with large damping. At present, when modeling the clamp-piping system, the clamps are mostly simplified into linear springs and dampers, which will affect the dynamic response analysis of the system. Therefore, the establishment of a clamp nonlinear model is of great significance.
Many people have done research on the mechanical characteristics analysis and modeling of clamps. Bezborodov et al. [2] established a solid model of the clamp pipeline using finite element software considering the influence of pipe joints, calculated the natural frequency of the The associate editor coordinating the review of this manuscript and approving it for publication was Jingang Jiang . pipeline, and verified the effectiveness of the clamp pipeline model through experiments.
The modeling method of clamp-piping system has been relatively mature. Gao et al. [3] simplified the stiffness of the clamp to spring support and proposed a method to reduce the dimension of the stiffness matrix based on finite element method, which was used to establish a liquid-filled pipeline model, and calculate the natural frequency of the pipeline. Liu et al. [4] discretized the clamps into springs and dampers, established a straight multi-clamp fluid-filled pipeline model, and performed dynamic response analysis. Chai et al. [5] based on the static stiffness experiment of the clamp, discretized the clamp into springs and dampers, and established the single-clamp and dual-clamp pipeline models by finite element method, calculated the natural frequency and dynamic response. The results are verified through experiments.
Regarding the high energy dissipation capacity, large damping and nonlinearity of metal rubber materials, many scholars have tried to establish mathematical models to describe them. Aiming at the hysteresis friction phenomenon VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ of metal rubber, Iwan [6] regarded the dry friction surface as a series combination of a spring and an ideal Coulomb friction pair, and proposed a bifold line hysteretic restoring force model. Zhang et al. [7] conducted theoretical research on metal rubber materials, and used semi-constitutive models of nonlinear elastic stiffness, nonlinear viscous damping and bilinear Coulomb damping to perform theoretical modeling and analysis on metal rubber material devices. Yu et al. [8] conducted comparative tests on two shock absorbers of metal rubber and rubber materials, and studied the deformation under different loading rates, different loading forces, quasistatic loading and other different conditions. To complete the modeling and solving of a large class of hysteretic systems under random excitation, Wen [9], [10] extended the nonlinear model proposed by Bouc [11] in 1976, and obtained the Bouc-Wen model. In 1986, a biaxial hysteretic restoring force model was proposed based on the Bouc-Wen model. The response of the hysteretic system under random earthquakes was studied, and the validity of the model was verified by the bidirectional force test of reinforced concrete columns. At present, the Bouc-Wen model has been widely used in the modeling of nonlinear hysteretic systems. Ok et al. [12] improved the Bouc-Wen model and got a more accurate rubber casing model. The VisualDOC program is used to find the optimal parameters of the bushing model to meet the static and dynamic characteristics caused by the sine and step excitation input. The results show that the model can predict the bushing force under sine, step and random input well. Yan et al. [13] used the Monte Carlo method to study the static response of the bilinear and Bouc-Wen model under the random excitation of Poisson white noise. Li et al. [14] used the stiffness degradation of the Bouc-Wen model to simulate the hysteretic characteristics of the tower structure, established a nonlinear equation for the connection structure, and studied the dynamic reliability of the connection structure at different locations using the pseudo-excitation method. Ding et al. [15] based on the background of atomic force microscope (AFM) driven by piezoelectric actuators (PEAs) proposes a sliding mode control coupled with an inverse Bouc-Wen hysteresis compensator to improve the positioning performance of PEAs. The effectiveness and feasibility of the proposed method are verified by experiments.
Regarding the parameter identification of Bouc-Wen model, scholars at home and abroad have conducted a lot of research. Peng [16] studied that the Bouc-Wen model cannot simulate the hysteresis loop with force hysteresis, and proposed an improved Bouc-Wen model BWBN to simulate the hysteresis loop. The parameters of the revised BWBN model were identified using MATLAB-SIMULINK design and optimization toolbox, and good results were achieved. Aiming at the problem that the Bouc-Wen model cannot reflects the hysteresis loop accurately at low speeds, lacks hysteresis characteristics at high speeds, and cannot express the deviation of the accumulator damping force in a single rod damper, Wang [17] proposed an improved Bouc -Wen model. The parameters of the new model were identified with experimental data, and finally the validity of the revised Bouc-Wen model was verified. Talatahari et al. [18] proposed an adaptive charged system search (CSS) optimization method to identify the parameters of Bouc-Wen model. This method has great robustness and accuracy, and can be successfully applied to strong nonlinear recognition problems.
How to choose the parameter identification process and set the optimization algorithm of Bouc-Wen model are always mentioned. Charalampakis and Dimou [19] uses two particle swarm optimization (PSO) algorithms to identify Bouc-Wen systems. This algorithm is used to identify a bolt welding hysteresis system that represents a full-size bolt-welded steel connection. The performance of the algorithms are compared with other evolutionary algorithms on highly nonlinear recognition problems with different computational budget levels. Ye and Wang [20] used the Bouc-Wen model to describe the softening and hardening of the material after reaching the yield limit, and used the improved PSO algorithm proposed to effectively complete the parameter identification of the hysteretic model. Ding et al. [21] used PSO method to identify the parameters of piezoelectric actuator based on Bouc-Wen hysteretic model, and also used it to adjust the PID parameters to improve the tracking accuracy of piezoelectric actuator at different input frequencies.
According to the above analysis, the clamp model is mostly established through complex finite element solid models. The mechanical properties of the clamp are mostly determined through experiments. In the process of establishing the clamp-piping model, the clamps are mostly simplified into spring supports and dampers, ignoring the nonlinearity of the clamps. In terms of the application of the Bouc-Wen model, few researchers apply the Bouc-Wen model to modeling the nonlinearity of clamps.
In order to describe the nonlinearity of the clamp more simply and effectively, this paper establishes a nonlinear hysteresis model based on the Bouc-Wen model. A multiobjective genetic algorithm is used to identify the parameters of the material softening and hardening hysteresis models. The effectiveness of the optimization algorithm for parameter identification of nonlinear models is verified. Then, based on the hysteresis curve obtained by the clamp static stiffness experiment, the proposed nonlinear hysteresis model is optimized to match the experimental curve. The results proved the applicability of the proposed model to the nonlinearity of the clamp, realized the parametric modeling of the clamp, and completed the identification of the stiffness and damping of the clamp.

II. NONLINEAR DYNAMIC MODEL OF CLAMP A. BOUC-WEN HYSTERESIS MODEL
Bouc-Wen hysteresis model is a set of nonlinear equations defined by [22] mẍ where, u (t) = B sin(ωt) which means the acceleration of the applied load. The parameters of Bouc-Wen model meet the following criteria [23] Define the following parameters [24] At this time, Eq. (1) can be rewritten as Bouc-Wen model is widely used to describe nonlinear hysteretic systems. By choosing appropriate parameters, the model can represent many different shapes of hysteresis loops (for example, Figure 1). Therefore, when using the Bouc-Wen model to describe the nonlinear hysteretic system, the parameter identification of the hysteretic model is an indispensable part.

B. PARAMETER IDENTIFICATION METHODS FOR BOUC-WEN MODEL
For the Bouc-Wen model in Eq. (4), the restoring force is defined as F (t) = k f x (t)+k w w(t) [24]. Given a set of appropriate parameters, the hysteresis loops shown in Figure 1 can be obtained.
When identifying the parameters of the hysteresis model, set the model parameters as 0 = [ρ 0 , σ 0 , n 0 , k f0 , k w0 , c 0 ], and the identified parameters are For the same input u(t), substituting and 0 will obtain different restoring force F and displacement x. Define the root mean square error [20] (RMSE) as The term 1/F 0max (1/x 0max ) in Eq. (6) is to eliminate the dimension. Take Eq. (6) as the objective function, and find its minimum value through optimization algorithm. In optimization algorithm, when the value of the objective function decreases, the corresponding parameter will gradually tend to 0 .

C. IMPROVED GENETIC ALGORITHM-GAMULTIOBJ ALGORITHM
The Gamultiobj function has some basic concepts [27]

1) DOMINATE AND NON-INFERIOR
When using the multi-objective optimization algorithm solves the minimum value of objective function, assuming that at least one target value of individual p is smaller than individual q, and all the target values of the individual p are not greater than the individual q, it is called p dominates q, or p is not inferior to q.

2) ORDER VALUE AND FRONT END
If individual p dominates individual q, the order value of p is considered to be lower than q. If the ordinal values of two individuals are the same, they are considered to be non-inferior to each other. Individuals with an ordinal value of 1 are at the first front end, and individuals with an ordinal value of 2 are at the second front end. Individuals at the first front end are not dominated, and individuals at the second front end are dominated by the first front end.

3) CROWDED DISTANCE
Only individuals within the same front end can calculate the crowded distance, which characterizes the degree of crowdedness between individuals. The large crowding distance indicates that the diversity of the population members is good.

4) OPTIMAL FRONT-END INDIVIDUAL COEFFICIENT
The optimal front-end individual coefficient refers to the ratio of the number of individuals in the first front-end to the number of members of the entire population, expressed by ParetoFraction, with a value between 0 and 1.
The idea of the algorithm is to initialize the members according to the set number of members of the entire population and the value range, and then allocate the members to each front end according to the order value. Each iteration will evolve the current members of the entire population to the next generation, and unqualified members will be eliminated. There are two criteria for the elimination of entire population members. One is that members with a high order value are eliminated. The other is that entire population members with the same order value and small crowding distance VOLUME 9, 2021 are eliminated to ensure the diversity of entire population members. When the iteration reaches the maximum, the algorithm stops running and returns the parameter value and the corresponding objective function value. The number of members returned is determined by Min {ParetoFraction × Number of members of the population, number of individuals in the first front end} (7) D. GAMULATIOBJ ALGORITHM The Gamultiobj optimization algorithm is used to identify the parameters of the hysteretic model in Figure 1. The feasibility of the optimization algorithm is evaluated through the shape matching of the hysteretic model, the restoring force and displacement error, the calculation efficiency and the accuracy of parameter identification. In order to perform numerical simulation on the Bouc-Wen model, the fourth-order Runge-Kutta method is used to solve Eq. (4), and the total time length is set to t = 4s, and the step length dt = 0.002s. Take the last four cycles of the time-domain displacement response as the prediction data, and the data length L = 201. In order to facilitate the observation of the convergence of the objective function, take the function value of Eq. (6) after each iteration and define  Figure 3 represents the relative error of displacement, and the ordinate represents the relative error of restoring force. The maximum value of the restoring force of hysteresis loop 1 is 83.89N, and the relative error is less than 0.6N. The maximum value of restoring force of hysteresis loop 2 is 35.6N, and the relative error is less than 0.5N, and the matching effect is perfect. Figure 4 shows the convergence of the RMSE of the Gamultiobj algorithm.
For the group search optimization (GSO) algorithm, the Refs. [28]- [30] has a detailed introduction. Setting the same number of population members, the number iterative and parameter value range for the GSO algorithm can compare   the parameter identification capabilities of the two optimization algorithms for the Bouc-Wen model. Therefore, for the hysteresis model of softening and hardening in Figure 1, optimization is performed separately based on the GSO algorithm and Gamultiobj algorithm. The comparison results are given in Table 1 and Table 2. It can be seen that the parameter identification ability of Gamultiobj algorithm is better than that of GSO algorithm.
In order to compare the optimization result of the two optimization algorithms more intuitively, for the two hysteresis curve models in Figure 1, compare the convergence of the RMSE values based on the Gamultiobj algorithm and the GSO algorithm as shown in Figure 5. The computational efficiency of the Gamultiobj algorithm and the ability to find the global minimum of the objective function are stronger than the GSO algorithm. In addition, Figure 6 shows the change process of each parameter with the iterative in the optimization process of the hysteretic loop 1 using the two algorithms.
In order to compare the optimization capabilities of the two optimization algorithms more accurately, the two optimization algorithms were run 10 times for the hysteresis loop 1,  and the RMSE value was recorded, and the box diagram are shown in Figure 7. The Gamultiobj algorithm has a better optimization effect on the model in this paper.
When using the Gamultiobj algorithm for optimization, the algorithm crossover rate is set to the default value of 0.8. In order to explore whether different crossover rates can improve the algorithm optimization ability, the algorithm crossover rates are set to 0.6, 0.7, and 0.9 to run 10 times respectively, and the results obtained are compared. Figure 8 shows different crossover rates and the convergence of RMSE values under different crossover rates when the results of 10 runs are averaged. The optimization effect is best when the crossover rate is 0.8, so it is reasonable to take the crossover rate 0.8.

E. NONLINEAR DYNAMIC MODEL OF CLAMP
At present, there are two types of clamps widely used in the external pipeline of aero-engines, steel clamps and gasketed clamps. The gasketed clamp, as shown in Figure 9, can adjust its own rigidity and has a shock absorption function. The material of the gasket is a nonlinear material with large damping. The friction between the gasket and the pipeline makes the clamp nonlinear.
The clamp-mass test system is shown in Figure 10(a). Since the stiffness of the mass is much greater than the stiffness of the clamp, it can be considered that the mass is absolutely rigid, which is simplified to the concentrated mass m. The supporting effect of the clamp on the mass is transformed into a damper and a spring support. Due to the characteristics of the metal rubber material, the elastic support of the clamp to the mass is not completely linear. When the clamp deformation reaches the yield limit, its stiffness will no longer be constant, but will show nonlinearity. In order to describe this characteristic, the Bouc-Wen model VOLUME 9, 2021   in Eq. (4) is introduced. The support of the clamp to the mass is decomposed into a damper with a damping of c, a spring support with a stiffness of k f , and the nonlinear stiffness k w .
The simplified model is shown in Figure 10(b), the corresponding governing equation is The form of Eq. (9) is completely consistent with Eq. (4). That is to say, given the mass of the mass m, acceleration excitationü g , and a hysteresis curve describing the nonlinear characteristics of the clamp as shown in Figure 1, a set of parameters = [ρ, σ, n, k f , k w , c] can be obtained by the optimization method described above. The group parameters include the inherent characteristics of the clamp, the stiffness and damping after yielding. Furthermore, by substituting this set of parameters into the dynamic equation of a certain  vibration system, the nonlinear constraint of the clamp can be added to the system. Figure 11(a) shows the horizontal static stiffness test rig of the double metal clamps. Turn the handle (driving device), the test device converts the rotary motion into linear motion, so that the clamp is stretched in the direction of the force. The force sensor shows the loading force of the clamp, and the displacement sensor shows the deformation of the clamp. By applying periodic loading forces of 0N, 60N, 100N, 60N, 0N, −60N, −100N, −60N, 0N to the clamp and recording the corresponding deformation, the hysteresis curve can be obtained as shown in Figure 11(b).

III. PARAMETRIC MODELING CLAMPS A. STATIC STIFFNESS TEST OF CLAMP
The test rig can test the static stiffness and damping of different types and specifications of clamps. Figure 12 shows the clamping methods of single clamp and double clamp respectively. Changing the joint rod of the fixed clamp can adapt to clamps with different inner diameters.

B. PARAMETER IDENTIFICATION OF CLAMP NONLINEAR MODEL BASED ON GAMULTIOBJ ALGORITHM
For Eq. (9), set m = 1kg,ü g = 10.2g sin(20πt), and use the Fourth order Runge-Kutta method to solve, set the time length t = 2s, and the step length dt = 0.001s. Since the excitation period T = 0.1s, a total of 20 periods can be calculated. The total data length is 2000, and one period contains 100 data points.
Eq. (6) is the objective function of the optimization algorithm, F is the restoring force, F 0 is the loading force, x is the displacement, x 0 is the clamp deformation, and L is the data length. Obviously, in the static stiffness experiment, there are only 9 data points in one period of the loading force,  which cannot meet the requirements of the objective function. Since the excitation is a sine function, the displacement and restoring force obtained in Section II(D) also meet the sine law in one period, so the cubic spline interpolation is used to increase the number of loading force and clamp deformation data points. Since one period contains starting point and end point should be 101 data points, so it is necessary to increase the load force and deformation data points to 101.
Five single clamps with an inner diameter of 8mm in Figure 12(a) were subjected to multiple static stiffness experiments, and 12 groups of experimental data were selected and shown in Table 3 and Table 4. It should be pointed out that the initial point of the experimental curve in Figure 9(b) is (0, 0) point, because no force is applied at this time, and naturally there is no deformation. However, such a hysteresis curve is not closed. In order to better match the model, the initial point of the experimental curve is modified to be consistent with the end point.
When using the Gamultiobj algorithm for parameter identification, the parameter range needs to be given first. The selection of the value range of each parameter is discussed below.
It can be seen from Fig. 13 that F el = α(F y /u y )u = αk i u = k f u, which represents the elastic force part of the restoring force; F y is the yield stress; u y is the yield displacement;   u represents displacement; k f is the post yield stiffness; k i is the pre yield stiffness; α is the ratio of post yield stiffness to pre yield stiffness. F h = (1 − α)k i z, which represents the hysteretic force part of the restoring force, and the restoring force R = F el + F h . It can be seen from the figure that k f is the post yield stiffness, and its range can be estimated from the loading force displacement experimental curve.
Then determine the value range of parameter c. c is damping, and its value range can be calculated by formula according to the area surrounded by the experimental curve. The following figure shows a typical restoring force displacement hysteresis curve. It is divided into upper and lower parts to calculate the energy consumed by loading and unloading one  cycle damping w and the stored maximum elastic energy w.
The solution formula of specific damping SDC is The relationship between specific damping and damping ratio is The calculation formula of equivalent viscous damping c eq is After the equivalent viscous damping is obtained, the damping range can be given. When the area of hysteretic curve is little different, the damping value range can remain unchanged.
It can be seen from Eq. (5) that n ≥ 1, and the larger the value of n is, the closer the hysteretic curve model is to the double broken line model. In this paper, the value range of n is [1,3]. What is difficult to determine is the value range of parameter k w , ρ and σ . These three parameters are strongly related to each other, and their values will affect each other.
In the example part of reference [24], it is mentioned that ρ is the reciprocal of yield displacement, which is wrong. Because the definition of ρ in this document is consistent with that in this paper, it is ρ = A/z 0 , should be dimensionless. In reference [25], however, the definition of ρ is ρ = A/z 0 /u y , its dimension is 1/displacement, u y is yield displacement. Therefore, there should be little difference between A and z 0 . In this paper, the range of parameter ρ is in the order of 10 0 or 10 ±1 .
According to the value taking experience of other relevant Bouc-Wen model documents, the order of magnitude of the value of A is generally 10 1 . According to the expression of parameter k w in Eq. (4), the order of magnitude of k w value is not different from k f . After determining the value range of these parameters, the value range of parameter σ is determined according to the curve fitting results of the simulation program.
It should be noted that the definitions of parameters in this paper don't conflict with those in reference [25]. Although their definitions of parameter ρ are inconsistent, after formula transformation, the overall vibration equation Eq. (1) is consistent.
After determining the order of magnitude of the parameters, the parameter range is adjusted according to the parameter identification results and fitting effect of the optimization algorithm. For parameters with clear physical meaning, such as k f and c, make them fit the actual situation. At this time, the values of other parameters are the identification results.
The flow chart of the parameter identification program is shown in Figure 15. The evaluation requirements for the operation of this program are as follows. From the fitting effect, the fitting curve should match the experimental curve. From the parameter identification results, it is necessary to keep the identified parameter results as far away from the value boundary as possible.
In addition, the Pareto solution set obtained by the Gamultiobj algorithm is a non-inferior solution of multiple objective functions, and the optimal solution should be selected according to one's own preferences and actual needs. For the model in this paper, the parameter identification results should conform to reality. For example, the smaller the displacement, the greater the stiffness; The larger the area of hysteresis curve, the greater the damping.

C. PARAMETER IDENTIFICATION RESULT
The parameters of the experimental data in Table 3 and Table 4 are identified according to the parameter identification program, and the identification results are shown in Table 5 and Table 6. At the same time, the matching effect of the fitted curve and the original curve is drawn, as shown in Figure 16.

1) PARAMETER DISTRIBUTION
Since the experimental data in Table 3 and Table 4 is obtained by experimental testing of 5 clamps with an inner diameter of 8mm, the parameters obtained in Table 5 and Table 6 should satisfy a certain probability distribution. Figure 17 analyzes the probability distribution that each parameter may satisfy. The '×' indicates the location of the parameter value. The points on the red dashed line are the fitting points to the parameter data, which meet the normal distribution, lognormal distribution, Weibull distribution and Gumbel distribution, respectively. The title has marked the two characteristic parameters of each distribution.
From the results, the parameters ρ, σ , and n tend to satisfy the normal distribution or log-normal distribution, and the parameters k f , k w , and c tend to satisfy the Gumbel distribution.

IV. CONCLUSION
In this paper, the nonlinear hysteretic restoring force expression of a clamp vibration system is determined based on the Bouc-Wen model, and the restoring force-displacement hysteretic curve model is established. For this model, a numerical example is used to compare the GSO optimization algorithm and the Gamultiobj multi-objective genetic algorithm parameter identification ability. It is found that the Gamultiobj algorithm not only makes the fitted curve fit the original curve better, but also has reliable parameter identification capabilities. In the numerical example, the maximum error of Gamultiobj algorithm parameter identification is 2.8%.
Based on the static stiffness experiment of the clamp, a hysteresis curve describing the nonlinear stiffness of the clamp is obtained. The proposed restoring force-displacement hysteresis model is used to fit experimental curves through the Gamultiobj algorithm, which verifies the effectiveness of the model.
The Bouc-Wen model parameters are obtained, and the identification of the stiffness and damping of the clamp are completed. From the identification results, the stiffness parameters obtained by the recognition of the smaller deformation of the clamp data are larger. For clamp data with a larger hysteresis curve area, the damping parameters obtained by identification are larger, which is consistent with the actual situation.
JUNZHE LIN received the Ph.D. degree in mechanical engineering from Northeastern University, Shenyang, China, in 2019. He is currently a Lecturer with Northeastern University. He is mainly engaged in research work in the fields of nonlinear dynamic and hydrodynamics.
ZHIHONG NIU is currently pursuing the master's degree in mechanical engineering with Northeastern University, Shenyang, China. His current research interest includes probabilistic modeling of pipe systems.
XUFANG ZHANG received the Ph.D. degree in civil engineering from the University of Waterloo, Waterloo, ON, Canada, in 2013. He is currently a Professor with Northeastern University. He is mainly engaged in research work on structural reliability methods and big data modeling.
HUI MA received the Ph.D. degree in mechanical engineering from Northeastern University, Shenyang, China, in 2007. He is currently a Professor with Northeastern University. He is mainly engaged in research work in the fields of rotor dynamic, nonlinear dynamic, and fault diagnosis.
YULAI ZHAO is currently pursuing the Ph.D. degree with Northeastern University, Shenyang, China. He is the author of eight articles. His current research interests include dynamic modeling, nonlinear dynamics, and fault diagnosis.
QINGKAI HAN received the Ph.D. degree in mechanical engineering from Northeastern University, Shenyang, China, in 1997. He is currently a Professor with Northeastern University. He is mainly engaged in research work in the fields of rotor dynamic, nonlinear dynamic, and multibody dynamic. VOLUME 9, 2021