Time-Variation Damping Dynamic Modeling and Updating for Cantilever Beams with Double Clearance Based on Experimental Identification

Abstract

The accuracy of a space manipulator’s end trajectory and stability is significantly affected by joint clearance. Aiming to improve the prediction accuracy of vibration caused by clearance, a dynamic clearance modeling method is developed based on parameter identification in this study. First, a dynamic model framework for manipulator arms is established based on the Hamilton principle and hypothetical mode method with time-variation damping. Then, a multi-resolution identification is performed for identifying the instantaneous frequency and damping ratio to estimate stiffness and damping by the sensors. The quantum genetic algorithm (QGA) is used to optimize the scale factor, which determines the identification accuracy and calculation efficiency. Finally, a case study is conducted to verify the presented model. In comparison with the initial dynamic model based on constant damping, the modal assurance criterion (MAC) of the proposed improved model based on time-variation damping is improved by 43.97%, the mean relative error (MRE) of the frequency response function (FRF) is reduced by 32.6%, and the root mean square error (RMSE) is reduced by 18.19%. The comparison results indicate the advantages of the proposed model. This modeling method could be used for vibration prediction in control systems for space manipulators to improve control accuracy.

1. Introduction

In order to improve the performance and service life of spacecraft, in-orbit service technologies such as fuel filling, maintenance, repair, and function module replacement are urgently needed. In complex space environments, the motion accuracy and stability of a space manipulator directly affect the in-orbit service technology performance. Due to wear and moving pairs, clearance exists widely in the structure of a space robot arm, causing nonlinear vibration of the robot arm, reducing the mechanism instability, positioning accuracy, and trajectory prediction [1,2,3].

As the main factor affecting nonlinear vibration, clearance dynamics have been studied extensively by scholars. Machado compared the application and accuracy of the Hunt–Crossley, Lee–Wang, Lankarani–Nikravesh, and Herbert–McWhannell contact force models [4]. Based on the Lankarani–Nikravesh and Winkler models, Bai and Zhao developed a mixed contact force model for studying the dynamic behavior of planar mechanical systems with joint clearance [5]. Ambrósio and Flores simulated the contact between the shaft neck and bearing based on the Lankaran–Nikravesh contact force model [6]. Chen established an improved contact force model considering clearances and studied the influence of clearance size and cylinder support stiffness on dynamic response [7]. Chen utilized the continuous force model proposed by Lankarani and Nikravesh to characterize the contact–collision phenomenon with the contact force model and examined the dynamic response under various driving speeds and clearance sizes [8]. Chen utilized the Lagrange multiplier method to formulate a dynamic model of a multi-link complex mechanism with multiple lubrication gaps. A comparative and analytical study was conducted on the impact of dry friction clearance and lubrication clearance on the dynamic characteristics of the mechanism [9]. Marques proposed a new formula to analyze the relationship between the spatial rotary joint and the radial and axial clearance. The test results show that the mechanical joint clearance is of great significance to the dynamic response of the system [10].

Considering the effect of contact collision and clearance, Xie developed a rigid–flexible coupling dynamic model with clearance to study the factors affecting angular velocity, which can provide a dynamic response with different clearances. It should be noted that angular velocity is affected by contact collisions, and the relationship between clearance and response is nonlinear [11]. To develop an improved clearance model for predicting end vibration, Liu established a cantilever beam dynamic model with double clearance. The relationship between the clearance and disturbance of the end was described quantitatively using dynamic information to identify the clearance parameters [12]. WANG established an improved non-lubricated nonlinear contact force model, an improved lubricated nonlinear transition force model, and studied the impact dynamics characteristics and parameters of the hinge mechanism with clearance [13].

In practice, clearances have often been filled with grease to ensure relative sliding. Therefore, to simplify the analysis, the influence of nonlinear lubrication factors has often been ignored. However, as the grease shape changes in the contact position, the damping obtains a nonlinear form over time. Therefore, considering nonlinear factors such as clearance lubrication, the time-varying damping parameter is introduced to describe varying characteristics of the clearance force more accurately. In this way, the shortcoming of the constant damping parameter is eliminated. In recent years, time-varying parameters have been widely used to describe nonlinear characteristics in many fields.

To describe the change in damping, Sun has developed a magnetorheological damper with adjustable stiffness and damping for automotive suspensions. Variable stiffness and damping suspensions can reduce the sprung mass acceleration effectively [14]. As a result, vehicle performance and comfort can be improved. Kim proposed an energy-shaping speed controller by injecting time-varying damping into the speed loop. In the proposed controller, offset errors were eliminated without using an integral tracking error, and a better suppression of interference was achieved. According to the experimental results, time-varying damping could enhance performance by 15% [15]. Matteo Scapolan introduced a time-varying damper to achieve parameter resonance in an electromechanical oscillator for energy collection [16]. To improve the bandwidth and performance of a vibration energy harvester, time-varying damping was used as a control strategy. Based on optimal control and a pseudo-spectral decomposition, Giuseppe Giorgi developed an algebraic formula to calculate the steady-state response by introducing time-varying damping. In addition, time-varying damping was described well by analyzing the problem of nonlinear factors [17]. Li focused on the effect of damping changes on the vibration isolation of a quasi-zero-stiffness vibration isolator. The time-varying damping can be equivalent to the addition of a stiffness term to the vibration system. The system declined in its vibration stability [18].

Due to the presence of clearances, there is mutual contact between the various components and joints of the mechanism. Friction, wear, and deformation of components at the hinge during use further increase the clearance size, leading to discrepancies with theoretical calculations, thus adding nonlinear characteristics to the system. Additionally, to ensure the stability of motion, substances such as lubricants are added within the clearances, which further exacerbate the nonlinear characteristics of the clearance position, resulting in a decrease in the control accuracy of the system.

Therefore, to improve control precision, it is essential to fully consider the influencing factors in the dynamic model. Reducing errors in the dynamic model is a prerequisite for establishing an accurate dynamic model.

Due to the existence of nonlinear factors, such as clearance and lubrication, and the presence of constant parameters in the current clearance impact model, it is difficult to accurately describe the dynamic characteristics of the real model through the established dynamic model with constant parameters. This paper proposes a three-parameter time-varying damping model by combining the Rayleigh and complex damping models to describe the dynamic characteristics for clearance. The dynamic model is decoupled, and modal information is retained. In addition, time-varying characteristics can accurately describe nonlinear vibrations caused by clearance and lubrication. Consequently, an error between the dynamic and real models can be avoided.

The paper is organized as follows. In Section 1, a dynamic model framework for a cantilever beam with clearance is developed based on the regularization theory using the modal analysis method. Also, a dynamic clearance model with time-varying stiffness and damping is constructed. In Section 2, based on wavelet transform, the time-varying parameters in the clearance model are identified, and the scale factors are optimized and verified to improve the identification accuracy and computational efficiency. In Section 3, a test model is constructed to verify the correctness of the modeling method through error comparison of the frequency response function (FRF), cross signature assurance criteria (CSAC), and cross signature scale factor (CSF). In addition, the effectiveness is confirmed by the model verification and validation theory. Finally, the conclusions are drawn in Section 4.

2. Materials and Methods

The space manipulator is abstracted as a cantilever beam model with double-clearance (CB-DC). In the existing studies, damping’s influence on dynamic characteristics has typically been ignored. However, a simplified model cannot describe the vibration of the CB-DC model well [19]. As shown in Figure 1, a dynamic model is constructed to reflect the vibration characteristics accurately. As a continuum constituted of particles along the x-axis, the CB-DC model is fixed at the left end, whereas the right end can move freely. The clearance is regarded as a spring damper. The Parameter List and Abbreviations are in Appendix E,

Table A13. Formula symbol.

