Nonlinear characterization of magnetorheological elastomer-based smart device for structural seismic mitigation

ABSTRACT Magnetorheological elastomer (MRE) has been demonstrated to be effective in structural vibration control because of controllable stiffness and damping properties with the effect of external magnetic fields. To achieve a high performance of MRE device-based vibration control, a robust and accurate model is necessary to describe nonlinear dynamics of MRE device. This article aims at realising this target via nonlinear modeling of an innovative MRE device, i.e. MRE vibration isolator. First, the field-dependent properties of MRE isolator were analysed based on experimental data of the isolator in various dynamic tests. Then, a phenomenal model was developed to account for these unique characteristics of MRE-based device. Moreover, an improved PSO algorithm was designed to estimate model parameters. Based on identification results, a generalised model was proposed to clarify the field-dependent properties of the isolator due to varied currents, which was then validated by random and earthquake-excited test data. Based on the proposed model, a frequency control strategy was designed for semi-active control of MRE devices-incorporated smart structure for vibration suppression. Finally, using a three-storey frame model and four benchmark earthquakes, a numerical study was conducted to validate the performance of control strategy based on the generalised current-dependent model with satisfactory results. Graphical abstract


Introduction
As a type of novel smart material, the magnetorheological elastomer (MRE) is composed of magnetic particles suspended in the elastomeric matrix [1][2][3]. Generally, two types of MRE, i.e. isotropic and anisotropic, are manufactured either without the magnetic fields or under the fields. When the external magnetic fields are applied to such materials, the material's mechanical properties such as shear modulus, stiffness, and damping will be varied accordingly [4,5]. Due to the unique characteristic of field-dependence, the MRE has been widely utilized in the design of smart devices for various engineering applications, such as vibration absorbers, vibration isolators, sensors, and actuators [6][7][8][9][10].
A large number of studies have been reported on the MRE-based devices. Li et al. proposed the first seismic isolator based on laminated MRE and steel sheets, which is able to provide the loading capacity of 375 kg in the vertical direction with maximum lateral displacement and no applied magnetic fields [11]. Under the effect of the magnetic field, the seismic isolator can contribute the increases of 45% on the maximum shear force and 38% on the lateral stiffness with 3 A applied current. Yang et al. developed another MREbased isolator using a similar multi-layered configuration [12]. Different from the device proposed by Li et al. with the stiffness hardening, this isolator has the feature of stiffness softening when supplied with external magnetic fields, which makes the device operate in the passive mode without any power. Ni et al. developed an adaptively tuned vibration absorber (ATVA) based on MRE with both squeeze and shear modes, in which the compressive loads to the MREs can be tuned by a mobile piezoelectric actuator [13]. Kallio designed a MRE-based adjustable spring component, which works in the squeeze mode [14]. This innovative device has the capacity against the compressive loads applied on it via tuning the magnetic fields, therefore presenting adjustable characteristics. Bose et al. utilized soft MRE to design a valve-type device to adjust the airflow [15]. The main principle of this device is that the gap between the MRE ring and the yoke contracts with the addition of an external magnetic field; therefore, the airflow that passes through the gap can be tuned in real-time via the field changing.
To well utilize the MRE-based devices in practice, it is absolutely necessary to model and simulate the dynamic behavior of the MRE under different excitation conditions. To date, there are two categories of models for MRE or MRE-based devices. One type is micromodel, which is on basis of law of conservation in continuum mechanics, micromechanism and mutual effect of magneto-dipole in particle-chain. The other type is macro-model, which is developed on basis of macroscopical elasticity property of the MRE, fielddependence feature of smart elastomer, and nonlinear relationship between force and deformation of the MRE. Research in modeling MRE or MRE-based device was started by Davis in 1999, focusing on characterizing magneto-mechanical property of the MRE without magnetic field and at saturation [16]. Li et al. proposed a linear phenomenal model with four parameters to predict viscoelasticity property of the MRE [17]. In this model, a spring component is added in parallel with a classical three-parameter solid model to represent the field-dependence of the shear modulus. To study viscoelastic behavior, field-induced property, and interface slip between particles and matrix of the MRE, Chen and Jerrams developed a rheological model consisting of a stiffness variable spring, a standard linear solid model, and a spring-Coulomb friction slider [18]. This model is able to forecast the dynamic response of the MRE with any matrix property and iron particle fraction through the parameter variation. Eem et al. designed a dynamic model to investigate the MRE's behavior under shear deformation, which is made up of the Maxwell model and the Ramberg-Osgood model [19]. In this model, the coefficients are the functions of the magnetic flux density. Norouzi et al. put forward a novel mathematical model based on revised Kelvin Voigt viscoelastic model to portray nonlinear relationships between shear stress and shear strain of the isotropic MRE [20]. Different from other models, this model has constant coefficients at different loading conditions and magnetic flux densities. Compared with the research on MRE models, the modeling studies on MREbased devices are rarely reported. To portray the dynamic responses of MRE regrading broad-band loading amplitude, frequency, and current level, Nguyen et al. developed a nonlinear hybrid model comprised of a magnetic interaction model, fractional Maxwell model, and adaptively smooth Coulomb friction model, in which magnetic dipole model was used for neighboring particles to calculate magnetic force with respect to particle volume, and fractional viscoelasticity model and adaptively smooth Coulomb friction were employed to characterize dynamic responses and current-dependent features of MRE in the broad ranges of loading frequency and amplitude [21]. Sun et al. proposed a magnetic dipole theory-based model to characterize the frequency shift of an MRE absorber working in squeeze mode [22]. Li and Li presented a phenomenal model to predict the shear force of an adaptive MRE isolator, in which a damping element and a power law function are added into standard three-parameter solid model to portray the strain stiffening phenomenon of the device [23]. They also developed another analytical model based on Koh Kelly model for the same device, which includes an elastic spring, a strain stiffening spring and a Coulomb damping dashpot [24]. A dynamic model, consisting of a smooth Coulomb component, a stiffness adjustable spring component and a three-parameter solid model, was also presented for evaluating the procedure of MRE device operation for further isolation system development in the practical application [25]. Besides, soft computing techniques were applied to the nonlinear behavior prediction of the MRE-based devices as well. Yu et al. developed a support vector machine (SVM)-based model, in which the device displacement and velocity as well as externally applied current are used as inputs and the model output is the shear force [26]. The model parameters are optimized via a pre-training procedure. Similarly, Fu et al. employed the nonlinear autoregressive exogenous (NARX) neural network with the three-layer architecture to estimate nonlinear dynamics of MRE isolator in shear-compression mixed mode [27]. This model also considers the historical displacement, velocity, current, and force as the inputs, therefore being more optimal and robust than other network-based models. However, the main disadvantage of this type of MRE model is that the model parameters cannot directly be related to the physical phenomenon of the device. If it is adopted for the controller design in practice, the model which is capable of connecting the model parameters with the nonlinear behavior of the device is definitely necessary.
On the other hand, a great number of control algorithms have been proposed and applied to semi-active control of MRE devices for suppressing structural vibrations in realtime. Brancati et al. developed a hybrid method consisting of regression neural networks and model-predictive control for vibration suppression of lightweight structures [28]. Guo et al. put forward the optimized proportional-integral-derivative (PID) control algorithm for the vibration mitigation of precision platforms incorporated with MRE isolators, in which the values of controller are optimally initialized by PSO and then updated by fuzzy algorithm [29]. Chen et al. designed an adaptive dual-loop controller for solving the time delay problem of MRE isolators [30]. In this work, a numerical investigation was conducted to illustrate the superiority of proposed controller via a comparison with Lyapunov-based controller on a seismically loaded five-storey benchmark building embedded with MRE isolators. Nguyen et al. proposed a robust control scheme, composed of three controllers of adaptive sliding mode, standard linearizing and single robust, for MRE isolator with robust stability and complete convergence of displacement responses, the performance of which was numerically validated by a two-storey building model subjected to earthquake loads [31]. Yang et al. utilized fuzzy logic (FL) controller to experimentally mitigate vibration of a scaled MRE devices-incorporated smart shear building excited by El-Centro earthquake [32]. Similar work can also be found in [33] on using FL to design semi-active controller for MRE-based variable stiffness isolation system. To improve the control performance, Fu et al. introduced genetic algorithm (GA) to design nonlinear selftuning law for FL controller to attenuate time-variant vibration in low frequency, the capacity of which has been validated by both numerical and experimental investigations [34]. Nevertheless, priori information is needed to design this type of controllers, and prompt tuning of controller's parameters is not easy to realize in the real application.
In conclusion, to make full use of MRE devices for structural vibration suppression, a generalized model, with a comprehensive capacity to forecast dynamic behavior under various excitation conditions, should be developed considering the field-dependent characteristic. Accordingly, this article attempts to achieve this target via designing an accurate and robust phenomenal model that is highly adaptable to changeable external loading. First, the experimental data of an MRE isolator are collected via the dynamic tests with various loading conditions. Following that, an analytical study is conducted to investigate the device output dependencies on applied current level, excitation frequency and amplitude. Then, a generalized field-dependent model based on the Bouc-Wen component, which is versatile in describing various characteristics of hysteresis behavior [35], is proposed to simulate nonlinear response of the device under external loads. During this process, a Kevin-Viogt element is adopted to illustrate the viscoelastic behavior of the device and the Bouc-Wen hysteresis element is used to describe the strainstiffening phenomenon in the response. The model parameters are also analyzed to find their correlations with the model outputs. Based on the collected experimental data with fixed excitation frequency and applied current, the model parameters are identified by improved particle swarm optimization (IPSO) via minimizing the deviations between the real values and model outputs. Moreover, the identification results with different loading conditions are summarized to obtain a generalized field-dependent model. Finally, a frequency control strategy based on generalized model is developed, and a numerical investigation on three-storey shear frame embedded with MRE devices as an isolation layer subjected to seismic loads is conducted to validate its performance with satisfactory results. The primary contributions of this study include: (1) developed a solid and highly accurate phenomenological model to characterize nonlinear and hysteretic responses of MRE device under the effect of magnetic field; (2) proposed an IPSO for estimating model parameters with high identification accuracy; and (3) designed a frequency control strategy based on the proposed model for structural seismic suppression. The outcomes of this research offer a potential solution to the infrastructure protection against external hazard loads in the actual situation.

