Frequency-Dependent Bouc–Wen Modeling of Magnetorheological Damper Using Harmonic Balance Approach

Abstract

Magnetorheological dampers (MRDs) are of great interest in engineering due to their continuously adjustable damping characteristics. Accurate models are essential for optimizing MRDs and analyzing system dynamics. However, conventional methods widely overlook the impact of excitation frequency and amplitude. To address this issue, this work proposes a modified Bouc–Wen model that can be adapted to various excitation conditions. The model’s parameters depend on the current, excitation frequency, and amplitude. The mechanical characteristics of the MRD were analyzed by the tests. The parameters in the Bouc–Wen model were identified by combining the harmonic balance method and the genetic algorithm. The modified Bouc–Wen model was established by analyzing the variation of each parameter with current, excitation frequency, and amplitude. Finally, the agreement between the modified prediction model and the test results was verified under sinusoidal excitation of 80 mm and 1 Hz. The average relative errors were 3.87%, 2.82%, 2.45%, 2.19%, and 3.27% for current excitations of 0 A, 0.5 A, 1 A, 1.5 A, and 2.0 A, respectively. Since the MRD in this paper operates from 0.5 Hz to 2 Hz, the modified model was validated in the same range. Experiments demonstrate that the modified Bouc–Wen model efficiently and accurately describes the mechanical properties of the MRD under various excitation conditions.

1. Introduction

Magnetorheological (MR) fluids are a new type of smart material that exhibits controllable rheological properties when subjected to an applied magnetic field [1]. These fluids have attracted considerable attention in various fields, including vehicle suspension, construction, and bridges [2] due to their continuously adjustable damping characteristics [3], fast response time, and low energy consumption [4]. However, accurately and reliably modeling the dynamics of MRDs is challenging due to their hysteresis and nonlinear characteristics [5]. Several models, such as the Bingham model [6], the Bouc–Wen model [7], and the hyperbolic tangent model [8], have been proposed to describe the mechanical properties of MRDs. The Bouc–Wen model effectively captures the nonlinear characteristics of the dynamic properties of MRDs at low speeds [9]. It is widely used in the nonlinear modeling of MRDs due to its ability to describe hysteresis curves of various shapes by adjusting the parameters [10]. However, the main challenge is to identify the characteristic parameters reasonably and effectively. Traditional methods often base their results on a particular excitation condition. When external excitation is altered, the model error can often be significant, requiring re-identification of parameter values.

Previous model analyses have typically treated model parameters as a function of the applied current. However, in reality, the damping output of a magnetorheological damper grows with increasing excitation frequency [11]. Therefore, it is necessary to propose a frequency-dependent dynamic model for MRDs, where the model parameters are functions of current, excitation frequency, and amplitude. The dynamics of the magnetorheological damper under different excitations can be predicted by considering the effects of the excitation frequency, amplitude, and current. Recently, numerous theoretical models have been proposed to predict the damping characteristics of MRDs. Dominguez [12] proposed a new Bouc–Wen model with current, amplitude, and frequency as variables. The new model can better represent hysteresis behavior under different working conditions. The new model adds to the original model a mathematical expression with excitation frequency and amplitude as independent variables. However, there are 16 parameters to be identified, which may decrease the efficiency. Ali [13] proposed a frequency-dependent MRD model based on the Spencer model. However, the model requires the identification of a total of 10 parameters, making the calculation less efficient. Boada [14] proposed a frequency-dependent MRD model based on trigonometric functions. The model utilizes a neural network (NN) to model the relationship between the parameters and frequency, which has fewer parameters and highly accurate model results. Nevertheless, this method may require large quantities of test samples and long-term training. Lei [15] established a predictive model by piecewise describing the damping characteristics of the MRD. The predictive model shows great agreement with experiments under different currents. However, there are some limitations when it comes to different amplitudes and frequencies. Lv [16] proposed a predictive MRD dynamics model based on the analysis of the rheological behavior of MR fluids at high shear rates. The predictive model performs well in describing the characteristics of MRDs, except for hysteresis.

Numerous strategies related to neural networks and machine learning have been used for modeling the dynamics of magnetorheological dampers. These methods require a large amount of data to be analyzed, and the computational cycle time is too long to obtain an accurate model in a short period of time. To describe the mechanical properties of MRDs under different excitation conditions efficiently, a modified model is necessary, as the traditional model requires re-identification when the excitation conditions change. The Bouc–Wen model explains the stiffness, kinematic viscosity, and other properties of the damper through various parameters that are related to the current, frequency, and amplitude of the excitation. This paper presents a dynamic model for MRDs using current, frequency, and amplitude as variables. The modified model accurately predicts the mechanical properties of MRDs under different external excitations.

This paper proposes a model based on the Bouc–Wen model that depends on current, frequency, and amplitude. The model predicts the output damping force of the MRD based on the excitation amplitude, frequency, and current. The model has fewer parameters to identify, ensuring high calculation efficiency. The MRD discussed in this paper is applied to vehicle suspension damping, with an operating frequency range of 0.5 Hz to 2 Hz. Since the main objective of the MRD in this paper is to improve the driving comfort of the vehicle by attenuating vibrations, the prediction model developed in this paper only applies frequencies up to 2 Hz, matching the natural frequencies of the body model. To achieve these objectives, the paper first conducted experiments to obtain the mechanical characteristic curve of the MRD. The analysis shows that the hysteresis force of the damper is influenced by the current, excitation frequency, and amplitude. Afterward, the MRD was modeled using the Bouc–Wen model. The harmonic balance method was used in combination with the genetic algorithm to determine specific parameter values under different excitation conditions. Then, the functions were established by analyzing the variation of each parameter with current, excitation frequency, and amplitude. A modified Bouc–Wen model is proposed to describe the mechanical properties of MRDs under various excitation conditions. The model includes current, excitation frequency, and amplitude as variables. Finally, the model was used to predict the mechanical properties of this damper under 80 mm and 1 Hz sinusoidal excitation. The simulations obtained by this method, as well as the other two methods, are compared with the experimental data. The findings suggest that the modified model can accurately predict the hysteresis force of the damper under different excitation conditions with an acceptable level of error.

2. MRD Mechanical Property Tests