Symbol	Parameter
$l_1$	Distance from a free beam’s fixed end $O$ to clearance 1
$l_2$	Distance from a free beam’s fixed end $O$ to clearance 2
$l_f$	Distance from the fixed end of a free beam to the excitation point
$l$	Length of a free beam
$E$	Young’s modulus of a free beam
$A$	Cross-sectional area of a free beam
$I$	Cross-section moment of inertia
$\rho$	Density of a free beam
$d_b$	Outer diameter of a free beam
$d_l$	Inner diameter of a free beam
$k_{d1}$	Equivalent contact stiffness at clearance 1
$k_{d2}$	Equivalent contact stiffness at clearance 2
$c_{d1}$	Equivalent contact damping at clearance 1
$c_{d2}$	Equivalent contact damping at clearance 2
$u\left(x,t\right)$	Vertical displacement of a particle of a free beam
$Y_i\left(x\right)$	Modal shape function
$p_i\left(t\right)$	Modal coordinates
$M_{pi}$	Modal mass of the ith free beam
$K_{pi}$	Modal stiffness of the ith free beam
$C_{pi}$	Modal damping of the ith free beam
$f_1\left(t\right)$	Nonlinear force of clearance 1
$f_2\left(t\right)$	Nonlinear force of clearance 2
${\eta }_1$	Rayleigh damping parameter 1
${\eta }_2$	Rayleigh damping parameter 2
$g_1\left(t\right)$	Displacement in clearance 1
$g_2\left(t\right)$	Displacement in clearance 2
${\dot{g}}_1\left(t\right)$	Velocity in clearance 1
${\dot{g}}_2\left(t\right)$	Velocity in clearance 2
${\alpha }_0$	Rayleigh damping model coefficient in clearance 1
${\alpha }_1$	Rayleigh damping model coefficient in clearance 2
${\alpha }_3$	Complex damping model coefficient of clearance
${\xi }_{d1}\left(t\right)$	Instantaneous damping ratio in clearance 1
${\xi }_{d2}\left(t\right)$	Instantaneous damping ratio in clearance 2
$H\left[x\left(t\right)\right]$	Hilbert transformation of a signal
$A\left(t\right)$	Signal amplitude
$\phi \left(t\right)$	Signal phase
$W\left(a,b\right)$	Wavelet coefficient
${\psi }_{a,b}\left(t\right)$	Morlet wavelet basis function
${\omega }_0$	Center frequency of a wavelet
$a$	Scale factor of the wavelet transform
$b$	Translation factor of the wavelet transform
${\epsilon }_a^{nor}$	Deviation degree between the identification and theoretical results
$T_a^{nor}$	Parameter identification time
$y_{best}^t$	Optimal fitness value under the current iteration number
$H_{Xi}\left(\omega \right)$	The ith FRF in the test of $\omega$
$H_{Ai}\left(\omega \right)$	The ith FRF in the simulation of $\omega$

and

Table A14. Abbreviations.

Abbreviation	Explanation
CB-DC	Cantilever beams with double clearance
QGA	Quantum genetic algorithm
FRF	Frequency response function
FRAC	Frequency response assurance criterion
CSAC	Cross signature assurance criterion
CSF	Cross signature scale factor
RMSE	Root mean square error
MRE	Mean relative error
MAC	Modal assurance criterion
ARA	Acceptable range accuracy

.

In Figure 1, $l_1$, $l_2$, and $l_f$ indicate the distances from the $O$ to clearance 1, 2, and the excitation point, respectively; $E$, $A$, $I$, $\rho$, and $l$ denote the Young’s modulus, cross-sectional area, cross-sectional moment of inertia, density, and length of the cantilever beam, respectively. The clearance values are the same, $d_1=d_2=\left(d_l-d_b\right)/2$, where $d_b$ and $d_l$ represent the beam diameter and the supporting inner diameter, respectively. The contact stiffness and damping in clearance 1 and clearance 2 are denoted by $k_{d1}$, $k_{d2}$ and $c_{d1}$, $c_{d2}$, respectively. The vertical displacement of a particle can be expressed as follows:

$$u\left(x,t\right)={\displaystyle \sum _1^n Y_i}\left(x\right)p_i\left(t\right),$$

(1)

where $Y_i\left(x\right)$ represents the modal shape function, $p_i\left(t\right)$ indicates the modal coordinates, and $i$ shows the vibration shape order.

The boundary condition is defined by

$${\left.{\displaystyle \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}}\right|}_{x=0}=0; {\left.{\displaystyle \frac{{\partial }^{3}u\left(x,t\right)}{{\partial }^{3}x}}\right|}_{x=0}=0; {\left.{\displaystyle \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}}\right|}_{x=l}=0; {\left.{\displaystyle \frac{{\partial }^{3}u\left(x,t\right)}{\partial x^3}}\right|}_{x=l}=0.$$

(2)

The dynamic equation of the CB-DC can be expressed as follows:

$$M_{pi}{\ddot{p}}_i+C_{pi}{\dot{p}}_i+K_{pi}{p}_i+Y_i\left(l_1\right)f_1\left(t\right)+Y_i\left(l_2\right)f_2\left(t\right)=Y_i\left(l_f\right)f\left(t\right),$$

(3)

where $M_{pi}$, $K_{pi}$, and $C_{pi}$ represent the ith modal mass, modal stiffness, and modal damping of the beam, and they are calculated by Equations (4)–(6), respectively; $f_1\left(t\right)$, $f_2\left(t\right)$, and $f\left(t\right)$ are the nonlinear forces in $l_1$, $l_2$, and $l_f$, respectively.

$$M_{pi}={\displaystyle {\int }_0^l \rho A{Y_i}^2}dx$$

(4)

$$K_{pi}={\displaystyle {\int }_0^l EI{\left({Y_i}^{″}\right)}^2}dx={{\displaystyle {\int }_0^l Y_{i}EI\left({Y_i}^{\prime }\right)}}^{″}dx$$

(5)

$$C_{pi}={\eta }_1{M}_{pi}+{\eta }_2{K}_{pi}$$

(6)

In Equations (4)–(6), the shape function $Y_i\left(x\right)$, which meets the boundary conditions, is specified as

$$Y_i\left(x\right)=\cos {\beta }_{i}x-ch{\beta }_{i}x+{\displaystyle \frac{\sin {\beta }_i-sh{\beta }_{i}l}{\cos {\beta }_{i}l-ch{\beta }_{i}l}} \left(\sin {\beta }_{i}x-sh{\beta }_{i}x\right),$$

(7)

where ${\beta }_i$ represents the ith solution of the frequency equation. The natural frequency can be defined by

$$\cos \beta lch\beta l=-1,$$

(8)

and ${\eta }_1$, ${\eta }_2$ represent the Rayleigh damping parameters.

Considering the damping loss of collision, the Kelvin–Voigt model describes normal collision characteristics [20], which include linear spring-damping elements:

$$f_1\left(t\right)=c_{d1}\left(t\right){\dot{g}}_1\left(t\right)+k_{d1}{g}_1\left(t\right),$$

(9)

$$f_2\left(t\right)=c_{d2}\left(t\right){\dot{g}}_2\left(t\right)+k_{d2}{g}_2\left(t\right),$$

(10)

where $g_1\left(t\right)$, $g_2\left(t\right)$ and ${\dot{g}}_1\left(t\right)$, ${\dot{g}}_2\left(t\right)$ denote displacements and velocity functions at clearances 1 and 2, respectively. Their specific expressions are as follows:

$$g_1\left(t\right)=\left\{\begin{array}{ll}u\left(l_1,t\right)-d_1\hfill & u\left(l_1,t\right) \ge d_1\hfill \\ 0\hfill & -d_1<u\left(l_1,t\right)<d_1\hfill \\ u\left(l_1,t\right)+d_1\hfill & u\left(l_1,t\right) \le -d_1\hfill \end{array}\right.,$$

(11)

$$g_2\left(t\right)=\left\{\begin{array}{ll}u\left(l_2,t\right)-d_2& u\left(l_2,t\right) \ge d_2\\ 0& -d_2<u\left(l_2,t\right)<d_2\\ u\left(l_2,t\right)+d_2& u\left(l_2,t\right) \le -d_2\end{array}\right.,$$

(12)

The clearance dynamics model established based on Hamilton’s principle and modal decomposition can accurately describe continuous dynamic behavior. At the same time, the oil film force generated by lubrication at the clearance is considered, which keeps the joint in a separated state, suppressing the strong vibrations and impacts caused by contact forces at the clearance. By simplifying the clearance unit into a linearized spring-damping model, it is possible to effectively analyze the influence of clearance factors on the dynamic system under conditions where the vibrations are stable. However, for highly complex and high-precision systems, as well as for cases where large clearances lead to intense collisions, the clearance contact model established in this paper has certain limitations.

The dynamic characteristics of the clearance position exhibit high nonlinearity, making it difficult for the model to accurately describe the clearance dynamics. This leads to the inability of the control system to precisely achieve the robotic arm’s motion trajectory and working position.

In the clearance dynamics modeling, by performing sensitivity analysis on the structural parameters, it was determined that the main dynamic parameter influencing the robotic arm vibration is the damping parameter between the clearances. Therefore, effectively and reasonably determining the clearance damping parameters is crucial for ensuring the correctness of the entire system.

In the clearance vibration analysis, the internal friction of solid materials and decoupling are considered. The time-varying damping model is constructed by other damping models [21]. This model is then decoupled by dynamic equations; it retains modal information. The constructed model can be expressed as follows:

$$c(t)={\alpha }_{0}m+{\alpha }_{1}k+{\alpha }_{3}jk$$

(13)

where $c(t)$ is the general expression of clearance parameters $c_{d1}(t)$ and $c_{d2}(t)$; $m$ and $k$ are the equivalent mass and stiffness; and ${\alpha }_0$, ${\alpha }_1$, and ${\alpha }_3$ denote the Rayleigh damping model coefficients. The parameters are expressed as follows:

$$\begin{array}{l}{\alpha }_0=2\omega \xi \\ {\alpha }_1=2\xi / \omega \\ {\alpha }_3=2\xi \end{array}$$

(14)

where $\xi$ and $\omega$ indicate the damping ratio and frequency, respectively.

To ensure stable movement, the clearance is filled with grease during the collision and recovery processes. In the process of collision, a change in the grease shape results in a complex interaction at clearance. In the Kelvin–Voigt model, the damping coefficient is constant, so this model cannot accurately model actual regular.

Lubrication in the clearance can alleviate the impact and vibration between moving parts during the motion, suppress high-frequency oscillations, and provide a good buffering effect. When oil film lubrication is considered, the dynamic characteristics of the mechanism are significantly improved compared to the non-lubricated condition. The high-frequency oscillation phenomenon of the contact force at the clearance is significantly suppressed at the extremal positions of the mechanism. The dynamic output characteristics, such as displacement, velocity, and acceleration, closely match the ideal values [13].

A time-varying spring-damping model is proposed to overcome the limitations, considering the time-varying characteristics of the damping ratio when clearances and lubrication are present. Therefore, a nonlinear damping model can be expressed by

$$c_{d1}\left(t\right)=\left(2\omega {\xi }_{d1}\right)m+\left(2{\xi }_{d1}/\omega \right)k+\left(2{\xi }_{d1}\right)jk,$$

(15)

$$c_{d1}\left(t\right)=\left(2\omega {\xi }_{d2}\right)m+\left(2{\xi }_{d1}/\omega \right)k+\left(2{\xi }_{d2}\right)jk,$$

(16)

where ${\xi }_{d1}$ and ${\xi }_{d2}$ represent the instantaneous damping ratios at clearances 1 and 2, respectively, and $\omega$ is the instantaneous frequency.

In this section, the dynamic model of the CB-DC is established based on the Hamilton principle. Using time-varying parameters, the clearance model is constructed to describe the effect of nonlinear factors on vibrations. In the time-varying damping model, the instantaneous damping ratio is crucial for achieving high accuracy. The damping ratio represents the attenuation form of the vibration after excitation. It is difficult to determine the accurate value of the instantaneous damping ratio through theoretical modeling, so a test method has often been used. To create a more accurate model, the instantaneous damping ratio will be determined using the parameter identification method in the next section.

3. Time-Varying Damping Identification

3.1. Identification Based on Wavelet Transform

In view of the typical time-varying system response of the CB-DC structure model, the damping parameters are time-varying using time–frequency analysis for identification in the time-varying system. Therefore, an identification method combining modal analysis and wavelet transform is proposed in this chapter. Wavelet transform can be thought of as a window that can be scaled to different sizes. It helps us observe signals at various scales, providing detailed information about the signal in both time and frequency domains. The vibration signals of clearance contact are non-stationary in nature, and wavelet transform is particularly effective in identifying signals of different frequencies. The identification process is shown in Figure 2.

The first part is signal denoising processing. The combination of the EEMD and wavelet threshold denoising achieves signal denoising. The second part is the identification of system parameters. By combining the analytical mode and Hilbert transform, the frequency aliasing and intermittent signal fluctuation are realized.

In view of the time-varying characteristics of clearance and lubrication, the collision force of clearance is nonlinear, and its details cannot be described well by the identified curve of vibration signals. In the parameter identification method, the wavelet analysis can provide an adaptive time–frequency window. Adaptive windows can automatically narrow and widen according to the high- and low-frequency signals, respectively. Thus, wavelet analysis can provide highly localized information in signal processing [22]. In this study, the wavelet ridge and wavelet skeleton are used to derive the instantaneous damping ratio and frequency to construct a time-varying damping model. Since the output characteristic of the constructed time-varying damping model is non-stationary, the time–frequency analysis of the response signal can be an effective method.

Assume that under the CB-DC excitation, the test signal at a point is denoted by $g\left(t\right)$, and wavelet transform is denoted by $W\left(a,b\right)$. A series of points $\left(a_i,b_i\right)$ form the ridges of the wavelet transform. At these points, the phase is stable, and the wavelet coefficients are maximum. The coefficient envelope on the wavelet ridge is defined as a skeleton. Changes in the frequency and damping ratio can be estimated by the wavelet ridge and skeleton [23].

First, the relationship between the wavelet ridge, signal instantaneous frequency, and amplitude is established. Assume the test signal is defined by

$$g\left(t\right)=A\left(t\right)e^{j\phi \left(t\right)}=x\left(t\right)+jH\left[x\left(t\right)\right]$$

(17)

where $H\left[x\left(t\right)\right]$ represents the Hilbert transform of $x\left(t\right)$, having an amplitude of $A\left(t\right)=\sqrt{x^2\left(t\right)+H^2\left[x\left(t\right)\right]}$ and a phase of $\phi \left(t\right)=\arctan \left({\displaystyle \frac{H\left[x\left(t\right)\right]}{x\left(t\right)}}\right)$.

Next, the relationship between the wavelet coefficients, instantaneous frequency, and amplitude is established. The analytic signal is transformed by the wavelet transform to solve the wavelet coefficient $W\left(a,b\right)$ as follows:

$$W\left(a,b\right)=〈x\left(t\right),{\psi }_{a,b}\left(t\right)〉={\displaystyle \frac{1}{2\sqrt{a}}}{\displaystyle \underset{+\infty }{\overset{-\infty }{\int }} A\left(t\right)e^{j\phi \left(t\right)}{\psi }^{\ast }\left(\frac{t-b}{a}\right)}dt,$$

(18)

where ${\psi }_{a,b}\left(t\right)$ is the Morlet wavelet basis function.

The Morlet wavelet basis function is adopted because its shape is similar to that of the vibration characteristics. By substituting the Morlet wavelet basis function $\psi \left(t\right)=e^{-{\displaystyle \frac{t^2}{2}}}{e}^{j{\omega }_{0}t}$ and Fourier transform of the Morlet wavelet $\mathrm{\Phi }\left(\omega \right)=\sqrt{2\pi }{e}^{-{\displaystyle \frac{{\left(\omega -{\omega }_0\right)}^2}{2}}}$ into Equation (18), the wavelet coefficients can be obtained as follows:

$$W\left(a,b\right)={\displaystyle \frac{\sqrt{a}}{2}}A\left(b\right)e^{-{\displaystyle \frac{N}{4}}{\left(a\dot{\phi }\left(b\right)-{\omega }_0\right)}^2}{e}^{j\phi \left(t\right)},$$

(19)

Finally, the relationship between the damping ratio and instantaneous frequency and amplitude is established. The response signal can be expressed by

$$x\left(t\right)=e^{-\xi \omega t}\cos \left(\omega t+\phi \right),$$

(20)

and the corresponding instantaneous damping ratio based on the response equation is given by

$$\xi =-{\displaystyle \frac{d\left[In\left(A\left(t\right)\right)\right]}{\omega dt}};$$

(21)

meanwhile, the corresponding instantaneous amplitude and frequency are, respectively, given by

$$\begin{gathered} \omega ={\displaystyle \frac{{\omega }_0}{a}}, \\ A\left(t\right)={\displaystyle \frac{2\left[W\left(a,b\right)\right]}{\sqrt{a}}}. \end{gathered}$$

(22)

In Equation (24), ${\omega }_0$ is the center frequency of the wavelet, $a$ and $b$ are the scale and translation factors, respectively, and $W\left(a,b\right)$ represents the wavelet coefficient. The function of the scale factor is to either stretch or compress a wavelet under defined similarity conditions. The time-varying damping model can be defined by

$$c_{d1}\left(t\right)={\displaystyle \frac{1}{\pi }}{\left(In(A_{d1}(t))\right)}^{\prime }m+\left({\displaystyle \frac{a^2}{{{\omega }_{d1}}^2}}{\left(In(A_{d1}(t))\right)}^{\prime }\right)k+j{\displaystyle \frac{1}{\pi }}\left({\displaystyle \frac{a}{{\omega }_{d1}}}{\left(In(A_{d1}(t))\right)}^{\prime }\right)k$$

(23)

$$c_{d2}\left(t\right)={\displaystyle \frac{1}{\pi }}{\left(In(A_{d2}(t))\right)}^{\prime }m+\left({\displaystyle \frac{a^2}{{{\omega }_{d2}}^2}}{\left(In(A_{d2}(t))\right)}^{\prime }\right)k+j{\displaystyle \frac{1}{\pi }}\left({\displaystyle \frac{a}{{\omega }_{d2}}}{\left(In(A_{d2}(t))\right)}^{\prime }\right)k$$

(24)

In this section, a time-varying damping model is constructed based on the identified instantaneous damping ratio using the wavelet transform. The scale factor represents a scale function coefficient in the identification process of design parameters. The scale factor defines the resolution ratio and calculation complexity of the frequency and time axis to facilitate parameter identification [24]. On the basis of interval selection and discretization, scale parameters are determined using an optimization algorithm to ensure calculation efficiency and accuracy.

3.2. Identification Optimization

Compared to the conventional path optimization algorithms, genetic algorithms are universal and adaptable [25]. The quantum genetic algorithm (QGA) [26] has a wide solution space coverage and high search efficiency for this study. The quantum genetic algorithm (QGA) is an optimization method that combines quantum computing with genetic algorithms. It leverages the parallelism and quantum superposition properties of quantum computing to improve the search efficiency of traditional genetic algorithms. Based on the QGA, the optimal scaling factor parameters can be searched in a shorter amount of time. According to the QGA, probability amplitudes represent chromosome genes, and each qubit may represent either one, zero, or a superposition state of one.

In the process of scale factor optimization, the objective function is weighted by the error and time indices: (1) deviation degree ${\epsilon }_a^{nor}$ between the identification and theoretical result, evaluated by the root mean square error (RMSE); (2) time for parameter identification $T_a^{nor}$. For a more uniform scale, two evaluation indices are normalized to [0, 1]. Therefore, the objective function can directly combine two weighted terms without considering the inconsistency in their dimensions. The optimized design variable is a scale factor of the wavelet transform. The design space $[a^{\min },a^{\max }]$ is determined by empirical value analysis. Mathematically, the optimization problem is defined by

$$\begin{array}{cc}obj.\hfill & \underset{a}{\min } Y_a=e_a^{nor}+T_a^{nor}\hfill \\ where\hfill & e_a^{nor}=\left|{\displaystyle \frac{e_a-e_a^{\min }}{e_a^{\max }-e_a^{\min }}}\right|,e_a=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N{({\tilde{c}}_i-c_i)}^2}},T_a^{nor}=\left|{\displaystyle \frac{T_a-T_a^{\min }}{T_a^{\max }-T_a^{\min }}}\right|\hfill \\ s.t.\hfill & a^{\min } \le a \le a^{\max }\hfill \end{array},$$