Proposed mathematical model
The MRE devices have controllable stiffness and damping properties, with nonlinear and hysteretic responses, under the effect of external magnetic fields. To optimize the performance of MRE device in structural vibration mitigation, a solid and highly accurate model should be established to depict these characteristics. The main challenge of characterizing nonlinear behavior of the MRE device is to accurately acquire the strain stiffening and softening phenomena of the MRE due to increasing magnetic fields. In this study, a novel phenomenal model is proposed and its physical configuration is shown in Figure 1. In this model, a Kevin-Voigt component, consisting of a linear spring element and a viscous dashpot element, is used to depict the viscoelastic behavior of the MRE isolator. Its strain stiffening behavior is captured by a Bouc-Wen element, which is commonly used in the material and structural fields to indicate various hysteresis behaviors due to simple mathematical expression and strong capacity. As can be seen in the figure, three elements are connected in parallel in the model. Consequently, the mathematical expression of the proposed model is given in Equation (1): where c 0 and k 0 denote the damping and stiffness coefficients, respectively. F 0 is the initial force of the device, which can be regarded as the offset. α is a scaling factor and y denotes the intermediate variable, given by In Equation (2), parameters A, n, γ, and β are related to the range and shape of the hysteresis loop. The value of A has a direct impact on the maximum boundary of the model. n denotes the factor that controls the smoothness of the hysteresis loop at the inflection point. γ and β are mainly responsible for the loop shape and different combinations of γ + β can generate different hysteresis shapes.