To create a nonlinear model of the MRD, it is essential to analyze its dynamic properties. The MRD used in the test was developed by the group. Fatigue testing machine input displacement excitation to the damper, and the transducer was used to obtain the output damping force. Figure 1 shows the MRD used in the test, which consists of a single outlet rod and double barrel construction. The primary and secondary barrels are cemented to each other using nuts. The primary barrel is an MR module that provides the damping force, while the secondary barrel has a pneumatic pocket to provide volume compensation. The test system is shown in Figure 2. It comprises a fatigue testing machine, a controlled current source, and a control computer. The MRD is mounted on the fatigue tester, and the control computer adjusts the amplitude and frequency of the sinusoidal excitation as needed. A current source is connected to the damper to provide controlled current. The displacement and velocity signals collected by the sensors connected to fatigue machines are used by the control computer to generate corresponding images.

The effect of friction on the damping force in the MRD should be analyzed prior to the mechanical property test. The MRD containing no fluid was tested by applying a small uniform stretch (0.01 mm/s), and the friction force was obtained, as shown in Figure 3. The friction fluctuates from 25 to 100 N, which is so small compared to the overall force that it can be completely ignored. This is due to the fact that the friction of the damping force in this paper comes mainly from the two seals, the end cap and the piston. Since the number of seals is small, the friction force is tiny.

Subsequently, two sets of tests were performed on the magnetorheological damper using a sinusoidal signal as input displacement excitation.

$$x=A\cdot \sin \left(2\pi \cdot f\cdot t\right)$$

(1)

where $x$ denotes the input displacement in millimeters, $A$ denotes the amplitude of the sinusoidal signal in millimeters, and $f$ denotes the frequency of a sinusoidal signal in Hz.

The MR damper was subjected to current excitations in the range of 0–2 A for different combinations of the sinusoidal signal frequencies (0.5, 1.0, 1.5, 2.0 Hz) and amplitudes (20, 50 mm).

Figure 4 shows that for the same current and displacement, the damping output of the magnetorheological damper varies at different frequencies. As shown in Figure 4a, an increase in excitation frequency results in an increase in both the maximum damping output and the enclosed area of the curve. This is because the high-frequency excitation signal has a larger instantaneous velocity, and the damping force is positively related to the velocity. The area enclosed by the force–displacement curve represents the energy dissipated by the damper during each vibration cycle. The area increases with frequency, indicating that the energy dissipated by the damper also increases. Figure 4b shows that the amplitude of the damping force at different frequencies for the same velocity is basically equal when the velocity and acceleration are in the same direction. However, when the velocity and acceleration are in the opposite direction, the damping force at the same velocity gradually decreases with increasing frequency.

Figure 4a shows a non-standard rectangular curve with a significant bulge in the region where the displacement is close to zero. This bulge is due to the elastic forces generated by the air spring used to compensate for the volume in the damper. Furthermore, the curve exhibits a discontinuity in the form of a step at the point where displacement reaches its maximum value. This corresponds to the region in Figure 4b, where velocity approaches zero and then reverses with acceleration. This is due to the fact that insufficient pressure in the air pocket prevents the flow of the MR fluid from keeping up with the movement of the piston, creating a vacuum cavity in the working area. The absence is more noticeable at higher frequencies and velocity amplitudes of the damper.

Figure 4b shows that the force–velocity curve of the MRD has a larger hysteresis loop in the low-speed region and a smaller hysteresis loop in the high-speed region. This is because the damping properties of MR fluid vary with speed, exhibiting elastic properties at low speed and plastic properties at high speed. Therefore, the low-velocity region of the force–displacement curve is referred to as the elastic region, and the high-velocity region is referred to as the plastic region.

Figure 5 and Figure 6 demonstrate that the output damping force of the magnetorheological damper increases with the input current. At 1.5 A, the damping output force gradually tends to saturate due to the magnetorheological fluid reaching its magnetic saturation strength. Simultaneously, as the current increases, the area of the force–displacement curve envelope also increases in both Figure 5 and Figure 6. This suggests that the damper’s energy dissipation also increases with the current.

In summary, the mechanical properties of the magnetorheological damper are not only related to the current but also to the frequency and amplitude of the excitation conditions. This is because the variation in frequency and amplitude affects the speed of movement of the piston.

3. Bouc–Wen Model for MRD Considering Inertial Effects

The Bouc–Wen model of MRD typically comprises a spring element, a linear damping component, an equivalent mass block, and a hysteresis system that characterizes the hysteresis properties [17], as shown in Figure 7.

The equations for the dynamics of the magnetorheological damper, taking into account inertial effects, are established using the Bouc–Wen model [18], as shown in Equation (2).