(25)

where ${\epsilon }_a$ is the RMSE of the constructed damping $\tilde{c}$ and theoretical damping $c$, $N$ is the sample number, and $T_a$ is the identification time.

As a result, a higher performance can be achieved under the condition of small populations. The optimization algorithm begins the iteration once the following conditions are met:

$$\left\{\begin{array}{l}\left|{\displaystyle \frac{y_{best}^t-y_{best}^{t-10}}{y_{best}^t}}\right|<1\times {10}^{-4}\hfill \\ k>k^{max}\hfill \end{array}\right.,$$

(26)

where $k$ is the current iteration number, $y_{best}^k$ represents the best fitness value under the current iteration number, and $k^{\max }$ denotes the maximum number of iterations; in this study, $k^{\max }$ is set to 100.

The iterative process is shown in Figure 3.

As shown in Figure 3, first, the population with scale $n$ is randomly generated. According to the probability amplitude value of population, $R(k)$ is constructed. Each $r_j$ is a binary string with a length of m. Then, the following strategy is adopted. A random number $\gamma$ from the range of $[0,1]$ is generated for each qubit of each chromosome. If the norm square value $r_i$ of the corresponding bit is larger than $\gamma$, the bit is set to “1”; otherwise, it is set to “0”. Then, the fitness value of each individual in $R(k)$ is calculated. If the stop condition is met, the algorithm flow is terminated. The population maintains the optimization and re-updates $a(k)$ if the population requires a catastrophe; otherwise, the population uses the quantum revolving gate $U(k)$ to update $a(k)$. Then, $a(k)$ is updated by quantum crossover and mutation; $R(k)$ is constructed for $k=k+1$, and $R(k)$ is constructed again. Finally, the optimization iteration continues until the termination conditions are met to finish the scale factor parameter optimization of wavelet transform.

3.3. Parameter Identification Verification

To verify the identification method validity, two types of time-varying damping models with two degrees of freedom were analyzed, which have a mass of 5 kg and a spring stiffness of 800 m/N. One of the time-varying damping models has the step-change damping function, and the other has the continuous-change damping function. By using these functions, the identification accuracy is further verified.

Step–change damping: $c(t)=0.06 \mathrm{Ns}/\mathrm{m}\left(t\le 50\right)$, $c(t)=0.01 \mathrm{Ns}/\mathrm{m}\left(t>50\right)$.

Continuous-change damping: $c(t)=0.06-0.0002t$.

The RMSE and convergence times were used to evaluate the identification efficiency and accuracy of the multiscale parameter identification method [27] and quantum genetic optimization identification method. The accuracy and efficiency of the parameter identification process were selected as an objective function. The optimal scale factor was set to $a=2^{-6.8}$, based on the multiscale parameter identification result, and $a=2^{-6.68}$, based on the quantum genetic optimization result. The identification results are presented in Figure 4a, where it can be seen that the identification result was more consistent with the theoretical value after parameter optimization. After identification, the full-time-domain RMSE was 0.0037, demonstrating a reduction of 16.49%. The convergence time was 2 s, which indicated a decrease of 42.17%. Under the sudden change in the damping coefficient, the identification method could realize fast-tracking. The tracking time was reduced by 26.24%, and the full-time-domain RMSE was 0.0017.

Further, the residual curve is introduced to illustrate the identification error intuitively at each point. The residual curves of the continuous and step damping processes are presented in Figure 4b. During the initial stages of parameter identification, the algorithm converged rapidly. The error increased when the damping ratio suddenly changed. After the algorithm converged to the parameter value after mutation, the overall mean square error was 0.0100. The results indicated that the parameter identification method based on the wavelet transform was highly accurate in identifying the parameters and capable of tracking. Compared with the multiscale parameter identification method, the mean and extreme values of the error curves were significantly reduced. Thus, obtaining an optimal scale factor by the QGA can significantly enhance the performance of the parameter identification method.

Figure 5a shows the stacked histogram of the normalized RMSE over time, and Figure 5b shows the QGA objective function curve with the design variables; blue and red colors indicate accuracy and time, respectively. A quantitative evaluation of time and accuracy is represented by the superposition. According to the multiscale identification method’s results, the interval between the scale factors was large in the stacked histogram. The scale factor curve obtained by the QGA was continuous. The method of QGA could provide better accuracy, efficiency, and scale factor value after weighting.

As an example, 10 experiments were conducted using step damping, and the optimized parameters and identified parameters are shown in

Table 1. Results of multiple experiments.

Test Number	Scale Factor	Damping Standard Value (t < 50 s)	Damping Identification Value	Damping Standard Value (t ≥ 50 s)	Damping Identification Value
1	2 × 106.680	0.06 Ns/m	0.061	0.01 Ns/m	0.012
2	2 × 10−6.679		0.058		0.010
3	2 × 10−6.678		0.062		0.009
4	2 × 10−6.680		0.061		0.011
5	2 × 10−6.680		0.057		0.010
6	2 × 10−6.681		0.059		0.008
7	2 × 10−6.678		0.058		0.009
8	2 × 10−6.680		0.061		0.008
9	2 × 10−6.681		0.060		0.01
10	2 × 10−6.679		0.061		0.01

. The experimental results indicate that, across multiple trials, the scaling factors optimized using the quantum genetic algorithm demonstrate good stability. Additionally, the identified parameters based on these scaling factors are also close to the true values, further proving the stability and reproducibility of the optimization algorithm, as well as the accuracy of the parameter identification.

4. Case Study

In the previous section, the scale factor in parameter identification was optimized by the QGA. By identifying the time-varying damping ratio, a dynamic equation was established to describe dynamic clearance characteristics. The effectiveness of the established model was evaluated by experiments on a test platform. The FRF of the constant damping model, time-varying damping model, and test model were analyzed using the relative error, MAC, CSAC, and CSF as evaluation indices.

4.1. Experimental Environment

The system calculates the time-varying stiffness and damping parameters of the clearance unit based on the dynamic parameters obtained from sensors using wavelet transform methods. A clearance dynamics model based on time-varying parameters is established, which can predict the variation pattern of dynamic parameters in contact collision states. The control system adopts a clearance dynamic compensation method, which compensates for errors based on the variation characteristics of the clearance’s dynamic parameters. This enables real-time adjustment of the joint position by the controller, reducing errors caused by the clearance. The space robotic arm is primarily connected by multiple joints, and the clearance at the base of the robotic arm is an important factor influencing its dynamic characteristics. Using similarity theory, the main arm structure of the space robotic arm is scaled to ensure that the dynamic characteristics of the experimental model are the same as those of the real model, thus demonstrating that the scaling model method is valid and reliable for real-world verification.

The CB-DC model was fixed at its left end and rested on two support sleeves. The effect of different clearances on dynamic characteristics was studied using six groups of contact rings. A = 0.05 mm, B = 0.1 mm, C = 0.5 mm, D = 1 mm, E = 1.5 mm, F = 2 mm, as illustrated in Figure 6.

The constructed test model is presented in Figure 7. In this model, the primary focus is on the clearance dynamics between the robotic arm’s upper arm and the base. The robotic arm and the rocking frame form a sliding pair, constrained by a tightening nut and a spring to limit axial movement. The base and shaft sleeve form a rotational pair, while the base and platform are connected by a high–low rod assembly. A seat ring, containing internal meshing teeth, is installed between the platform and the base, with the directional machine having external meshing teeth. The engagement of these teeth adjusts the direction of the robotic arm. The model particularly focuses on the rotational pair clearance formed by the rocking frame and bushing, as well as the clearances between the base components.

The excitation and control equipment used in the experiment is the ET-40 electric vibration table system. The system consists of a controller, power amplifier, vibration table body, and fan. The data acquisition system uses Bruel Kjaer LAN-XI 3053, with gyroscope and accelerometer models VG91, PCB 352C66, and PCB 356A02, respectively.

The experimental data were obtained by applying two different excitations within the same clearance by adjusting the clearance. One part of the identification data was used to obtain time-varying damping parameters, while the other served as experimental comparison data.

The joint optimization simulation based on Isight software is implemented using the established dynamic model, with the dynamic simulation module and parameter calculation module working in tandem. The clearance unit is modeled using a spring-damping model, and motion parameters are obtained through the dynamic simulation module. These parameters are then passed in real time to the parameter calculation module, which computes the stiffness and damping parameters using a parameter identification method. The computed parameters are then sent back to the dynamic simulation software as input parameters for real-time simulation. This process involves multiple parameter transfers, enabling the real-time acquisition of dynamic characteristics.

4.2. Experimental Result

Using different clearance values of 0.05 mm, 0.1 mm, 0.5 mm, 1 mm, 1.5 mm, and 2 mm, the FRF comparison between the test and simulated model was performed, as shown in Figure 8. Figure 8a–d show the FRF results for clearance 1, clearance 2, the front, and the end at the clearance value of 0.5 mm, respectively; the blue, black, and red represent the results of the initial model with constant damping, the updated model with time-varying damping, and the test model, respectively. For the first and second orders, the FRF results obtained under constant damping agreed well with the results of the test model. However, there were certain errors for the third and fourth orders. These errors were mainly due to nonlinear factors such as clearance lubrication that affected the precision of dynamic parameters. Considering nonlinear factors, the FRF results obtained under time-varying damping had a smaller error in the first three orders, but there were certain errors for the fourth order because of the parameter identification errors. Overall, the error was within acceptable limits.

To further compare the differences in natural frequency, the model correlation analysis was performed by the relative error (RE) and modal assurance criterion (MAC), as shown in

Table 2. Natural frequency errors.

Order	Test Frequency (Hz)	Initial Model			Updated Model
		Frequency (Hz)	RE (%)	MAC	Frequency (Hz)	RE (%)	MAC
1	9.25	9.15	1.0811	0.7338	9.20	0.5405	0.9183
2	39.75	38.40	3.3962	0.8352	40.25	1.2579	0.9049
3	91.15	98.10	7.6248	0.7391	92.35	1.3165	0.8506
4	189.30	192.95	1.9282	0.6507	187.45	0.9773	0.8785

The error of the initial model with constant damping was less than 7.7%, while the error of the updated model based on time-varying damping was less than 1.4%; thus, the maximum error was reduced by 6.3%. The MAC of the updated model was over 85%, which was higher than that of the initial model. The mean MAC of the first four orders of the updated model increased by 14.84% compared to the initial model. The error comparison of other clearance values is provided as shown in Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6 and

Table A1. Natural frequency errors (clearance = 0.05 mm).

Order	Test Frequency (Hz)	Updated Model			Initial Model
		Frequency (Hz)	Error (%)	MAC	Frequency (Hz)	Frequency (Hz)	Error (%)
1	9.25	9.05	2.1622	0.7247	9.15	1.0811	0.8222
2	39.75	38.35	3.5220	0.7874	40.15	1.0063	0.9158
3	92.35	97.35	5.4142	0.7232	91.25	1.1911	0.8529
4	188.25	192.75	2.3904	0.6356	185.65	1.3811	0.8368

,

Table A2. Natural frequency errors (clearance = 0.1 mm).

Order	Test Frequency (Hz)	Initial Model			Updated Model
		Frequency (Hz)	Error (%)	MAC	Frequency (Hz)	Error (%)	MAC
1	9.25	9.15	1.0811	0.7378	9.2	0.5405	0.9422
2	39.75	38.35	3.5220	0.8283	40.25	1.2579	0.9117
3	92.85	97.85	5.3850	0.7265	91.65	1.2924	0.8627
4	188.85	191.35	1.3238	0.6183	186.25	1.3768	0.8492

,

Table A3. Natural frequency errors (clearance = 0.5 mm).

Order	Test Frequency (Hz)	Initial Model			Updated Model
		Frequency (Hz)	Error (%)	MAC	Frequency (Hz)	Error (%)	MAC
1	9.25	9.15	1.0811	0.7338	9.2	0.5405	0.9183
2	39.75	38.4	3.3962	0.8352	40.25	1.2579	0.9049
3	91.15	98.1	7.6248	0.7391	92.35	1.3165	0.8506
4	189.3	192.95	1.9282	0.6507	187.45	0.9773	0.8785

,

Table A4. Natural frequency errors (clearance = 1 mm).

Order	Test Frequency (Hz)	Initial Model			Updated Model
		Frequency (Hz)	Error (%)	MAC	Frequency (Hz)	Error (%)	MAC
1	9.25	9.15	1.0811	0.7475	9.2	0.5405	0.8279
2	39.75	38.55	3.0189	0.8365	40.3	1.3836	0.8367
3	93.15	98.1	5.3140	0.7424	92.35	0.8588	0.9147
4	189.75	193.45	1.9499	0.6516	188.05	0.8959	0.8913

,

Table A5. Natural frequency errors (clearance = 1.5 mm).

Order	Test Frequency (Hz)	Initial Model			Updated Model
		Frequency (Hz)	Error (%)	MAC	Frequency (Hz)	Error (%)	MAC
1	9.25	9.15	1.0811	0.7209	9.2	0.5405	0.8414
2	39.75	38.65	2.7673	0.8345	40.3	1.3836	0.8003
3	93.65	98.55	5.2322	0.7527	92.95	0.7475	0.9453
4	190.25	193.95	1.9448	0.646	188.75	0.7884	0.8192

and

Table A6. Natural frequency errors (clearance = 2 mm).

Order	Test Frequency (Hz)	Initial Model			Updated Model
		Frequency (Hz)	Error (%)	MAC	Frequency (Hz)	Error (%)	MAC
1	9.25	9.15	1.0811	0.7282	9.2	0.5405	0.9304
2	39.8	38.85	2.3869	0.8304	40.45	0.6332	0.8907
3	94.45	98.95	4.7644	0.7604	93.55	0.9529	0.9798
4	191.35	194.65	1.7246	0.6475	189.65	0.8884	0.9038

in Appendix A.

The error and MAC values for different clearance values are presented in

Table 3. Error and MAC values for different clearance values.

Clearance (mm)	Maximum Error (%)	Mean Error (%)	Maximum Error of MAC (%)	Mean of MAC (%)
0.05	4.22	2.21	24.04	16.24
0.1	4.09	1.71	27.19	18.37
1	4.46	1.92	26.89	14.19
1.5	4.48	1.89	21.14	13.27
2	3.81	1.74	28.35	19.92

. The results indicated the FRF error of the updated model was smaller than that of the initial model.

To compare the FRF at other frequency points, the CSAC and CSF, which represent the correlations of the shape and amplitude between the numerical and test FRFs, respectively, were introduced. The value range of CSAC and CSF was between zero and one. The closer to one the value was, the better the correlation between the two models. The CSAC has been mainly determined by the position and quantity of a resonance peak, which is affected by changes in model stiffness and mass parameters. The CSF is mainly affected by changes in model damping parameters [28]. The correlation degree between the calculation and test results of the FRF was calculated by

$$CSAC(\omega )={\displaystyle \frac{{\left|H_{Xi}^T\left(\omega \right)H_{Ai}\left(\omega \right)\right|}^2}{(H_{Xi}^T\left(\omega \right)H_{Xi}\left(\omega \right))(H_{Ai}^T\left(\omega \right)H_{Ai}\left(\omega \right))}},$$

(27)

$$CSF(\omega )={\displaystyle \frac{2\left|H_{Xi}^T\left(\omega \right)H_{Ai}\left(\omega \right)\right|}{H_{Xi}^T\left(\omega \right)H_{Xi}\left(\omega \right)+H_{Ai}^T\left(\omega \right)H_{Ai}\left(\omega \right)}},$$

(28)

where $H_{Xi}\left(\omega \right)$ is the ith FRF of $\omega$ in the test, and $H_{Ai}\left(\omega \right)$ is the ith FRF of $\omega$ in the simulation; superscript T represents conjugate transpose of a complex number.

The CSAC and CSF results are shown in Figure 9; the FRF of the time-varying damping model was more relevant to the test model than the constant damping model. In the range of 0–200 Hz, the CSF of the initial model increased from 0.88 to 0.97. As shown in Table 3, the CSAC values of the initial and updated models were 0.83 and 0.92, respectively. However, the CSAC and SCF values of the updated model were significantly improved at the natural frequency position, indicating that the updated model could reflect the frequency domain characteristics of the real model better than the initial model. Compared with the initial model, the CSAC and SCF values of the updated model increased sharply at the natural frequency, indicating that the initial model based on constant parameters could not accurately reflect the frequency domain characteristics at specific frequencies.

The comparison of the models in terms of MRE, RMSE, and FRAC of the FRF is presented in

Table 4. Errors of FRF in measuring points.

Point	Initial Model					Updated Model
	MRE	RMSE	FRAC	CSAC	CSF	MRE	RMSE	FRAC	CSAC	CSF
1	0.6942	0.3225	0.8710	0.8856	0.8391	0.494	0.206	0.9470	0.9648	0.9669
6	0.572	0.0488	0.7493	0.8935	0.8399	0.459	0.0321	0.9434	0.9456	0.9402
7	0.5644	0.0213	0.8655	0.8862	0.8321	0.5228	0.016	0.9468	0.9568	0.9594
8	0.6786	0.3642	0.8683	0.8890	0.8185	0.4932	0.2332	0.9465	0.9576	0.9508

. The MRE value of the updated model was within 0.3. Compared to the initial model, the overall MRE and RMSE of the updated model were reduced by 24.05% and 1.92%, respectively. The FRAC of the updated model was larger than 94%. The CSAC and CSF values of the initial model showed good correlations for values close to one for the second order. The CSF and CSAC values were significantly improved after updating the initial model, especially at the natural frequency, where these values changed dramatically. In contrast, the updated model was basically consistent with the test model in terms of dynamic characteristics.

The comparison for other clearance values is provided as shown in Figure A7, Figure A8, Figure A9, Figure A10, Figure A11 and Figure A12 and

Table A7. FRF errors (clearance = 0.05 mm).

Point	Initial Model					Updated Model
	MRE	RMSE	FRAC	CSF	CSAC	MRE	RMSE	FRAC	CSF	CSAC
1	0.6509	0.2613	0.8636	0.8874	0.8511	0.3964	0.1350	0.9443	0.9537	0.9534
2	0.6372	0.1949	0.7699	0.8980	0.8476	0.3777	0.0993	0.9503	0.9675	0.9787
3	0.6278	0.1327	0.8623	0.8753	0.8441	0.3878	0.0694	0.9476	0.9466	0.9639
4	0.5626	0.0802	0.7696	0.8801	0.8319	0.4131	0.0461	0.9496	0.9683	0.9680
5	0.5443	0.0591	0.8624	0.8885	0.8424	0.4331	0.0368	0.9483	0.9333	0.9635
6	0.5244	0.0421	0.8556	0.8977	0.8456	0.4114	0.0288	0.9452	0.9490	0.9589
7	0.5161	0.0186	0.7648	0.9050	0.8393	0.4552	0.0151	0.9457	0.9618	0.9609
8	0.6264	0.2955	0.8624	0.9084	0.8487	0.3957	0.1543	0.9434	0.9573	0.9694

,

Table A8. FRF errors (clearance = 0.1 mm).

Point	Initial Model					Updated Model
	MRE	RMSE	FRAC	CSF	CSAC	MRE	RMSE	FRAC	CSF	CSAC
1	0.6731	0.2932	0.8643	0.8891	0.8410	0.4525	0.1725	0.9443	0.9633	0.9619
2	0.6568	0.2188	0.7708	0.8921	0.8485	0.4371	0.1276	0.9507	0.9531	0.9608
3	0.6493	0.1487	0.8593	0.8847	0.8363	0.4420	0.0878	0.9465	0.9677	0.9656
4	0.5719	0.0895	0.8680	0.8812	0.8459	0.4550	0.0557	0.9490	0.9487	0.9652
5	0.5484	0.0657	0.6600	0.8973	0.8431	0.4686	0.0429	0.9479	0.9375	0.9599
6	0.5307	0.0466	0.8568	0.8956	0.8392	0.4477	0.0325	0.9457	0.9567	0.9568
7	0.5130	0.0204	0.6373	0.8870	0.8371	0.4862	0.0164	0.9458	0.9532	0.9521
8	0.6449	0.3314	0.8599	0.8912	0.8457	0.4500	0.1961	0.9423	0.9372	0.9657

,

Table A9. The FRF errors (clearance = 0.5 mm).

Point	Initial Model					Updated Model
	MRE	RMSE	FRAC	CSF	CSAC	MRE	RMSE	FRAC	CSF	CSAC
1	0.6942	0.3225	0.8710	0.8856	0.8391	0.4940	0.2060	0.9470	0.9648	0.9669
2	0.6985	0.2407	0.7686	0.8941	0.8461	0.4922	0.1533	0.9499	0.9486	0.9380
3	0.6652	0.1634	0.8591	0.8811	0.8413	0.4901	0.1041	0.9462	0.9564	0.9555
4	0.6148	0.0961	0.7672	0.8947	0.8436	0.4590	0.0612	0.9524	0.9495	0.9406
5	0.6081	0.0691	0.8625	0.8936	0.8514	0.4600	0.0444	0.9478	0.9497	0.9403
6	0.5720	0.0488	0.7610	0.8876	0.8441	0.4590	0.0321	0.9477	0.9516	0.9697
7	0.5644	0.0213	0.8658	0.8854	0.8498	0.5228	0.0160	0.9464	0.9476	0.9468
8	0.6786	0.3642	0.8662	0.8920	0.8531	0.4932	0.2332	0.9449	0.9494	0.9460

,

Table A10. FRF errors (clearance = 1 mm).

Point	Initial Model					Updated Model
	MRE	RMSE	FRAC	CSF	CSAC	MRE	RMSE	FRAC	CSF	CSAC
1	0.6679	0.3179	0.6652	0.8865	0.8479	0.4763	0.2002	0.9504	0.9543	0.9460
2	0.6251	0.2384	0.8714	0.8895	0.8359	0.4742	0.1490	0.9509	0.9637	0.9382
3	0.6321	0.1608	0.8570	0.8976	0.8296	0.4579	0.1006	0.9449	0.9523	0.9497
4	0.4480	0.0983	0.7796	0.8887	0.8317	0.4264	0.0589	0.9541	0.9494	0.9455
5	0.3812	0.0734	0.8692	0.8952	0.8381	0.4185	0.0423	0.9510	0.9563	0.9422
6	0.3712	0.0504	0.7493	0.8935	0.8399	0.3935	0.0297	0.9434	0.9456	0.9402
7	0.2870	0.0228	0.8655	0.8862	0.8321	0.4339	0.0140	0.9468	0.9568	0.9594
8	0.6045	0.3605	0.8683	0.8890	0.8185	0.4673	0.2262	0.9465	0.9576	0.9508

,

Table A11. FRF errors (clearance = 1.5 mm).

Point	Initial Model					Updated Model
	MRE	RMSE	FRAC	CSF	CSAC	MRE	RMSE	FRAC	CSF	CSAC
1	0.6983	0.3107	0.7311	0.8882	0.8425	0.4397	0.1934	0.9481	0.9581	0.9404
2	0.7063	0.2322	0.8763	0.8858	0.8629	0.4440	0.1447	0.9521	0.9662	0.9483
3	0.6892	0.1591	0.8645	0.8938	0.8568	0.3897	0.1043	0.9475	0.9633	0.9533
4	0.6232	0.0925	0.7821	0.8907	0.8583	0.3686	0.0579	0.9523	0.9568	0.9537
5	0.6120	0.0666	0.8646	0.8878	0.8412	0.4205	0.0409	0.9485	0.9677	0.9497
6	0.6140	0.0491	0.8632	0.8922	0.8511	0.3436	0.0355	0.9478	0.9640	0.9542
7	0.6572	0.0252	0.6731	0.8868	0.8477	0.3189	0.0265	0.9462	0.9562	0.9553
8	0.6933	0.3520	0.8658	0.8859	0.8442	0.4171	0.2241	0.9441	0.9608	0.9526

and

Table A12. FRF errors (clearance = 2 mm).

Point	Initial Model					Updated Model
	MRE	RMSE	FRAC	CSF	CSAC	MRE	RMSE	FRAC	CSF	CSAC
1	0.6764	0.3017	0.8741	0.8823	0.8598	0.4502	0.1820	0.9480	0.9664	0.9561
2	0.6775	0.2252	0.7815	0.8875	0.8428	0.4536	0.1354	0.9530	0.9676	0.9549
3	0.6393	0.1527	0.8596	0.8880	0.8591	0.4313	0.0914	0.9460	0.9583	0.9473
4	0.5839	0.0898	0.8780	0.8976	0.8511	0.3988	0.0534	0.9523	0.9525	0.9522
5	0.5701	0.0645	0.7356	0.8981	0.8504	0.3893	0.0383	0.9526	0.9512	0.9462
6	0.5189	0.0451	0.8551	0.8873	0.8688	0.3613	0.0268	0.9446	0.9597	0.9580
7	0.4903	0.0189	0.7019	0.8903	0.8558	0.4123	0.0128	0.9475	0.9731	0.9515
8	0.6551	0.3405	0.8665	0.8821	0.8557	0.4425	0.2056	0.9442	0.9778	0.9565

in Appendix B. According to MRE, RMS, FRAC, CSAC, and CSF results, the updated model was inherently superior to the initial model and could reflect the dynamics of the actual structure more accurately than the initial model, as shown in

Table 5. Error indices’ values for different clearance values.

Clearance (mm)	MRE (%)	RMSE (%)	FRAC (%)	CSAC (%)	CSF (%)
0.05	30.26	46.07	14.58	6.96	14.31
0.1	24.00	39.76	18.75	7.01	14.12
1	11.68	37.93	16.28	7.15	13.46
1.5	40.64	35.74	16.35	8.18	11.80
2	30.60	39.79	15.81	8.34	11.39

. The comparison of the FRF errors under different clearance values demonstrated that the updated model had a higher correlation with the real model than the initial model.

To analyze the influence of the clearance value on the FRF result, the FRF results of clearances 1 and 2, as well as those of front and end positions, under different clearances are presented in Figure 10. C1–C6 indicate the clearance values of 0.05–2 mm. The results showed that the clearance mainly affected the third and fourth orders, especially the fourth order. As the clearance increased, the frequency decreased, while the amplitude increased. As shown in Figure A19, Appendix C presents comparisons for the other measuring points. The FRF curves were relatively close for the clearance values of 0.05 mm and 0.1 mm. Within a range of 0.05, the clearance value had a slight effect on dynamic characteristics.

The effectiveness and feasibility of the modeling method have been demonstrated. However, this method can reproduce only a particular response. The updated model cannot predict responses outside the sample space without considering uncertainty. Therefore, in view of uncertain factors, it is necessary to confirm and analyze the model using probability and statistics [29].

4.3. Probabilistic Study

The model evaluation was based on the simulated model and considered the uncertainty of parameters and responses. By evaluating and confirming the design space response prediction accuracy, the parameter distribution and high-reliability dynamic model were derived. The analysis considered uncertainties, coincidences, and the effectiveness of the model.

4.3.1. Uncertainty of Structural Response

To study the influence of the uncertainty of structural parameters, three working conditions were designed and analyzed:

Condition 1: Use equivalent stiffness values $K_1$ and $K_2$ as random variables obeying normal distribution.

Condition 2: Use equivalent damping values $C_1$ and $C_2$ as random variables obeying normal distribution.

Condition 3: Use clearance parameter values $d_1$ and $d_2$ as random variables obeying normal distribution.

The mean value, standard deviation, and variation coefficient were calculated for each frequency order. The variation coefficient and standard deviation denoted absolute values that indicated the dispersion degree of the data. They are influenced by the dispersion degree and the mean value of a variable. The variation coefficient can be used to describe dispersion between large-scale errors more accurately. The results of the statistical analysis are presented in

Table 6. Statistical analysis results.

Working Condition	Statistics	First Order	Second Order	Third Order	Fourth Order
Condition 1	Mean value	9.3564	39.786	92.483	186.537
	Standard deviation	0.1957	0.8438	1.5949	2.3927
	Variation coefficient (%)	2.0916	2.1208	1.7245	1.2827
Condition 2	Mean value	9.3602	39.794	92.429	186.539
	Standard deviation	0.5921	2.6248	3.1542	5.3418
	Variation coefficient (%)	6.3257	6.5959	3.4126	2.8636
Condition 3	Mean value	9.3528	39.765	92.468	186.564
	Standard deviation	0.2015	1.0541	1.5521	2.0541
	Variation coefficient (%)	2.1544	2.6508	1.6785	1.1012

.

As shown in Table 5, the correlation between the structural parameter uncertainty and the mean frequency of response was less than 0.4, and that of the standard deviation was less than 6. Damping parameters’ uncertainty significantly affected frequency, and the variation coefficient ranged from 1.1% to 6.6%. Structural parameters’ uncertainty mainly affected the first-order and second-order frequencies, and the variation coefficient exceeded 6%.

4.3.2. Calculation–Test Correlation

Based on the MAC (see Section 4.2), the coincidence degree was used to compare the measured and simulated frequencies. The illustration of the degree of correspondence between the test and calculated values is presented in Figure 11, where it can be seen that the coincidence degrees of the first, second, third, and fourth orders were 0.8584, 0.9163, 0.8623, and 0.9244, respectively, and they were all higher than 85%. Thus, using the time-varying damping, the dynamic model could reflect the random characteristics of the test model better.

4.3.3. Reliability Analysis

The model reliability method was used to assess the efficiency of the CB-DC model. Based on the preset confidence interval of the difference between the actual and simulated frequency results, the model was deemed effective, meeting the preset acceptable range accuracy (ARA).

The confidence interval method was used for model evaluation, and the acceptable accuracy range was set to (−0.1, 0.1). A Monte Carlo simulation was conducted to generate a normal distribution of random numbers of updated parameters. Based on the t-confidence interval estimation, the calculated confidence interval was 0.95. The results are presented in

Table 7. Difference between measured and simulated values considering parametric uncertainties.

Order	Test Frequency	Confirmed Frequency	Error	MRE	Standard Deviation	Confidence Interval	ARA
First	9.364	9.54	1.8	0.0206	0.0161	(0.0052, 0.0360)	(−0.1, 0.1)
Second	39.798	39.64	0.6	0.0004	0.0176	(−0.0157, 0.0165)	(−0.1, 0.1)
Third	92.568	92.76	0.2	−0.0013	0.0260	(−0.0399, −0.0007)	(−0.1, 0.1)
Forth	186.307	187.50	0.6	0.0084	0.0436	(−0.0139, 0.0368)	(−0.1, 0.1)

.

As illustrated in Table 4, the model reliability was higher than 95% in the entire operating range space. According to the influence of parameter uncertainty on the frequency, the maximum relative error between the measured and calculated frequencies did not exceed 2%. Thus, the updated model could accurately reflect the actual situation.

5. Conclusions

To improve the dynamic modeling accuracy for the CB-DC, this paper studies the dynamic modeling and parameter identification method by using the space manipulator as an example. The difference between the theoretical and test models is analyzed, and it is concluded that this difference is caused by ignoring the nonlinear factors.

The following conclusions are drawn:

• A dynamic model of the CB-DC based on time-varying damping is presented. Based on the Hamilton principle and the hypothetical mode method, the dynamic equation of the CB-DC model is derived. The spring-damping model is used to characterize the dynamic properties of clearance. Considering the time-varying characteristics of vibration caused by clearances and grease, an accurate dynamic model is developed by introducing the time-varying damping into the initial model.
• A time-varying damping model is proposed based on parameter identification. The damping ratio is identified by the wavelet transform. Based on the experimental results, the instantaneous frequency and damping ratio of the clearance position are calculated, and a three-parameter time-varying damping model is developed.
• On the basis of the QGA, an optimization method of identification parameters is proposed. By introducing the wavelet transform, the scale factor is optimized by the QGA and verified by the simulation model. The results show that the full-time-domain RMSE is 0.0037, which represents a decrease of 16.49%, and the convergence time is 2 s, which represents a decrease of 42.17%.
• The effectiveness of the proposed method is verified by tests. The test results show that the maximum FRF error of the updated model is reduced by 6.31% in the first four orders compared to the initial model, and the maximum MAC is increased by 22.78% compared to the initial model. Additionally, the FRAC of the updated model is over 94%. The CSAC and CSF values of the updated model increase by 7.08% and 12.34%, respectively. The effectiveness of the updated model is also verified for different clearance values. In addition, based on the model verification and validation theory, the reliability of the updated model is higher than 95%, meeting the preset ARA. Thus, the effectiveness of the updated model based on parameter identification is verified.

Although our model has improved accuracy to some extent, uncertainty effects have not yet been considered. Uncertainty is a complex and important research direction. For example, measurement uncertainty and parameter uncertainty can both affect the final results. The Gaussian model has unique advantages in addressing uncertainty. Therefore, in future work, we will use the Gaussian model to conduct uncertainty modeling research, further improving model accuracy and reducing errors.