Model parameter analysis
To better employ the proposed Bouc-Wen-based model to predict the dynamic response of the MRE isolator caused by external loadings, a parametric analysis should be conducted to investigate the effects of the model parameters on the model output responses. In this study, take a single degree-of-freedom (DOF) MRE isolator system as the example, and the reference values of the model parameters are set as follows: k 0 = 20, c 0 = 0.3, α = 20, A = 3, β = 3, γ = −2, and n = 1. Besides, an external harmonic loading is used to drive the system, with the following expression: where Amp denotes the loading amplitude and f denotes the excitation frequency. In this simulation, Amp = 4 mm and f = 2 Hz. The sampling frequency is set as 200 Hz, and 50 points (one cycle) are collected for each response loop. Figure 2 shows the resultant hysteresis loops with different values of k 0 , when other model parameters are set as reference values. The value of k 0 changes from 10 to 40 with an interval of 10. It is seen from the results that the maximum force and effective stiffness (slope of curve) change almost linearly with k 0 . Another phenomenon worth noting is that all the loops intersect at two points, as pointed out in Figure 2(a). Besides, the forcevelocity loops join together at two points as well, as shown in Figure 2  this unique behavior tends to be more evident from these points. Furthermore, in the force-velocity responses, the enclosed area of each curve enlarges as the value of k 0 increases. The effect of the value of parameter c 0 on the model output is described in Figure  3. Unlike parameters k 0 , c 0 has little influence on the maximum shear force and effective stiffness of the system, which indicates that all the generated forces under different values of c 0 reach the maximums with the same value at two endpoints. In the meantime, the force-velocity loop represents that peak shear force occurs when the velocity is zero and decreases to minimum value at the maximum velocity, which indicates that the damping force has little contribution to the total shear force. Such a conclusion can be supported by the fact that the value of c 0 is generally used as a small value. It can also be observed that the variation in c 0 can result in slight differences among corresponding hysteresis loops, not to mention the maximum shear force. Nevertheless, it can still be noticeable from Figure 3 that the forcedisplacement loop grows wider with the increasing value of c 0 , which can cause a larger energy dissipation capacity. Figures 5 demonstrate the effects of parameters α and A, respectively, on the model dynamic output. An obvious similarity between both figures can be observed, which means that any one of parameters can be replaced by the other one if the hidden relationship between two parameters can be found. It will not be discussed in this work, but the Bouc-Wen model can be improved from this point. As aforementioned, α is post-to-pre-yielding stiffness ratio, which represents the linearity level of the hysteresis loops. As can be seen in    to be more obvious as the value of β gets smaller. On the contrary, when β grows, the hysteresis shapes are inclined to be linear ellipses. It is noteworthy that an effective β should be kept positive. In other words, β should fluctuate in an effective range for reproducing reasonable hysteresis shapes on the condition that all the other referenced parameters are fixed. As for the influences of parameter γ, the variation trend of hysteresis loops subjected to the tuning of parameter value is quite similar to that of β. In the case when γ is negative, the nonlinearity level increases with the increasing absolute value of γ. Otherwise, the hysteresis loops progress linearly. It is also observed that even though the nonlinearity of the model reduces with the increase of both values of β and γ, the influence of parameter value's variation on the hysteresis shape becomes rather limited when the parameter reaches a certain value. Figure 8 shows the comparison of model outputs with different values of n. Obviously, n is a parameter to just tune the smoothness of the hysteresis loop. Its value has little effect on the loop shape and maximum shear force.