$$\left\{\begin{array}{l}{F}_d=k\cdot x+c\cdot \dot{x}+\alpha \cdot z+f_0+m\cdot \ddot{x}\\ \dot{z}=\rho \cdot \left[\dot{x}-\sigma \cdot \left|\dot{x}\right|\cdot {\left|z\right|}^{n-1}\cdot z+\left(\sigma -1\right)\cdot \left|\dot{x}\right|\cdot {\left|z\right|}^n\right]\end{array}\right.$$

(2)

where $F_d$ denotes the damping output force, $k\cdot x$ denotes the elastic component of the accumulator in the magnetorheological damper, $c$ denotes the viscosity coefficient of the magnetorheological material after yielding, $x,\dot{x},\ddot{x}$ denote the excitation displacement, velocity, and acceleration of the system, $\alpha$ denotes the adjustment coefficient of the hysteresis force to the total damping output force, $f_0$ denotes the zero-point drift of the damping force at the initial position, $m$ denotes the equivalent mass of the inertial effect, $z$ denotes the hysteresis displacement, $\dot{z}$ denotes the derivative of the hysteresis displacement with respect to time, and $\rho ,\sigma ,n$ are the parameters controlling the shape of the hysteresis loop, and the system is passively stabilized when $\rho >0,\sigma >0.5$.

4. Parameter Identification Method for the Bouc–Wen Model of MRD

The Bouc–Wen model can effectively describe the hysteresis characteristics of magnetorheological dampers, but it contains numerous unknown parameters and exhibits strong nonlinearity. Directly using the genetic algorithm to identify the model parameters may result in identification that matches the test data but still has some issues, as follows:

• The identification results do not take into account the physical significance of the parameters and disregard the coupling between the parameters and the input frequency and amplitude.
• The model has an excessive number of unknown parameters, leading to a prolonged identification process.
• The model identified is only suitable for a single working condition. It is necessary to identify a new model for other working conditions.

Therefore, an efficient parameter identification method is proposed. In this process, parameters with physical significance are analyzed by other methods, while the remaining parameters are identified using the genetic algorithm.

For a given displacement input, the Bouc–Wen model requires identification of eight parameters. The sensitivity analysis of these parameters found that $\sigma ,n$ were less sensitive and had less impact on the model output [19]. During the identification process, the initial setting can be established as $\sigma =1,n=2$. Therefore, the parameters to be identified are reduced to 6. $k,c,m$ have clear physical meanings. $\alpha ,\rho$ describe the various patterns of hysteresis curves. Therefore, the parameter identification process should be separated, and the physical meaning of the specific parameter should be analyzed.

4.1. Identification of $k$

The air spring provides the elastic force during the operation of the damper. This force is typically expressed as the product of the elastic coefficient and the displacement in previous models of MRDs. However, these models treat the air spring as an elastic medium with constant stiffness, which is not entirely accurate. In reality, the air spring exhibits an asymptotic stiffness property, which means the stiffness increases progressively with displacement. To ensure modeling accuracy, it is necessary to consider stiffness asymptotic properties when modeling elastic forces.

According to ideal gas state equation [20]:

$$P\cdot V^r=const$$

(3)

When the floating piston is in any position, the volume of the gas chamber can be expressed as [21]:

$$V_1=V_0+A_S\cdot x$$

(4)

where $x$ denotes the displacement of the piston; $A_S$ denotes the cross-sectional area of the piston rod; $A_S\cdot x$ denotes the compensation volume provided by the air spring; $V_1$ denotes the current volume of the gas chamber of the air spring; $V_0$ denotes the volume of the gas chamber in static equilibrium.

The expression for the pressure of the gas chamber when the floating piston is in any position can be obtained as follows [21]:

$$P_1={\displaystyle \frac{P_0\cdot V_0^r}{{\left(V_0+A_S\cdot x\right)}^r}}$$

(5)

where $P_1$ denotes the pressure of the air spring in its current state; $P_0$ denotes internal pressure of the gas spring in its initial state; $r$ denotes the gas polytropic exponent. Since the process of filling and deflating the gas chamber can be approximated as an adiabatic process, $r$ can be set as 1.4.

Equation (5) includes fractional order terms, which can complicate the analysis and solution of the system response. Set $h={\displaystyle \frac{V_0}{A_S}}$, and Taylor expands the above Equation (5) at $x$ as follows:

$$F_S=P_1\cdot A_S=P_0\cdot A_S+{\displaystyle \sum _{i=0}^n{k}_i\cdot x^i}$$

(6)

where $k_i$ denotes the stiffness coefficient of the air spring, and $i$ denotes the order of expansion.

$$k_i=P_0\cdot A_S\cdot {\displaystyle \frac{{\displaystyle \prod _{m=1}^i\left(-r-m+1\right)}}{h^i\cdot i!}}$$

(7)

The accuracy and complexity of the solution are affected by the value of $i$. When the expansion contains a few terms, the results obtained may not accurately reflect the nonlinear characteristics of the system. However, when the expansion contains numerous terms, the solution process becomes significantly more difficult. By taking the first three orders of Equation (6) as the system stiffness, the elastic force can be expressed as follows:

$$F_S=P_0\cdot A_S+k_1\cdot x+k_2\cdot x^2+k_3\cdot x^3=k\cdot x$$

(8)

The system stiffness is determined by selecting the first three orders. This is because including primary and tertiary terms in the stiffness expression ensures the asymptotic stiffness characteristic of the air spring.

4.2. Identification of $m,c,f_0$

When addressing nonlinear vibration problems, numerical analysis can provide an approximate solution for the system. The harmonic balance method is a concise and user-friendly numerical analysis method that can solve both weak and strong nonlinear systems [22]. This method expands the input excitation and response of the system into the same level of the Fourier series. The coefficients of the harmonic terms of the same order on both sides of the equation must be equal. The coefficients of the Fourier series can be determined by solving the series of equations obtained.

The equation for the nonlinear system under forced vibration can be expressed in universal form as follows:

$$\ddot{x}+f\left(x,\dot{x}\right)=F\left(t\right)$$

(9)

Assuming that $F(t)$ is an even function and does not contain a constant subterm, $f(x,\dot{x})$ is a nonlinear term containing the elastic and damping forces of the system. When the system is in periodic motion ($T=2\pi /w$), the right-hand side of the Equation (9) can be expanded into a Fourier series with period $T$ as follows:

$$F(t)={\displaystyle \sum _{n=1}^{\infty }{f}_n\cos (nwt)}$$

(10)

where

$$f_n={\displaystyle \frac{1}{T}}{\displaystyle {\int }_{-T/2}^{T/2}F(t)\cos (nwt)dt}\hspace{1em}(n=1,2,\cdots )$$

(11)

Expanding the left side of Equation (9) to a Fourier series with period $T$ as follows:

$$x(t)=A_0+{\displaystyle \sum _{n=1}^{\infty }[A_n\cos (nwt)+B_n\sin (nwt)]}$$

(12)

By substituting Equations (10) and (12) into Equation (9), the equation that contains the harmonics with the order of $n$ can be obtained. The coefficients preceding the harmonics of the same order are strictly equal. Therefore, the equation for $n$ sets of harmonic coefficients can be obtained. By solving the equation of finite order, the relationship between the amplitude and frequency of the harmonics for each order of the system can be obtained [23].

Periodic inputs to the Bouc–Wen model result in outputs of the same period [24]. Therefore, the analytical solution can be obtained by using the harmonic balance method for the Bouc–Wen model.

This paper presents a set of equations that contain the model parameters. These equations are obtained by processing the inputs and outputs of the Bouc–Wen model using the harmonic balance method. The purpose is to establish the mathematical relationships between the model parameters.

The first order harmonic balance method is used to analyze the Bouc–Wen model. Firstly, the hysteresis term equation containing the hysteresis displacement $z$ is processed. $z$ is expressed as follows:

$$\dot{z}=\rho \cdot \dot{x}\cdot [1-z^2\cdot (\sigma \cdot \mathrm{sgn}\left(\dot{x}\right)\cdot \mathrm{sgn}\left(z\right)+\left(\sigma -1\right)\cdot \mathrm{sgn}\left(\dot{x}\right)]$$

(13)

The equation includes a sign function. Therefore, the sign of the velocity and hysteresis displacement during the input variation with time needs to be considered when using the harmonic balance method of processing, which make it difficult to obtain accurate parametric equations.

It is proposed that the hysteresis curve of the Bouc–Wen model has a plastic region when the elastic force term is removed [25]. In the plastic region, there exists $\left|z\right|\approx 1$ in the Bouc–Wen model. Therefore, the expression can be obtained as follows:

$$F=\left\{\begin{array}{l}k\cdot x+c\cdot \dot{x}+m\cdot \ddot{x}+\alpha +f_0\\ k\cdot x+c\cdot \dot{x}+m\cdot \ddot{x}-\alpha +f_0\end{array}\right.\hspace{1em}\begin{array}{l}z>0\\ z<0\end{array}$$

(14)

By substituting Equation (8) into Equation (14), the following is obtained:

$$F=P_0\cdot A_S+k_1\cdot x+k_2\cdot x^2+k_3\cdot x^3+c\cdot \dot{x}+m\cdot \ddot{x}\pm \alpha +f_0$$

(15)

The expressions of the Fourier expansion for $x$ and $F$ are as follows:

$$\left\{\begin{array}{l}F=F_1\cos wt+F_2\sin wt=F_0\sin (wt+{\theta }_0)\\ x=x_1\cos wt+x_2\sin wt=x_0\sin (wt+{\theta }_x)\end{array}\right.$$

(16)

When the input signal x to the system is a standard sinusoidal signal that determines the excitation and frequency, the value of ${\theta }_x$ is zero. By substituting Equation (16) into Equation (15), the following is obtained:

$$\begin{array}{l}{F}_0\cdot \sin (w\cdot t+{\theta }_0)=P_0\cdot A_S+k_1\cdot x_0\cdot \sin (w\cdot t)+k_2\cdot x_0^2\cdot {\sin }^2(w\cdot t)+\\ k_3\cdot x_0^3\cdot {\sin }^3(w\cdot t)+c\cdot x_0\cdot w\cdot \cos \left(w\cdot t\right)-m\cdot x_0\cdot w^2\cdot \sin \left(w\cdot t\right)\pm \alpha +f_0\end{array}$$

(17)

Expanding Equation (17), the coefficients in front of $\sin wt$ and $\cos wt$ are equal. The following expression can be obtained:

$$\left\{\begin{array}{l}{F}_0\cdot \cos {\theta }_0\cdot \sin \left(w\cdot t\right)=\left(k_1\cdot x_0+{\displaystyle \frac{3}{4}}\cdot k_3\cdot x_0^3-m\cdot x_0\cdot w^2\right)\cdot \sin \left(w\cdot t\right)\\ F_0\cdot \sin {\theta }_0\cdot \cos \left(w\cdot t\right)=c\cdot x_0\cdot w\cdot \cos \left(w\cdot t\right)\end{array}\right.$$

(18)

Solving the above Equation (18), the following expression can be obtained:

$$\left\{\begin{array}{l}m={\displaystyle \frac{k_1\cdot x_0+\frac{3}{4}\cdot k_3\cdot x_0^3-F_0\cdot \cos {\theta }_0}{x_0\cdot w^2}}\\ c={\displaystyle \frac{F_0\cdot \sin {\theta }_0}{x_0\cdot w}}\end{array}\right.$$

(19)

The mathematical expressions for the parameters m and c in the Bouc–Wen model can be derived from Equation (19) mentioned above.

$$\left\{\begin{array}{l}{k}_1=P_0\cdot A_0\cdot {\displaystyle \frac{\left(-r\right)}{h}}\\ k_2=P_0\cdot A_0\cdot {\displaystyle \frac{-r(-r-1)}{2h^2}}\\ k_3=P_0\cdot A_0\cdot {\displaystyle \frac{-r(-r-1)(-r-2)}{6h^3}}\end{array}\right.$$

(20)

When the cross-sectional area of the piston rod $A_s$, the cross-sectional area of the air spring chamber $A_0$, the initial volume of the air spring $V_0$, and the initial pressure $P_0$ are known, the values of $k_1$, and $k_3$ in Equation (20) can be obtained. When the frequency of the sinusoidal excitation $f={\displaystyle \frac{w}{2\pi }}$, the amplitude of the displacement $x_0$, the amplitude of the damping force $F_0$, and the phase difference ${\theta }_0$ are determined, the values of $m$ and $c$ are obtained. Equation (19) shows that both m and c are influenced by frequency and amplitude, with a clear frequency dependence.

Parameter $f_0$ represents the offset force of the MRD during testing. This is due to the inconsistency of the initial position. By analyzing the force–time curve of the MRD, the amplitude of the offset force $f_0$ can be determined based on the maximum and minimum damping force values of the image.

4.3. Identification of $\alpha ,\rho$

According to the analysis above, the model parameters can be obtained, except for $\alpha$ and $\rho$. Since there are only two unknown parameters, they can be identified by genetic algorithms. Genetic algorithm is an optimization algorithm that mimics natural evolution. First, an objective function is created to determine the deviation between the calculated and actual values. The results are then continuously corrected through crossover and mutation to search for solutions from the global domain, with the aim of minimizing the fitness function. The obtained solution can be considered optimal or the best within the given range [26]. The genetic algorithm’s basic flowchart is illustrated in Figure 8.

The algorithm requires basic settings prior to identification. This paper sets the initial population size of the genetic algorithm to 100, the crossover probability to 0.8, and the variance probability to 0.2, and the selection algorithm adopts random uniform distribution. The number of iterations is set to 500.

The fitness function is set as follows:

$$\mathrm{Fitness}={\displaystyle \sum _{i=1}^{m_0}({\widehat{F}}_t-F_t}{)}^2$$

(21)

where $m_0$ denotes the number of sampling points, ${\widehat{F}}_t$ is the damping force calculated according to the Bouc–Wen model, and $F_t$ is the actual damping force obtained. The specific values of $\alpha$ and $\rho$ can be obtained by genetic algorithm through multiple identification.

Before targeting parameter identification, it is necessary to define the parameter range. The upper and lower limits of the parameters are set to [0, 1 × 10⁸] to ensure a sufficiently large range. During the identification process, the range is gradually reduced until the optimal or better solution is obtained.

5. Establishment of Modified Bouc–Wen Model

Following the above process, all the unknown parameters in the Bouc–Wen model can be identified. At this point, Equation (22) can be rewritten as follows:

$$\left\{\begin{array}{l}{F}_d=P_0\cdot A_S+k_1\cdot x+k_2\cdot x^2+k_3\cdot x^3+{\displaystyle \frac{F_0\cdot \sin {\theta }_0}{x_0\cdot w}}\cdot \dot{x}+\\ {\displaystyle \frac{k_1\cdot x_0+\frac{3}{4}\cdot k_3\cdot x_0^3-F_0\cdot \cos {\theta }_0}{x_0\cdot w^2}}\cdot \ddot{x}+\alpha \cdot z+f_0\\ \dot{z}=\rho \cdot \left(\dot{x}-\left|\dot{x}\right|\cdot \left|z\right|\cdot z\right)\end{array}\right.$$

(22)

Table 1. Parameter identification results (20 mm).

	Current (A)	$\mathit{m} (\mathbf{kg})$	$\mathit{c} (\mathbf{N}\cdot \mathbf{s}/\mathbf{m})$	$\mathit{\alpha } (\mathbf{N})$	${\mathit{f}}_0 (\mathbf{N})$	${\mathit{F}}_0$	${\mathit{\theta }}_0$	$\mathit{\rho } ({\mathbf{m}}^{-1})$	$\mathit{\sigma }$	$\mathit{n}$
20 mm 1.0 Hz	0	134	12,322	2381	−271	1550	1.515	7550	1	2
	0.5	154	17,473	5653	−220.5	2198	1.515	1050	1	2
	1.0	167	20,724	8073	−140	2607	1.514	850	1	2
	1.5	191	24,281	9104	−90	3055	1.516	800	1	2
	2.0	195	25,315	9725	−265	3185	1.515	720	1	2
20 mm 1.5 Hz	0	34	11,061	2450	−267	2085	1.538	5650	1	2
	0.5	52	16,265	5888	−238	3066	1.539	850	1	2
	1.0	64.7	19,124	8573	−110	3605	1.535	755	1	2
	1.5	80.3	22,705	9539	−45	4280	1.537	680	1	2
	2.0	88.3	23,267	10,443	−10	4386	1.538	650	1	2
20 mm 2.0 Hz	0	12.3	10,839	2559	−269.5	2723	1.557	3090	1	2
	0.5	28	13,761	6532	−160	3457	1.558	750	1	2
	1.0	43	15,126	8750	−162	3800	1.558	700	1	2
	1.5	60.2	17,745	10,671	90	4458	1.557	650	1	2
	2.0	64.6	18,111	11,526	15	4550	1.557	630	1	2

and

Table 2. Parameter identification results (50 mm).

	Current (A)	$\mathit{m} (\mathbf{kg})$	$\mathit{c} (\mathbf{N}\cdot \mathbf{s}/\mathbf{m})$	$\mathit{\alpha } (\mathbf{N})$	${\mathit{f}}_0 (\mathbf{N})$	${\mathit{F}}_0$	${\mathit{\theta }}_0$	$\mathit{\rho } ({\mathbf{m}}^{-1})$	$\mathit{\sigma }$	$\mathit{n}$
50 mm 0.5 Hz	0	482	9912	828.1	−280	1672	1.429	15,080	1	2
	0.5	502.8	14,376	3968	−432	2280	1.429	2280	1	2
	1.0	518.5	17,989	6706	−459.5	2853	1.429	1180	1	2
	1.5	531.3	20,694	8405.6	−615	3282	1.428	870	1	2
	2.0	540.7	21,596	9351	−655	3425	1.429	820	1	2
50 mm 1.0 Hz	0	136.6	11,158	812.8	−347	3514.2	1.494	7575	1	2
	0.5	153	12,492	4395.2	−411	3934	1.495	1080	1	2
	1.0	169.3	13,825	7367.5	−500	4354	1.497	880	1	2
	1.5	189.3	15,458	9379	−555	4868	1.494	820	1	2
	2.0	191.6	15,642	9889	−695	4926	1.496	750	1	2
50 mm 1.5 Hz	0	32.3	12,201	1203.3	−329	5750	1.538	5780	1	2
	0.5	50.5	14,281	4350	−350.5	6730	1.539	980	1	2
	1.0	68.46	15,083	7588	−570	7108	1.537	880	1	2
	1.5	81.34	15,720	9676	−575	7408	1.539	790	1	2
	2.0	90.72	15,868	10,127	−850	7478	1.538	740	1	2

present the identification results for each parameter. Based on these results, a comparison between the identification and the test under different excitations is plotted. Figure 9 shows the comparison between the identification and the test results for different currents at an excitation amplitude of 20 mm and three frequencies (1.0 Hz, 1.5 Hz, and 2.0 Hz). Figure 10 shows the comparison between the identification and the test results for different input currents at an excitation amplitude of 50 mm and three frequencies (0.5 Hz, 1.0 Hz, and 1.5 Hz).

In the following, the variation of each parameter with current, excitation frequency, and amplitude will be analyzed to establish the functions.

5.1. Functions of c

The mathematical relationship between c and $F_0,{\theta }_0$ is established by Equation (19). Therefore, the relationship between $F_0,{\theta }_0$ and the current, the frequency, and the amplitude will be analyzed first. Figure 11 shows the variation curve of $F_0$ with current under different working conditions. It is observed that the trend of increasing $F_0$ flattens out with increasing current. This phenomenon can be described using quadratic functions. Meanwhile, altering the working conditions has no impact on the trend of $F_0$ with current. Therefore, the fitted curves merely shift up or down for different working conditions. This is because $F_0$ represents the damping force of the MRD, which comprises the viscous damping force and the Coulomb damping force. The impact of various working conditions on the MRD is primarily due to changes in velocity, which only affect the viscous damping force component. As a result, the effect of current on the damping output remains consistent despite changes in working conditions.

The function expressing the relationship between $F_0$ and current I for different working conditions obtained by fitting is given by Equation (23). Figure 12 and Figure 13 show the fitting curves and the comparison of R-square, respectively.

$$\left\{\begin{array}{l}{F}_0=-351.505\cdot I^2+1609.6\cdot I+1437\hspace{1em}20 \mathrm{mm} 1.0 \mathrm{Hz}\\ F_0=-351.505\cdot I^2+1609.6\cdot I+2402\hspace{1em}20 \mathrm{mm} 1.5 \mathrm{Hz}\\ F_0=-351.505\cdot I^2+1609.6\cdot I+2715\hspace{1em}20 \mathrm{mm} 2.0 \mathrm{Hz}\\ F_0=-351.505\cdot I^2+1609.6\cdot I+1600\hspace{1em}50 \mathrm{mm} 0.5 \mathrm{Hz}\\ F_0=-351.505\cdot I^2+1609.6\cdot I+3237\hspace{1em}50 \mathrm{mm} 1.0 \mathrm{Hz} \\ F_0=-351.505\cdot I^2+1609.6\cdot I+5812\hspace{1em}50 \mathrm{mm} 1.5 \mathrm{Hz}\end{array}\right.$$

(23)

where the unit of $F_0$ is Newton, and I denotes the current and its unit is ampere.

In the following, the effect of the change in working conditions on the parameter $F_0$ will be considered. Equation (24) shows the relationship between $F_0$ and current.

$$F_0=-351.505\cdot I^2+1609.6\cdot I+b$$

(24)

where b is the constant term of the quadratic function, which is affected by the frequency and amplitude of the excitation.

It is observed that there is a positive correlation between the constant term of this quadratic function and the product of excitation frequency and amplitude. A quadratic function was used for fitting. The fitted curve was obtained, as shown in Figure 14. The functional relationship is shown in Equation (25).

$$b=\mathrm{543,300}\cdot {(Af)}^2+\mathrm{25,764.2}\cdot (Af)+781.7095$$

(25)

where A denotes the amplitude of the excitation, and the unit is meter; f denotes the frequency of the excitation, and the unit is Herz.

It is observed that parameter ${\theta }_0$ remains constant under the same working conditions, indicating that ${\theta }_0$ is a function of frequency and amplitude. The identification results show that for both the 20 mm and 50 mm conditions, the value of ${\theta }_0$ falls between 1.4 and 1.5, while the value of $\sin {\theta }_0$ is greater than 0.99. Therefore, it is believed that any change in the value of ${\theta }_0$ resulting from a change in conditions has a negligible effect on c.

Equation (26) provides the expression for parameter c in terms of current, frequency, and amplitude in summary.

$$c={\displaystyle \frac{-351.505\cdot I^2+1609.6\cdot I+\mathrm{543,300}\cdot {(Af)}^2+\mathrm{25,764.2}\cdot (Af)+781.7095}{2\pi \cdot Af}}$$

(26)

5.2. Functions of m

The trend of m with current was fitted using a quadratic function. The function expressing the relationship between m and current I for different working conditions obtained by fitting is given by Equation (27). Figure 15 and Figure 16 show the fitting curves and the comparison of R-square, respectively.

$$\left\{\begin{array}{l}m=-6\cdot I^2+43.8\cdot I+133.4\hspace{1em}20 \mathrm{mm} 1.0 \mathrm{Hz}\\ m=-4.886\cdot I^2+37.15\cdot I+34.04\hspace{1em}20 \mathrm{mm} 1.5 \mathrm{Hz}\\ m=-5.829\cdot I^2+39.02\cdot I+11.35\hspace{1em}20 \mathrm{mm} 2.0 \mathrm{Hz}\\ m=-7.32\cdot I^2+43.84\cdot I+482.2\hspace{1em}50 \mathrm{mm} 0.5 \mathrm{Hz}\\ m=-7\cdot I^2+43.26\cdot I+135.2\hspace{1em}50 \mathrm{mm} 1.0 \mathrm{Hz} \\ m=-6.549\cdot I^2+42.67\cdot I+31.79\hspace{1em}50 \mathrm{mm} 1.5 \mathrm{Hz}\end{array}\right.$$

(27)

The relationship between the parameters m and current can be expressed by Equation (28). It is noted that the trends of the fitted curves for m and current are similar under different working conditions, with the working conditions only affecting the y-axis intercept of the curve. Based on Equation (27), it can be observed that changes in $a_m,b_m$ are negligible across different working conditions, while changes in working conditions primarily affect the amplitude of $d_m$. Therefore, the average value of $a_m,b_m$ is calculated for each working condition, and only the relationship between parameter $d_m$ and working conditions is analyzed. When the excitation amplitude is constant, $d_m$ decreases with increasing frequency. Simultaneously, altering the amplitude of the excitation has an insignificant impact on the parameters.

$$m=a_m\cdot I^2+b_m\cdot I+d_m$$

(28)

where $a_m,b_m,d_m$ are the constant terms.

Equation (29) is used to fit the relationship between $d_m$ and frequency f.

$$f(\lambda )=e_m\cdot {\lambda }^{h_m}$$

(29)

where $e_m,h_m$ are the constant terms, and $\lambda$ is the function independent variable.

Equation (30) is the fitting function of $d_m$ to frequency f. Figure 17 shows the curve of the fitting function.

$$d_m=114.8\cdot f^{-2.076}$$

(30)

Equation (31) provides the expression for m in summary.

$$m=-6.264\cdot I^2+41.62\cdot I+114.8\cdot f^{-2.076}$$

(31)

5.3. Functions of $\alpha$

The trend of $\alpha$ with current was fitted using a quadratic function. The function expressing the relationship between $\alpha$ and current I for different working conditions obtained by fitting is given by Equation (32). Figure 18 and Figure 19 show the fitting curves and the comparison of R-square, respectively.

$$\left\{\begin{array}{l}\alpha =-1829.4\cdot I^2+7904.6\cdot I+1927\hspace{1em}20 \mathrm{mm} 1.0 \mathrm{Hz}\\ \alpha =-1829.4\cdot I^2+7904.6\cdot I+2318\hspace{1em}20 \mathrm{mm} 1.5 \mathrm{Hz}\\ \alpha =-1829.4\cdot I^2+7904.6\cdot I+2867\hspace{1em}20 \mathrm{mm} 2.0 \mathrm{Hz}\\ \alpha =-1829.4\cdot I^2+7904.6\cdot I+811\hspace{1em}50 \mathrm{mm} 0.5 \mathrm{Hz}\\ \alpha =-1829.4\cdot I^2+7904.6\cdot I+1128\hspace{1em}50 \mathrm{mm}1.0 \mathrm{Hz}\\ \alpha =-1829.4\cdot I^2+7904.6\cdot I+1535\hspace{1em}50 \mathrm{mm} 1.5 \mathrm{Hz}\end{array}\right.$$

(32)

Based on the above analysis, Equation (33) shows the relationship between $\alpha$ and current I, considering the effect of working condition variation.

$$\alpha =-1829.4\cdot I^2+7904.6\cdot I+d_{\alpha }$$

(33)

where $d_{\alpha }$ is the constant term of the quadratic function, which is affected by the frequency and amplitude of the excitation.

It is observed that $d_{\alpha }$ increases with frequency for the same amplitude. When the amplitude is increased, $d_{\alpha }$ decreases for different frequencies. It is also found that the variation of $d_{\alpha }$ with frequency and amplitude is linear. A primary function was used to fit the relationship between $d_{\alpha }$ and frequency. The function obtained by fitting between $d_{\alpha }$ and frequency is presented in Equation (34). Figure 20 shows the curve of the fitting function.

$$\left\{\begin{array}{l}{d}_{\alpha }=837\cdot f+1130\hspace{1em}20 \mathrm{mm}\\ d_{\alpha }=837\cdot f+341.5\hspace{1em}50 \mathrm{mm}\end{array}\right.$$

(34)

By fitting the constant term in Equation (34) to the amplitude, the expression for the variation of $d_{\alpha }$ with frequency and amplitude at any working condition is shown in Equation (35).

$$d_{\alpha }=837\cdot f-\mathrm{26,280}\cdot A+1656$$

(35)

Equation (36) provides the expression for $\alpha$ in summary.

$$\alpha =-1829.4\cdot I^2+7904.6\cdot I+837\cdot f-\mathrm{26,280}\cdot A+1656$$

(36)

5.4. Functions of $\rho$

The relationship between $\rho$ and current is more affected by working conditions when the current is small. However, when the current is large, the relationship between $\rho$ and current remains consistent across different working conditions. This is because changes in working conditions primarily affect the mechanical characteristics of the damper through the speed of the damper piston. The damping force is affected by speed primarily through the viscous force, while the Coulomb force reflects the effect of current on the damping force. When the current is small, the damping force is primarily viscous force, making it more susceptible to changes in working conditions and having a greater impact on $\rho$. Conversely, when the current is large, the damping force is mainly Coulomb force, resulting in a smaller impact on $\rho$ with changes in working conditions.

The relationship between $\rho$ and the current is fitted. The function expressing the relationship between $\rho$ and current I for different working conditions obtained by fitting is given by Equation (37). Figure 21 and Figure 22 show the fitting curves and the comparison of R-square, respectively.

$$\left\{\begin{array}{l}\rho =689.6+6825.6\cdot {0.0098}^I\hspace{1em}20 \mathrm{mm} 1.0 \mathrm{Hz}\\ \rho =699.4+4951.9\cdot {0.0104}^I\hspace{1em}20 \mathrm{mm} 1.5 \mathrm{Hz}\\ \rho =709.4+2430.5\cdot {0.0112}^I\hspace{1em}20 \mathrm{mm} 2.0 \mathrm{Hz}\\ \rho =720.7+\mathrm{14,330.3}\cdot {0.0094}^I\hspace{1em}50 \mathrm{mm} 0.5 \mathrm{Hz}\\ \rho =698.2+6840\cdot {0.0102}^I\hspace{1em}50 \mathrm{mm} 1.0 \mathrm{Hz}\\ \rho =677.8+5059\cdot {0.0089}^I\hspace{1em}50 \mathrm{mm} 1.5 \mathrm{Hz}\end{array}\right.$$

(37)

The function used for fitting is expressed in Equation (38), where $b_{\rho }$ decreases as the frequency increases.

$$\rho =700+b_{\rho }\cdot {0.01}^I$$

(38)

Observing the fitted curves, it is found that the curves at the same frequency basically overlap, and the parameters of the fitted function are very similar. Therefore, it can be assumed that the parameters in the fitted function are only influenced by the frequency. The relationship between $b_{\rho }$ and frequency is fitted using the function presented in Equation (39), and the expression of the fitted function is obtained, as shown in Equation (40). Figure 23 displays the fitted curve of $b_{\rho }$ versus frequency.

Where $a_0,b_0$ are the constant terms, and $\lambda$ is the function independent variable.

$$b_{\rho }=6772\cdot f^{-1.087}$$

(39)

Equation (40) provides the expression for $\rho$ in summary.

$$\rho =700+6772\cdot f^{-1.087}\cdot {0.01}^I$$

(40)

5.5. The Modified Bouc–Wen Model

The modified Bouc–Wen model can be established by using the functional relationships of its parameters with respect to current, frequency, and amplitude, as shown in Equations (41)–(43).

$$\left\{\begin{array}{l}{F}_d=P_0\cdot A_S+k_1\cdot x+k_2\cdot x^2+k_3\cdot x^3+c\cdot \dot{x}+m\cdot \ddot{x}+\alpha \cdot z+f_0\\ \dot{z}=\rho \cdot \left(\dot{x}-\left|\dot{x}\right|\cdot \left|z\right|\cdot z\right)\end{array}\right.$$

(41)

$$\left\{\begin{array}{l}{k}_1=P_0\cdot A_0\cdot r\cdot h^r\\ k_2=P_0\cdot A_0\cdot {\displaystyle \frac{r(r+1)}{2}}\cdot h^r\\ k_3=P_0\cdot A_0\cdot {\displaystyle \frac{r(r+1)(r+2)}{6}}\cdot h^r\end{array}\right.$$

(42)

$$\left\{\begin{array}{l}c={\displaystyle \frac{-351.505\cdot I^2+1609.6\cdot I+\mathrm{543,300}\cdot {(Af)}^2+\mathrm{25,764.2}\cdot (Af)+781.7095}{2\pi \cdot Af}}\\ m=-6.264\cdot I^2+41.62\cdot I+114.8\cdot f^{-2.076}\\ \alpha =-1829.4\cdot I^2+7904.6\cdot I+837\cdot f-\mathrm{26,280}\cdot A+1656\\ \rho =700+6772\cdot f^{-1.087}\cdot {0.01}^I\end{array}\right.$$

(43)

6. Discussion

Section 4 presents the modified Bouc–Wen model. The aim of the following work is to predict the mechanical characteristics of the MRD under different working conditions. To achieve this, the sinusoidal excitation (80 mm, 1 Hz) is substituted into the modified Bouc–Wen model, resulting in the values of each parameter under this working condition, as shown in

Table 3. Parameter results in the modified Bouc–Wen model (80 mm 1 Hz).

$\mathit{I} (\mathbf{A})$	$\mathit{c} (\mathbf{N}\cdot \mathbf{s}/\mathbf{m})$	$\mathit{m} (\mathbf{k}\mathbf{g})$	$\mathit{\alpha } (\mathbf{N})$	$\mathit{\rho } ({\mathbf{m}}^{-1})$
0	12,573.17	114.8	390.6	8578.3
0.5	13,999.44	134.04	3885.6	1377.2
1.0	15,076.07	150.16	6465.8	767.7
1.5	15,803.04	163.14	8131.4	706.8
2.0	16,180.37	172.98	8882.2	700.7

. The mechanical characteristic curves of the MRD under different currents are obtained and compared with the experimental results, as shown in Figure 24. The comparison of the simulation identified by the traditional Bouc–Wen model and the test results are shown in Figure 25.

Meanwhile, the errors of simulation and test results at different currents are quantitatively analyzed. They present the average relative errors of predicted and test forces, as shown in Equation (44). The relative errors of the identified results at different currents are given in

Table 4. Comparison of the average relative error between the modified and the traditional Bouc–Wen model.

Current (A)	${\mathit{E}}_{\mathit{m}}$
	Modified Bouc–Wen Model	Original Bouc–Wen Model	Ali et al. [24]
0	3.87%	4.56%	/
0.5	2.82%	2.57%	/
1.0	2.45%	2.52%	/
1.5	2.19%	2.21%	3.13%
2.0	3.27%	2.96%	2.61%
3.0	/	/	2.35%
3.5	/	/	2.57%
4.0	/	/	2.63%
4.5	/	/	3.35%
5.0	/	/	3.75%
Mean $E_m$	2.92%	2.964%	2.913%

. The results are also compared with the average of the results reported [13].

$$E_m={\displaystyle \frac{{\displaystyle \sum _{i=1}^j\left|\frac{F_{pre}-F_{\exp }}{F_{\exp }}\right|}}{j}}\times 100\%$$

(44)

where $F_{\exp }$ denotes the test-measured damping force, $F_{pre}$ denotes the predicted force obtained by modified model, and $j$ is the number of test sites.

Figure 26 shows the accuracy and the identification rate comparison by using different identification methods at different currents for sinusoidal excitation at 1 Hz and 80 mm. The recognition method used by Ali [13] requires more parameters to be recognized as the recognition rate is slow, whereas the identification using the genetic algorithm directly is time-consuming and prone to falling into local optimal solutions.

Based on the comparison in Figure 26, it is obvious that the method adopted in this paper has a high identification rate and accuracy comparable to the original Bouc–Wen model. The accuracy of the modified model in this paper is close to Ali’s model but with fewer unknown parameters and a faster identification rate [13]. The proposed method is more efficient within acceptable error limits.

Table 5. Comparison of different kinds of predictive models.

Model	Parameters to Be Identified	Advantages	Drawbacks
This paper	2	High efficiency Acceptable accuracy	Requires experiments under different excitations
Dominguez [12]	16	High accuracy	Low efficiency
Ali [13]	10	High accuracy	Low efficiency Only the effect of frequency is considered
Boada [14]	7	High accuracy	Low efficiency Requires large quantities of test samples and long-time training
Lei [15]	/	No need to identify parameters	Exists some limitations at the pre-yield stage when applied to different amplitudes and frequencies
Lv [16]	/	No need to identify parameters	Exists some limitations in the descriptions for hysteresis

compares the modified Bouc–Wen model in this paper with other methods. The models proposed in [12,13] can accurately describe the characteristics of MRDs, but the increase in the parameters reduces the efficiency of the model. Although the model in [14] has fewer parameters and highly accurate model results, the neural network used requires a large number of samples. Refs. [15,16] propose the model by describing the rheological behavior of MR fluids. However, both models show some limitations.

7. Conclusions

MRDs are commonly used in various fields due to their excellent damping controllability. However, accurately and effectively modeling the MRDs has always been a challenge in research due to the obvious nonlinear and hysteresis characteristics. The Bouc–Wen model can describe a series of hysteresis curves by adjusting parameters. However, the traditional methods neglect the relationship between the model parameters and the frequency and amplitude and can only be adapted to a single working condition. This paper proposes a modified MRD Bouc–Wen model. The model is designed to be adaptable to different excitation conditions. The physical significance of each parameter in the Bouc–Wen model is analyzed, and the mathematical relationship between the parameters and the frequency and amplitude is obtained using the harmonic balance method. Additionally, the modified model is proposed based on the functions between the parameters and the current, frequency, and amplitude. The accuracy of the model is verified by comparing the simulation with test results under a different condition. The MRD in the paper is used to improve vehicle driving comfort. Hence, the model has been analyzed under 0–2 Hz excitation (the natural frequency of the body model). This paper presents research on efficiently and accurately predicting the mechanistic characteristics of MRDs under different excitation conditions. The new model proposed in this paper can be applied to different frequencies and amplitudes of excitation. It can predict the dynamics of MRDs under different excitations. The error of the predictions is within the acceptable range, while the prediction model is highly efficient with fewer identification parameters. The modified model can be used for the research of damping force control strategy. The inputs of the modified model contain the frequency and amplitude of the excitation, which can be obtained by real-time frequency domain analysis through FFT. Therefore, the relevant control strategy can be established to achieve the damping force control under different working conditions. The main findings and conclusions are outlined below:

(1)

The mechanical properties of the MRD were tested, and the causes of the variations of these properties were analyzed under different currents, frequencies, and amplitudes. The results indicate that the damping output of the MRD is affected not only by the current but also by changes in the excitation frequency and amplitude.

(2)

The Bouc–Wen model was used for the preliminary modeling of the MRD, and a method of identifying the model parameters based on the harmonic balance method was proposed. Firstly, the unknown parameters in the model were analyzed for their specific physical meanings. Different methods were then adopted for each parameter. $\sigma ,n$ were taken as 1 and 2, respectively, for their low sensitivity. The mathematical relationships between m, c with the frequency and amplitude were obtained by the harmonic balance method. $\alpha ,\rho$ were directly identified by genetic algorithm. The simulation results were compared to the original data and were found to be in good agreement. The model parameters consider the coupling relationship between the excitation frequency and amplitude, resulting in a smaller number of parameters to be identified. The improved method greatly improves computational efficiency compared to previous identification methods.

(3)

The model parameters were established as the functions of current, excitation frequency, and amplitude, and a modified nonlinear model of the magnetorheological damper was presented. The agreement of the fitted curves was compared by the R-square under different excitation conditions. The results show that the established functional relationship matches the identification results well.

(4)

The simulations of the modified Bouc–Wen model for the MRD were obtained and validated against experiments under the excitation condition of 80 mm and 1 Hz. The simulation results’ average relative errors were calculated and compared with those of the conventional Bouc–Wen model, which was built using the genetic algorithm and other modeling approaches in the literature. The results indicate that the modified Bouc–Wen model proposed for MRDs can accurately and efficiently describe the mechanical properties of the dampers under various working conditions.