Model parameter identification
Parameter identification of the proposed model is a global optimization procedure, in which the optimization objective is to minimize the errors between the experimental data and model outputs. In the proposed model, the offset F 0 can be calculated by averaging the minimum and maximum values of the force in one loading cycle. Hence, k 0 , c 0 , α, A, β, γ and n are parameters to be identified. However, the parameters in the model for the MRE-based device are field-dependent, which means that these parameter values may be different under different loading conditions. In this case, for each condition, the model parameters should be identified based on corresponding experimental data. Then, the correlations of parameter values with the loading frequency, amplitude and current level are explored to acquire a generalized field-dependent model that is appropriate for changeable loading conditions. Besides, the main challenge for parameter identification is to solve the nonlinear differential equation (Equation (2)) in the model expression. Here, 4-order Runge-Kutta algorithm is employed to calculate the evolutionary variable y in an iteration way. The mathematical expressions are given in Equations. (4)-(8) as follows: where ∆t denotes the time-slot of data sampling in the dynamic test and its value is 0.00390625 s (1/256) in this study. After the values of y are obtained at each time-slot, the parameter identification of the proposed model can be deemed as solving the following optimization problem.
min Obj ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi s:t:all the parameters are positive: In Equation (9), Obj denotes the objective function, which is the root mean square error (RMSE) between real values and model outputs in one loading cycle. N denotes the total number of data that are used for parameter identification. Obviously, this objective function is highly complicated and difficult to be solved because of a large number of identified parameters, especially for the exponential parameter n. For such a problem, the common direct research methods may not be effective due to the fact that the gradient information is difficulty represented. Therefore, an improved particle swarm optimization (IPSO) algorithm is proposed in this work for the parameter identification of the MRE device model.
In the standard PSO, each particle generally updates its velocity and location information according to itself optimum Pb and swarm optimum Gb, which indicates that all the particles will move toward the optimal location. Previous studies have shown that the particles in the swarm will gather around the optimal location Gb finally regardless of the algorithm premature convergence or global convergence [36,37]. Nevertheless, if Gb is the local optimum rather than global optimum, the algorithm will fall into this local solution and will be difficult to jump out the constraint of Gb to search for global optimal solution. To fix this problem, the proposed IPSO tries to improve the convergence rate and accuracy via effective control of Pb and Gb. As is well known, in addition to summarizing its experience, the individual observes the behaviors of other individuals to enhance itself in the biological evolution procedure. Hence, a new mechanism is developed to update the Pb and Gb during the algorithm iteration. First, for ith particle, randomly select individual optimum from other particles in the swarm, denoted by Pb j . The new individual optimum of ith particle can be expressed as follows: where λ is a random number between 0 and 1. Then, sort and select K best individuals (Pb � 1 , Pb � 2 , . . ., Pb � K ) and calculate mean value as new global optimum Gb* using the following equation: The primary objective of two modifications is to improve the diversity of particle swarm and to reduce the possibility that the algorithm traps into the local optimum. Hence, even in the later stage of evolutionary iteration, the algorithm still can conduct fine-tuning of optimal particle location to polish the solution accuracy. This is also the innovation of proposed IPSO. Based on it, the identification procedure of the model parameters can be summarized as the following steps.
Step 1. Determine the objective function and parameters to be identified. In this study, the objective function is given in Equation (9) and identified parameters include k 0 , c 0 , α, A, β, γ and n.
Step 2. Confirm the parameter values of IPSO, including swarm size N p , maximum iteration number N m , inertia weight w, learning factors c 1 and c 2 , and K.
Step 3. Initialize the velocity and location of each particle in the swarm using the following equations: where v i and x i denote the velocity and location of ith particle, respectively; v max denotes the maximum value of the velocity; x l and x u denote the lower and upper boundaries of the location of ith particle. Here, v max = 2, x l = [0 0 0 0 0-2 1] and x u = [3, 50 5,5 5,10].
Step 4. Calculate the value of the objective function based on current particle location and record the individual optimum Pb of each particle and swarm optimum Gb.
Step 5. Compare the current iteration number n c with maximum iteration number N m . If n c >N m , go to Step 7. Otherwise, go to Step 6.
Step 6. For all the particles in the swarm, conduct the following operations: Step 6.1. Update the velocity and location of ith particle using Equations. (14) and (15).
Step 6.2. Calculate the value of the objective function based on updated location, and record Pb and Gb through the comparison between current and previous values.
Step 6.3. Update the Pb i using Equation (10).
Step 6.4. Update the Gb using Equation (11). If Gb* < Gb, replace Gb with Gb*. Otherwise, keep current Gb unchanged.
Step 7. Output optimal Gb and x (parameters), and terminate the algorithm. The procedure of the IPSO to identify the parameters of the proposed Bouc-Wen model is also shown in Figure 9.

Setup and data collection
To obtain sufficient data for the task of interest in this research, an MRE isolator prototype was developed and tested dynamically using a shake table. The MRE isolator employs a design of laminated architecture, in which 26 thin steel layers and 25 MRE sheets are alternatively placed with the diameter of 120 mm and the thickness of 1 mm. All the MRE and steel layers are arranged into the cylindrical shape with a 120 mm diameter. The top and bottom of laminated structure are rigidized by two steel plates, which make up the cole of the device. An electromagnetic coil is placed outside the core to produce the magnetic fields across smart materials. A nonmagnetic support, made of epoxy, is used to contain both the core and coil. The MRE material utilized in the device development is comprised of silicone oil, rubber and carbonyl iron particles, to accomplish better MR effect. A complete description of MRE material fabrication and device design could be found in [11]. The working principle of MRE isolators can be summarized as follows. They are generally deployed underneath the protected structures to form a base isolation layer. With the seismic loads, the stiffness of MRE isolator can be adjusted in real-time, so the smart system parameter can be altered, changing the fundamental frequency of the structure to an ideal value to avoid the resonance accordingly. Hence, the dynamic tests of MRE isolator are conducted, the responses of which should be further investigated for better application in structural seismic mitigation. The experimental setup is shown in Figure 10. The shake table was able to provide horizontal loadings to the device in the dynamic mode. The bottom plate of the isolator was fixed on the table and moved along with the motion of the table while the top plate of the isolator is fixed to a reference frame so as to make sure that the displacement of the shake table equals the horizontal deformation of the device. In the meantime, a load cell was amounted to the fixed reference frame to measure the horizontal load from the device. The key to this setup rests with that the top plate of the isolator and load cell remain relatively static during the test, therefore avoiding undesirable inertia force in the measurements. A DC power supply was adopted to provide the DC current to energize the isolator and a slider was utilized to tune the current level applied to the device, which are shown in Figure 10.
Three types of excitations are considered in the dynamic tests: harmonic, random and seismic excitations. For the harmonic excitations, the amplitudes are 2 mm, 4 mm, and 8 mm, the frequencies are 1 Hz, 2 Hz, and 4 Hz, and the applied current levels are 0 A, 1 A , 2 A, and 3 A, which corresponds to different magnetic fields. For the random excitation, the maximum amplitude is 4 mm, and the frequency varies from 1 Hz to 10 Hz. For the seismic excitation, the scaled El-Centro earthquake waves are considered and the maximum Update velocity and location of particle using Eqs. (14) and (15) Calculate objective function based on updated location and update Pb and Gb Update Pb using Eq. (10) Update Gb using Eq. (11) n c = n c + 1 Output optimal Gb and parameter values End Input experimental data Yes No Figure 9. Flowchart of IPSO to identify parameters of the MRE isolator model.
amplitude is 4 mm. The applied current randomly varies between 0 A and 3 A for both random and seismic excitations. For all the tests, the sampling rate is 256 Hz and sampling time is 60 s for obtaining the stable results. To avoid the performance affected by the rising temperature due to changeable current, a temperature adjustment mechanism is utilized to keep the temperature nearly stable in each test. The force and displacement data of the device can be obtained directly from the sensor readings, and the velocity data can be obtained via the differential calculation from time-displacement responses.

Data analysis
In this part, the experimental results from the harmonic excitations are used to analyze the dynamic behavior of the MRE isolator. First, we try to investigate the effect of the magnetic field on the device response. Figure 11 plots the hysteresis loops of the MRE isolator under different applied currents (magnetic fields), when excitation frequency is 2 Hz. In this figure, the slope of the loop denotes the equivalent stiffness of the isolator. Therefore, it is clearly seen that the device stiffness is dependent on magnetic field strength and ascends with the increase of the current level. This phenomenon is deemed as stiffness hardening. Via the response comparison of three figures, it can be also found that the increase of stiffness is more prominent in small  Figure 10. Experimental setup and testing equipment. loading amplitude. Besides, the elliptical loop area represents energy dissipation in each cycle, which corresponds to the damping property of the device. With the addition of a magnetic field, the equivalent damping is increased as well. Accordingly, it can be concluded that the stiffness and damping properties of the MRE isolator are the functions of the applied current level. Figure 12 shows the force responses of the device under different excitation amplitudes, when the loading frequency and applied current are 2 Hz and 1 A, respectively. It can be observed that different from stiffness hardening phenomenon, the slope of response loop is slightly declined with the increase of loading amplitude, the effect of which is stiffness softening. This phenomenon is pretty obvious when the device is under high applied current input. On the other hand, when the loading frequency and applied current are fixed, the elliptical loop area will be enlarged with increase of loading amplitude, which indicates that the increasing amplitude of the excitation will contribute to more energy dissipation by the isolator in each cycle. Figure 13 gives the force responses of the MRE isolator under different loading frequencies, when the excitation amplitude and applied current are 4 mm and 1 A. Apparently, different from applied current and loading amplitude, the loading frequency has little effect on the equivalent stiffness, energy dissipation and maximum shear force of the device.

Modeling results and discussions
In this section, the experimental data captured from the harmonic loads are used to set up the generalized field-dependent model for the MRE device. Then, the data from random and seismic excitations are used to validate the performance of the proposed model, because in the practice the amplitudes and frequencies of the external excitations are always changeable. Hence, it is necessary to evaluate the performance of the proposed field-dependent model to trace the force variation caused by changeable loading Simulink, which is shown in Figure 14. It is clearly seen that in the model there are a total of seven parameters to be identified, denoted by x(1) to x (7). Besides, the model inputs are displacement and velocity of the device in this identification procedure, since the maximum amplitude, frequency, and current are always fixed at each harmonic loading condition. Before the model parameter identification, the values of IPSO parameters should be set appropriately. Among parameters in IPSO, the inertia weight w is an important factor that can directly affect the convergence rate as well as identification accuracy. If w is assigned with a small value, the algorithm will converge very slowly. On the contrary, if w is assigned with a big value, the algorithm will have a fast convergence rate but may fall into the local optimum. To avoid this problem, this work adopts a self-adaptive inertia weight, the value of which is adaptively adjusted during the procedure of algorithm iteration. Generally, in the early stage of iteration, a big value of inertia weight is selected to accelerate the convergence, while a small value of inertia weight is used in the late stage of iteration to improve the identification accuracy. Hence, we should design the inertia weight as a decreasing way with the addition of iteration number. Two types of decreasing function can be considered for possible solutions: linear and nonlinear decreasing solutions, the compassion of which is shown in Figure 15. It is apparent that the nonlinear decreasing function is better than linear one because it can keep a big inertia weight for a longer time which is beneficial to global search and will rapidly decrease to the small value which is beneficial to local search. Accordingly, the nonlinear self-adaptive inertia weight will be employed in this case with the following expression:

Responses of 4mm amplitude and 1A current
where w min and w max denote the minimum and maximum inertia weights, respectively. In this work, the values of w min and w max are set as 0.4 and 0.9. Besides, the population size N p and maximum iteration number N m are another two significant parameters, which are related to final identification accuracy. In this work, the trial-and-error method is employed to select an optimal combination of N p and N m via comparing the solutions of the objective function. The parameter range settings are given as follows: N p ranges between 10 and 60 and N m is in the range of [200,800]. Furthermore, other parameters of IPSO are set as: K = 5 and c 1 = c 2 = 1.5. Figure 16 compares the RMSE values of the model identified by IPSO with different combinations of N p and N m . It is noticeable that with the increase of N p and N m , the value of RMSE decreases gradually. However, after they arrive at certain values, the change of RMSE will not remarkable and RMSE value is even kept at a stable level. Moreover, it can be observed that the RMSE is more sensitive to the population size than maximum iteration number. When the population size is 60, the RMSE shows little change even though the N m is increased from 200 to 800. From the figure, the optimal combination of N p and N m is (50, 500) due to the smallest RMSE error. Based on optimal algorithm parameter setting, the parameters of the proposed MRE device model are identified at each harmonic loading condition. Figures 17-19 illustrate the comparison between experimental results (straight line) and model prediction results   (17)- (19), where F BW and F real denote the predicted and real forces, respectively, and N is the data number of one loading cycle. In this study, the numbers of N are 256, 128 and 64 for the cases of excitation frequencies of 1 Hz, 2 Hz,   Table 1 displays the evaluation indices of model performance under different loading scenarios. It can be observed that with the increments of loading amplitude and applied current, the errors between predicted and real forces become larger, which results in increasing RMSE and MAE. However, regardless of excitation conditions, the values of R 2 always keep at a stable level, which is above 0.98, further demonstrating the high performance of the proposed model.

RMSE ¼
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 N Then, based on the identification results, a parametric study should be conducted to find the relationships between model parameters and excitation conditions, such as loading amplitudes, frequencies and current levels. Nevertheless, since the information of excitation frequencies and amplitudes are always unknown and unpredicted in practice, there is no need to investigate their dependences on model parameter values. The generalized model with the parameters of current-dependency property is imperative to realize the adjustment of the shear force generated by the device against external loading. In this part, the results of model parameters identified from different loading conditions are grouped in accordance with applied current levels. In each level, the parameters' values identified from different loading frequencies and amplitudes are averaged. Then, the mathematical expressions are developed to describe the relationships between parameters' averaged values and current levels via the way of curve fitting, which is shown in Figure 20. It is clearly seen that except α, all the model parameters have almost linear relationships with applied current level, which indicates that the first-order polynomial functions can be used to express these relationships. For parameter α, a second-order polynomial function is suitable to portray its current-dependent property. The proposed IPSO is used to estimate optimal values of the coefficients in the polynomial functions, the results of which are given in Equations (17)   Based on the field-dependent parameters, a generalized Bouc-Wen-based model can be obtained and its mathematical expression is given as follows: The effectiveness of this generalized field-dependent model is evaluated using the data captured from the MRE isolator under a random excitation with maximum amplitude of 4 mm and a magnitude scaled El-Centro earthquake with randomly variable current inputs between 0 A and 3 A. Figures 21 and Figures 22 show time-historical comparisons between experimental data and model predictions, when the device is loaded with random excitation and scaled earthquake, respectively. It is apparent that the predicted force responses outstandingly agree with the experimental data for the random case. Even though there are obvious errors in several areas for the scaled earthquake case, the model is capable of well tracing the changing trend of the force with the external loading, which is also acceptable for the practical application. According to the fitting results in Figures 21 and Figures 22, the proposed generalized field-dependent model has proved its capacity to dynamically predict the shear force of the MRE isolator under instantaneously changeable loading conditions, which are similar to strong winds or earthquakes in real situation. Figure 23 quantitatively appraises the capability of the generalized model using three evaluation indices RMSE, MAE and R 2 , in terms of radar plot. Although scale El-Centro earthquake has smaller RMSE and MAE, random excitation case has higher R 2 . Overall, the model performance is acceptable, since the R 2 is above 0.9, satisfying the requirement of the modeling study. The outcomes of this result can be perceived as a good solution for the design of adaptive controllers for structural vibration control using MRE isolators.

Application in structural seismic mitigation
On the basis of the proposed field-dependent generalized model, a frequency-control strategy is developed in this part to regulate the current applied to the MRE devices in real-time for suppressing structural hazard vibrations due to earthquakes. The effectiveness of the control algorithm based on the proposed model is investigated via the numerical simulation conducted on a three-storey frame under four scaled benchmark earthquakes.

System model of MRE device-incorporated smart structure
Suppose there is a smart structure consisting of an n-level building configured with several MRE smart devices installed underneath as the isolation level. The motion equations of this (n + 1) DOFs system could be expressed as where n m denotes the number of MRE devices used as the isolation layer; s is the position vector of seismic load applied to the structure; l s is the position vector of MRE device force applied to the structure; K, C, and M are matrices of stiffness, damping, and mass, respectively, and the corresponding expressions are shown in Equations (27), (28), and (29): . . . 0 � � � À k nÀ 1 k nÀ 1 þ k n À k n 0 À k n À k n . . .

Control strategy
To avoid the predominant earthquake frequency coinciding with the fundamental frequency of a structure, i.e. the occurrence of resonance, a frequency control strategy is developed based on [38] in this section. Its fundamentals lie in that the stiffness property of the MRE device can vary among different values by switching the command current according to the control algorithm requirement. In this way, the system is endowed with real-time controllable non-resonant natural frequencies.
The entire procedure of the designed frequency control strategy includes four main processes that are measuring, anticipating, selecting and changing, which will be executed for every time interval. Firstly, the installed accelerometers measure the acceleration of seismic excitation at each instant. The earthquake signals are then forwarded to a motion anticipating analyzer to approximately predict the system's responses in every stiffness scenario. According to the forecasted results, the stiffness scenario that can achieve the smallest responses will be selected by a control decision processor based on the judgment algorithm. Finally, the control signal flow is converted to the corresponding current and sent to the MRE device to change the natural frequency of the control system.
As the core of the frequency control strategy, the motion anticipating analyzer can forecast the system's responses of next state in every stiffness scenario based on the current state. Here, the inter-storey drift, relative displacement and absolute acceleration responses of the system are adopted to assess the performance of each stiffness scenario, according to which an evaluation index can be calculated as follows: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P N f n¼1 SR ist;jsr t À n À 1 where ist denotes the type of stiffness scenario. Two stiffness scenarios are considered in this research, where ist = 1 represents the scenario without current and ist = 2 indicates the scenario with maximum current (3A). jsr corresponds to different types of system responses, where jsr = 1, 2 and 3 represent inter-storey drift, relative displacement and absolute acceleration responses, respectively. τ ist;jsr is a coefficient that stands for the maximum resonant amplitude of jsr-th response in ist-th stiffness scenario.SR ist;jsr represents the jsr-th predicted response in ist-th stiffness scenario. Δt denotes the time interval. t ist reflects the fundamental period of the system in ist-th stiffness scenario. N f ¼ t ist = 2Δt ð Þ is the number of interval samplings during half of the fundamental period.
In Equation (31), both the average and maximum responses are considered, hence, EV ist;jsr t ð Þ has relatively high robustness. To comprehensively evaluate the predicted performance of the control system in each stiffness scenario and determine the desired current for the next state, an appraisal index is defined as described in Equation (32) where t is a small-time increment. In this study, the values of both Δt and t are set to 5 ms, corresponding to the time interval of seismic excitation data.
There are two terms in Equation (32). The former reflects the ratio between the evaluation index of each system's response in every stiffness scenario and the maximum evaluation index among all the responses in all stiffness scenarios. The latter can indicate the changing tendency of the evaluation index is increasing or decreasing. The stiffness scenario achieving the smallest judgment index will be selected as the optimal state for the next instant accordingly. Therefore, the command current for the next time interval and the corresponding control force provided by the MRE device can be determined by Equation (33) and Equation (24), respectively. ; The proposed frequency control strategy can not only protect the structure via diminishing the impact of possible huge force due to the earthquake but also maintain the integrity and guarantee the function and facility of the structure. The benefit of this control method is that less power is needed to drive the MRE isolator to tune the system stiffness to achieve a non-resonant status, which is able to decrease the input of earthquake energy to the structure. Hence, it provides an ideal solution to power issues existing in practical control systems.

Case study
To assess the performance of frequency control strategy based on proposed model, a numerical study is carried out using a structural model incorporated with two MRE devices as the isolation level subjected to four benchmark earthquakes including El-Centro, Kobe, Hachinohe and Northridge earthquakes. In this case study, a threestorey shear frame model, designed by Gu [39], is employed and embedded with MRE devices. Four scenarios are considered for the comparison, including bare frame, smart frames without current input (passive-off), smart frames with maximum current input (passive-on), and smart frame with frequency control. Here, the maximum current is set to 3A, because higher current may lead to an unstable performance of MRE device. The structural parameters used for this simulation are listed in Table  2. Furthermore, since the designed gap between coils and laminated layers is 15 mm, the amplitudes of four earthquake records are scaled by 20% to confine the maximal deformation of MRE devices. Figure 24 displays the time histories of change of those two evaluation indicators together with control current of smart building with frequency control subjected to four scaled earthquake loads. Figure 25 describes the comparison of shifts of fundamental frequencies between bare frame and MRE isolators-embedded frame under different seismic excitations. Compared to other floors of the frame, the MRE isolation layer is relatively 'softer,' so the superstructure (three-storey frame) can be considered as a rigid body for the analysis. Since the maximum displacement responses always occur on first floor than second and third floors, the frequency responses in Figure 25 are the FFTs calculated from first-floor displacements of different structure systems (bare, passive-off, passive-on and frequency control) subjected to four scaled earthquakes. It can be observed that with the introduction of MRE-based isolation layer, the fundamental frequency of the frame has been changed to deviate from that of scaled earthquakes, which can effectively avoid the resonance and reduce the vibration responses of the structure. Figure 26 portrays the displacement comparison of base floor of different structural scenarios, i.e. bare structure, passive-off system, passive-on system and proposed frequency control system, subjected to four scaled earthquake excitations. It is clearly seen that the passive-off system conducts the worst performance in reducing displacement of base floor because of low stiffness. In detail, the passive-off system remarkably augments the displacements of base floor of the building during the periods of 4 s to 15 s for the case of El-Centro, 8 s to 13 s for the case of Kobe, 4 s to 8 s for the case of Northridge, and the whole excitation period for the case of Hachinohe. Especially, for the case of Kobe, the maximum base displacement of passive-off system is around 40 mm, far exceeding maximal allowable deformation of MRE device, which can lead to the failure of base isolation system in the practical application. Both passive-on and frequency control systems are effective in mitigating the displacement of base floor of the structure. Compared to passive-on system, the proposed frequency control system is capable to significantly decrease the base drift, indicating that it can be considered as a valid solution to excessive base drift problem existing in MRE-based passive isolation systems. In addition to base displacement, top floor acceleration is another commonly used indicator to evaluate the earthquake mitigation capacity of different structural systems. Figure 27 depicts time histories of top floor acceleration of different structure scenarios under four scaled earthquakes. It is noticeable that compared to the bare building, all three types of isolation systems can obviously attenuate acceleration response of top floor, not only decreasing maximal acceleration value but also keeping the acceleration of top floor at a relatively low level throughout whole loading period. Specifically, for the case of Hachinohe, the passive-off system generates an impulse at time 4.5s, which results in the top floor acceleration surpassing that of bare building. Obviously, this phenomenon is similar to that of base displacement under Hachinohe earthquake, indicating that the passive-off system fails to provide a satisfactory performance of seismic response suppression for Hachinohe earthquake. Overall, the proposed frequency control system has the best capacity in reducing structural vibration among three isolation systems.
On the other hand, maximal responses of each floor are of great interest, because they are perhaps main indicators to evaluate structural safety. Two maximal responses of different systems are chosen in this investigation, i.e. maximal interstorey drift ratio and maximal acceleration. Figure 28 shows the comparison of maximal inter-storey drift ratio against structural height of different structural scenarios under different earthquake loads, respectively. By and large, the maximal interstorey drift ratio declines with the increment of structural floor for all the seismic loads. It should be noted that the passive-off system has the largest value of peak inter-storey drift ratio at the base level for all the earthquake scenarios. For the earthquakes of El-Centro, Kobe, and Hachinohe, the values even exceed that of bare structure, which further highlights the drawback in traditional isolation systems. Even though passive-on system has smaller maximal inter-storey drift ratios of base floor than frequency control system for the cases of scaled Kobe and Hachinohe earthquakes, its overall performance is worse than that of frequency control systems for all four scaled earthquakes and even worse than passive-off system for scaled El-Centro and Northridge earthquakes. Figure 29 describes maximal acceleration of each floor of four structure scenarios under different earthquake excitations. Unlike maximal inter-storey drift ratio, the floor maximal acceleration ascends with the increment of structural floor. Apparently, the performances of bare structure and passive-off systems are worse than that of passive-on and frequency control systems in terms of acceleration attenuation. Even if the passive-on system outperforms frequency control system in reducing base acceleration for the cases of scaled El-Centro and Kobe earthquakes, the acceleration changes of different floors of the frequency control system are less significant than passive-on system. It is worth noting that the fundamental of MRE-based isolation is to decouple the structure from the base, so the base response (either inter-storey drift or acceleration) is not a primary concern compared with the responses of three floors of the frame. Moreover, the differences of peak responses of base between frequency control system and passive-on system (with optimal base responses) are very small for all four scaled earthquake cases. Accordingly, from the perspective of comprehensive analysis, the proposed frequency control strategy has the best performance in seismic mitigation among four structure scenarios. To comprehensively evaluate three controlled systems and elaborate the superiority of proposed frequency control method, multiple assessment indicators, comprising of maximum peak inter-storey drift, maximum RMS of inter-storey drift, maximum peak acceleration, maximum RMS of acceleration, maximum peak base shear and maximum RMS of base shear, are considered in this research. The relevant mathematical expressions of these indicators are shown in Equationss.   Table 3. Overall, the performance of passive-off system is better than that of passive-on system for the cases of scaled El-Centro and Northridge earthquakes, while passive-on system outperforms passive-off system under the loads of scaled Kobe and Hachinohe earthquakes. Once more, the proposed frequency-controlled system exhibits the best capacity in seismic response attenuation. Moreover, the frequency-controlled structure possesses the lowest value of peak base shear of three control systems, indicating that minimal control force is needed for semi-active control of the MRE devices incorporated smart system      Maximum RMS of base shearJ 6 ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 N P N t¼1 P 4 k¼1 m k acc c k À � 2 q ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi In this research, MRE is employed to develop adaptive vibration isolator, and the studies of modeling and semi-active control application of this innovative smart device have been comprehensively investigated for mitigating hazard vibration of structures subjected to different earthquakes. In practical applications, in addition to be installed underneath the protected structure to form the isolation system, the developed MRE devices can also be embedded between arbitrary two neighboring floors in the building to form smart storey isolation system. Such smart systems of both base isolation and storey isolation can be applied in building, bridge, electricity transmission tower, offshore platform or wind turbine, etc, to against strong wind, impact, blast, wave loads or a combination of different types of loads. Further studies are required to demonstrate these potential applications of smart MRE devices in protecting civil infrastructure.

Conclusions
This article proposed a Bouc-Wen-based phenomenological model to identify the unique behaviors of the MRE isolator. The proposed model consists of a Kevin-Voigt component and a nonlinear Bouc-Wen component connected in parallel to clarify the viscoelastic behavior and strain stiffening behavior, respectively. A numerical study that investigates the effects of model parameters on hysteresis outputs of the model was also conducted to better associate the physical meaning of each parameter with the nonlinear behavior of the MRE isolator. For the model parameter identification, an IPSO, with new individual and global optimums selection mechanism as well as nonlinear decreasing inertia weight, was developed to quickly find optimal parameter values via minimizing the errors between experimental results and model predictions. The results demonstrate that in addition to the good agreements between real and predicted shear forces, the proposed model is able to perfectly capture both strain softening and hardening phenomena caused by increasing excitation amplitude and applied current. Then, a generalized model was developed based on field-dependent parameters as the polynomial functions with current, which was validated using the experimental data of the device under random and scaled earthquake excitations. Based on this model, a frequency control strategy was proposed to tune the properties of stiffness and damping of MRE devices in real-time for structural vibration attenuation. Finally, a numerical investigation was carried out for elaborating the control performance via the comparison of four seismic loaded structural systems: bare frame, passive-off controlled frame, passive-on controlled frame and frequency controlled frame. The analysis results illustrated that the frequency control strategy based on the proposed model has outstanding ability in reducing both acceleration and inter-storey drift of structure under different seismic loads. This promising result provides a potential solution to the problem of structural vibration mitigation using MRE devices.

Highlights
• Developed a novel field-dependent phenomenological model for hysteretic responses of MRE isolator • Proposed an improved PSO algorithm for model parameter identification based experimental measurements • Designed frequency-controlled based strategy to regulate MRE isolators for structural seismic mitigation • Significantly diminished structural responses in terms of inter-storey drift, acceleration and base shear.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by the Australian Research Council.