Magneto-Sensitive Rubber in a Vehicle Application Context – Exploring the Potential

The application of magneto-sensitive (MS) rubber in a vehicle vibration control area is likely to be expected. This conclusion is based on the following two reasons: the maturity of fabrication of MS rubber which meets the application requirement and the feasibility of the constitutive model of MS rubber that accurately reflects its mechanical performance. Compared with the traditional rubber, small ferromagnetic particles are embedded in the elastomer of MS rubber, leading to a change of mechanical properties when an external magnetic field is applied. Therefore, devices with MS rubber, can be viewed as a semi-active actuator. In this study, MS rubber with a relative high increase in the magneto-induced modulus is fabricated and characterized. Furthermore, a one-dimensional constitutive model to depict the magnetic field-, frequency-, and strain amplitude-dependent dynamic modulus of MS rubber is applied. Finally, simulations of a MS rubber semi-active suspension under a bump and a random ground excitation with different control strategies on a quarter vehicle model are conducted to illustrate the feasibility of the MS rubber in the vehicle vibration control application context.

The application of magneto-sensitive (MS) rubber in a vehicle vibration control area is likely to be expected. This conclusion is based on the following two reasons: the maturity of fabrication of MS rubber which meets the application requirement and the feasibility of the constitutive model of MS rubber that accurately reflects its mechanical performance. Compared with the traditional rubber, small ferromagnetic particles are embedded in the elastomer of MS rubber, leading to a change of mechanical properties when an external magnetic field is applied. Therefore, devices with MS rubber, can be viewed as a semi-active actuator. In this study, MS rubber with a relative high increase in the magneto-induced modulus is fabricated and characterized. Furthermore, a one-dimensional constitutive model to depict the magnetic field-, frequency-, and strain amplitude-dependent dynamic modulus of MS rubber is applied. Finally, simulations of a MS rubber semi-active suspension under a bump and a random ground excitation with different control strategies on a quarter vehicle model are conducted to illustrate the feasibility of the MS rubber in the vehicle vibration control application context.

INTRODUCTION
The performance of vehicle suspensions is closely related to the automotive noise, vibration, and harshness (NVH) issue (Heißing and Ersoy, 2011). In order to improve the ride comfort and the road holding ability and to reduce suspension deflection, three types of suspensions, namely, passive, active, and semi-active, are utilized. The passive suspension has the advantage of cost-effectiveness and design simplicity. However, its stiffness cannot be adjusted according to the external loadings. Therefore, the performance of the passive suspension is limited. Regarding the active suspension, it can exert control force independently without relying on the structural response, which makes it effective to reduce vibration. Nevertheless, its high cost, design complexity, and the lack of stability when the power input is cut off under extreme conditions prohibit its practical application (Martins et al., 1999). The recent progress in smart material development makes the semiactive vehicle suspension to have a superior suspension measure. In particular, desired vibration-reduction performance is more likely to be achieved by semi-active suspensions, compared to the results using passive suspensions. Furthermore, less energy is consumed by semi-active suspensions when compared with active ones. In the past few years, magneto rheological (MR) fluid (MRF) Biglarbegian et al., 2008;Wilson and Abdullah, 2010;Phu and Choi, 2019;Du et al., 2020) was widely used in the semi-active control area due to its rapid response, temperature insensitivity of the mechanical property, and its high MR effect. However, it is worth noting that although the motion of the suspension can hinder the settlement of ferromagnetic particles in MRF, the increased zero-field viscosity due to the wear of ferromagnetic particles after a long-term operation will reduce the vibration control effect of MRF suspensions (Hu et al., 2012). For magneto-sensitive (MS) rubber, ferromagnetic particles are fixed in the elastomer matrix, thus avoiding the problem of ferromagnetic particle wear and material leakage (Bastola et al., 2019). Therefore, semi-active suspensions based on MS rubber provide a feasible and practical solution to the vibration control problem for vehicles.
The research of MS rubber can be divided into three categories, including fabrication optimization to improve the mechanical property, constitutive modeling of MS rubber, and exploring of possible applications of MS rubber. Davis (1999) analyzed that the optimum particle volume fraction for MS rubber to achieve the highest magnetic modulus is 27%. This conclusion was verified experimentally by Lokander and Stenberg (2003) and Zhou (2003). Subsequently, the influence of the matrix type (Shen et al., 2004), coupling agent (Chen et al., 2007a), pre-structuring (Li et al., 2008), and the surface modification  of particles on the MR effect and damping property of MS rubber were investigated. Tong et al. (2018) fabricated a kind of MS rubber with an improved fieldinduced storage modulus by mixing flow-like particles into the elastomer matrix. It can be foreseen that it will be possible to fabricate MS rubber with a good mechanical performance and a high relative MR effect to meet the engineering application needs, due to the rapid development of MS rubber, including the optimization of their manufacturing process.
In order to predict the mechanical performance of MS rubber, constitutive models with strain and magnetic field strength as inputs and stress as output are needed. Initially, Jolly et al. (1996) developed a magnetic dipole model to simulate the shear modulus increase of MS rubber under a magnetic field. The dipole model was further extended by taking the viscoelastic behavior of MS rubber (Chen et al., 2007b), the interaction of multi-chains (Zhu et al., 2006), and the distribution (Zhang et al., 2010) of chain angles into account. In parallel, constitutive models based on the rheological theory were developed to depict the magnetic fieldrelated non-linear viscoelastic behavior of MS rubber (Li et al., 2010;Chen and Jerrams, 2011). Wang et al. (2018) developed a constitutive modeling-electromagnetic analysis combining method to describe the mechanical behavior of MS rubber. By using the methods mentioned, the mechanical behavior of MS rubber under a small deformation can be predicted accurately. In order to predict the mechanical performance of MS rubber under a large deformation, Dorfmann and Ogden (2003) as well as Kankanala and Triantafyllidis (2004) derived the fundamental balance principles of magneto-elastic materials based on the continuum mechanical theory. Subsequently, the continuum mechanical-based model was further expanded to include the viscoelasticity (Saxena et al., 2014), anisotropy (Bustamante, 2010;Danas et al., 2012), and the strain amplitude-dependent non-linear behavior of MS rubber (Wang and Kari, 2020). Based on the work of Kankanala and Triantafyllidis (2004) and Zhao et al. (2019) developed a model for hard-magnetic particle-based MS rubber. However, for the application of MS rubber in the field of semi-active vibration control, normally, the required force and structural responses are known values where the magnetic field required to achieve the control force needs to be calculated. Accordingly, an inverse model with stress and strain as inputs and the magnetic field strength as output is needed. Current research (Jung et al., 2011;Yang et al., 2016;Wang and Kari, 2019) showed that a better vibration reduction effect can be achieved by MS rubber-based devices, as compared to the traditional passive ones, by controlling the stiffness of MS rubber-based devices discontinuously through switching the magnetic field between zero and saturation state (on-off control). Compared with the onoff control strategy, the development of the inverse model makes the continuous adjustment of the stiffness of the MS rubber-based device a reality. Therefore, it can be foreseen that a semi-active control strategy based on the inverse model of MS rubber can achieve an even better vibration control effect than the previously mentioned on-off control strategy.
To explore the application of MS rubber, Liao et al. (2013) and Sun et al. (2018) designed MS rubber vibration absorbers with the function of tracking external excitation frequency and verified the vibration control effect experimentally. Wang et al. (2017) designed a conical shaped MS rubber-tuned mass damper with an improved stiffness changing ability. Blom and Kari (2012) designed an MS rubber bushing and demonstrated that an optimal isolation effect can be obtained by changing the magnetic field under different frequencies. Subsequently, Alberdi-Muniain et al. (2012 verified the previous conclusion experimentally. Wang and Kari (2019) simulated the vibration control effect of a MS rubber isolator under a random loading case by using a fuzzy logical control algorithm. Jin et al. (2020) explored the possible application of MS rubber in highspeed railway vehicles. Fu et al. (2019) investigated the design method of an adaptive fuzzy controller for a MS rubber vibration isolator under time-varying sinusoidal excitation. In the field of automobile suspension, Du et al. (2011) deduced the design method of the control gain of the MS rubber semi-active seat suspension H-infinity controller under the conditions of stiffness saturation and parameter uncertainty. Liu et al. (2020) used an adaptive neutral network control design method to simulate the control effect of the MS rubber-based seat isolator under a bump and a random road excitation. However, the magnetic fielddependent non-linear viscoelastic behavior of MS rubber has not been fully taken into account. Meanwhile, Liu et al. (2020) applied a polynomial expansion fitting method to determine the relationship between the shear modulus of the MS rubber isolator and the magnetic field strength. However, the method is only suitable for some specific working conditions which is difficult to extend directly to practical application situations.
In this study, a silicone rubber-based MS rubber with a relatively high MR effect is fabricated and characterized on a commercial rheometer. Subsequently, a constitutive model which depicts the magnetic field-dependent non-linear viscoelastic behavior of MS rubber is applied to depict the mechanical performance of MS rubber. Furthermore, an inverse model based on the constitutive model of MS rubber with stress and strain as input and the magnetic field strength as output is proposed. On the basis of the constitutive and inverse model, the H-infinity control strategy is used to simulate the vibration control effect of the MS rubber-based vehicle suspension under a bump and a random road excitation.
The research in this study includes material fabrication, characterization, forward and inverse modeling, and simulation analysis, covering a broad application range for MS rubber in vibration control area. The research performed in this article is meaningful for advancing the application of MS rubber in the vehicle suspension area. Furthermore, the forward constitutive model and the inverse model proposed can also be used in other semi-active control scenarios.

Material Fabrication and Characterization
The MS rubber sample in this study was fabricated by mixing silicon rubber (type multi-purpose silicone sealant, Selleys Pty Ltd., Australia), silicon fluid (XIAMETER PMX-200, Dow Corning GmbH, United States), and carbonyl iron particles (CIPs, type CN, BASF, Germany diameter 7 µm on average) with a mass fraction ratio of 2:1:7 at room temperature. After mixing for 5 min, the mixture was placed in a vacuum chamber with a pressure of 0.06 MPa for 10 min to remove air bubbles. Subsequently, the mixture was put into a mold under a magnetic field strength of 1 T for 30 min. After the CIPs were chained, the mixture was left for 48 h at room temperature for post-curing. The images of the MS rubber sample under the scanning electron microscope (Gemini 500, Carl Zeiss, Jena, Germany) at different magnification rates are in Figure 1, where arrows represent the chain direction of the iron particles. It can be seen that most iron particles are distributed along the chain direction in the matrix and that there is no large-scale clustering of the iron particles.
After MS rubber was fabricated, the sample with a diameter of 20 mm and a thickness of 2.12 mm were put into a commercial rheometer (Physica MCR 301, Anton Paar Co., Austria) to test the dynamic performance of MS rubber. The test was performed at room temperature 22 ± 1 • C. Throughout the test, the external magnetic field was parallel to the iron particle chain direction, while the twist direction of the commercial rheometer was perpendicular to the chain direction. During the measurement, three kinds of conditions were altered: strain amplitude (1, 2, and 5%), magnetic field strength (0, 0.17, 0.34, 0.51, 0.68, and 0.83 T), and frequency (1,2,5,10,20,30,40,50,60,70,and 80 Hz). All combinations of the above three types of conditions were measured, recorded, and averaged over 20 periods. In order to avoid the influence of the Mullins effect (Mullins, 1969), a repeated shear strain where the strain amplitude is slightly larger than the maximum test shear strain was applied before the measurements started. The repeated preloading was terminated when the test results of two consecutive cycles are the same.
Test results for the magnitude and loss factor of the shear modulus vs. frequency with different strain amplitudes and magnetic field strengths are shown in the dot symbols in Figures 2-4. Similar to the measurement result by Blom and Kari (2005), the measured shear modulus of MS rubber shows a magnetic field, frequency, and strain amplitude dependence. Specifically, a higher frequency, a stronger magnetic field strength, and a lower strain amplitude led to an increased shear modulus magnitude, while the loss factor almost keeps constant under different frequencies, magnetic field strengths, and strain amplitudes. Furthermore, it can be found that the increasing rate of the shear modulus magnitude, with respect to the increasing magnetic field strength, slows down gradually until a magnetic saturation is reached for the material. Although the shear modulus magnitude of MS rubber seems to increase after 0.83 T, the upper limit of the magnetic field strength in the test was nevertheless set to be 0.83 T due to the sharp rise of the test object temperature that is generated by even stronger magnetic fields.

Constitutive Modeling
An elastic spring with a parameter µ e , a viscoelastic element with parameters a and b, and a frictional element with parameters τ fmax and γ 1/2 in parallel are used to depict the dynamic mechanical behavior of MS rubber. Initially, a fractional derivative viscoelastic element is utilized to depict the frequencydependent shear modulus of MS rubber and to simplify the parameter identification. The relation between the viscoelastic stress τ ve and strain γ for the fractional derivative element in frequency domain isτ where (·) represents the Fourier transform, j is the imaginary unit, and ω is the angular frequency. From Equation (1), it can be observed that the modulus magnitude bω a increases with increasing frequency and the loss factor tan (aπ/2) is a constant value for the fractional derivative element. For more details regarding the fractional derivative model, please refer to Lubich (1986), Lion and Kardelky (2004), and Kari (2017Kari ( , 2020. In order to simulate the strain amplitude dependence of the shear modulus of MS rubber, a smooth frictional stress model is used (Berg, 1998). The relation between the friction stress (τ f ) and  strain γ is where τ fmax and γ 1/2 are model parameters, representing the maximum friction stress and the strain where one half of τ fmax is reached andγ represents the strain rate. Initially, τ fs = 0 and γ s = 0. When the loading direction is changed, τ fmax and γ 1/2 are updated (ifγ = 0, τ fs = τ f |γ =0 and γ s = γ |γ =0 ). For more details regarding the smooth frictional stress model, the reader is referred to Berg (1998). By Equations (1) and (2), the total dynamic modulus G * in frequency domain is After establishing the basic frame of the constitutive model, a non-linear least square fit method is implemented to obtain the basic material parameters by using the measurement results at zero magnetic field. The function lsqnonlin in MATLAB R (MATLAB Release 2015b, The MathWorks, Inc., Natick, Massachusetts, United States) is utilized for parameter identification. In the process of parameter identification, the objective function is set as where G * mea (i) and G * sim (i) are the measured and simulated shear modulus of MS rubber, respectively. The symbol i represents the number of loading cases ranging from 1 to 33 since there are 3 kinds of strain amplitude and 11 kinds of frequency during the test under zero magnetic field. The identified material parameters are µ 0 e = 3.4587 × 10 4 N m −2 , a = 0.2071, b = 3.4587 × 10 4 N s a m −2 , τ 0 fmax = 3.2583 × 10 2 N m −2 , and γ 0 1/2 = 2.6820 × 10 −2 , where symbols with a superscript zero represent the material parameters at zero magnetic field.  Compared with the generalized Maxwell model, less material parameters are needed for the fractional derivative element to depict the viscoelastic behavior of rubber materials (Haupt and Lion, 2002). Thus, the material parameters are easier to identify for the fractional derivative element. The solution of the viscoelastic stress of the fractional derivative model at the current time step τ ve (n t) in time domain is where c 1 (α) = 1, c j (α) = j − 1 − α c j−1 (α) j and t is the time interval for the numerical implementation. As a comparison, the solution of the viscoelastic stress of the Maxwell element at the current time step (Simo and Hughes, 2006) is where G and η are the modulus and the relaxation time of the Maxwell model, respectively. The equivalence of Equation (6) Frontiers in Materials | www.frontiersin.org By comparing Equations (5) and (6), it can be seen that in order to obtain the viscoelastic stress at the current time step, only stress and strain data at the previous time step are needed for the Maxwell model. However, data of all strain history are needed to obtain the current viscoelastic stress for the fractional derivative element. Therefore, for numerical implementation concern, the Maxwell model is more effective than the fractional derivative element as the time step increases. Taking into account the need for the subsequent time domain simulation, three Maxwell models in parallel are used to substitute the previous fractional derivative model. The identified parameters τ 0 fmax and γ 0 1/2 are fixed for the second round of parameter identification. The new material parameters identified are µ 0 e = 6.7583 × 10 4 N m −2 , G 1 = 1.8256 × 10 4 N m −2 , G 2 = 3.1283 × 10 4 N m −2 , G 3 = 6.2925 × 10 4 N m −2 , η 1 = 1.8256 × 10 −1 s, η 2 = 1.5884 × 10 −2 s, and η 3 = 2.3592 × 10 −3 s with a relative error between the measurement and simulation result of 2.769%. The comparison between the measurement and simulation results are in Figure 5.
After modeling the frequency and the strain amplitude dependence, the magnetic field dependence of the shear modulus of MS rubber depicted by a hyperbolic tangent function is included in the constitutive model. According to Berg (1998), the area of the hysteresis loop for the smooth frictional stress model increases with increasing τ fmax and decreasing γ 1/2 . Following the path of Blom and Kari (2011), the magnetic field dependence is included in τ fmax and γ 1/2 by and where B s is a material parameter used to reflect the magnetic saturation behavior of MS rubber and δ 1 and δ 2 are parameters to reflect the magnetic enhancement of MS rubber. To guarantee that the loss factor is a relatively constant value with respect to the magnetic field strength, the hyperbolic tangent function is included into the elastic part as well where the meaning of δ 3 is similar to δ 1 and δ 2 . The third round of parameter identification was conducted and the magneticrelated parameters are B s = 0.4263 T, δ 1 = 16.2172, δ 2 = 0.6538, and δ 3 = 2.0776, with a relative error of 8.707 %. The comparison between the simulation and measurement results are in Figures 2-4. It can be found that the magnetic field-, frequency-, and the strain amplitude-dependent modulus of MS rubber can be well depicted by the model developed. Since three Maxwell elements without magnetic dependency are utilized to simulate the frequency dependency of MS rubber, there is a certain degree of deviation between the simulation and measure results at high frequency and high magnetic field. Due to the wide distribution of the relaxation time of polymer material in the time domain and the possible particle-elastomer interaction within MS rubber, more Maxwell elements with magnetic dependency are favorable to further improve the fitting accuracy between the simulation and experimental results.

Bisection Method-Based Inverse Model of MS Rubber
As a semi-active actuator, the desired control force cannot be applied by MS rubber suspension directly. The implementation of the control force depends on the interaction between the structure response and the controlled magnetic field. Therefore, it is of great importance to obtain the required magnetic field by the desired control force and the structure response. However, for the constitutive models developed in section 2.2, there is no corresponding direct inverse model with an explicit form. In order to obtain the magnetic field, methods using polynomial fitting of the constitutive models (Choi et al., 2001;Meng et al., 2018) and using an adaptive neuro-fuzzy inference system (Wang et al., 2011;Zong et al., 2012;Yang et al., 2019) are often used. The magnetic field applied to the MS rubber device in practical applications has a lower limit (zero magnetic field strength) and an upper limit (saturated magnetic field strength). Therefore, a bisection-based method can be used to obtain the required magnetic field strength. The the constitutive model developed in section 2.2 can be expressed as The difference between the target stress and the actual stress is expressed by a function Then, the bisection-based algorithm can be applied to obtain the required magnetic field. Firstly, the strain response γ , the strain rateγ , stress needed τ needed , and the tolerance ε are known. The upper limit of the magnetic field B upper is set to be 0.83 T and the lower limit of the magnetic field B lower is set to be 0 T. Next, if then B required is selected as the value of f γ ,γ , B upper and f (γ ,γ , B lower ) that is closest to τ needed . Else, the tolerance of the magnetic field B tol is set to be If |B tol | is larger than the tolerance ε, then After updating the magnetic field upper limit, if else, After the upper and the lower limit of the magnetic field is updated again, the value of B required is compared with ε. If the convergence criteria is not met, new rounds of iterations from Equations (15) to (18) are applied until the convergence criteria is met. Finally, for the case when the condition in Equation (13) is met, B applied is selected as the value of f γ ,γ , B upper and f (γ ,γ , B lower ) that is closest to τ needed ; otherwise, B applied = B upper . In order to verify the bisection-based inverse model, a Gaussian white noise-based strain signal, ranging from 0 to 25 Hz, with an RMS value of 0.03 m and a time duration of 2 s, and a Gaussian white noise-based magnetic field signal, ranging from 0 to 100 Hz, with magnitudes vary from 0 to 1.5 T, are applied. By the constitutive model developed in section 2.2, the total stress (target stress) corresponding to the strain and magnetic field can be obtained. Subsequently, by the bisectionbased inverse model, the predicted magnetic field and the corresponding predicted stress can be obtained. The comparison of the target and predicted values (magnetic field and stress) by the constitutive model and the inverse model is given in Figure 6. The overlapping between the target and predicted magnetic field/stress in Figure 6 demonstrates the effectiveness of the bisection-based inverse model. The values of the target magnetic field, which is larger than 0.83 T cannot, be tracked by the inverse model, since the upper limit of the magnetic field strength is set to be 0.83 T.

Vehicle Suspension Model and Formulation of H-Infinity Controller
A quarter-car suspension model consists of a one-fourth car body, suspension components (spring plus MS rubber-based semi-active actuator), a wheel, and a tire are investigated in this study. For more details of the quarter-car suspension model, please see Figure 7. Symbols z s , z t , and z g represent the car body, wheel, and the ground displacement, respectively. The symbol u represents the control force that can be adjusted by the MS rubber suspension through changing the magnetic field. The model parameters for the quarter-car suspension model are m s = 504.5 kg, m t = 62 kg, k s = 10,100 N m −1 , c s = 200 N s · m −1 , and k t = 252,000 N m −1 .
In order to increase the driving performance, it is required that the transfer function from the road disturbance to the vertical acceleration of the car body (z s ), the tire deflection (z t − z g ), and the suspension deflection (z s − z t ) should be small enough. Therefore, the controlled output is defined as z = [z s z s − z t z t − z g ] T . The measured output is defined as y = [ z t − z gżs ] T due to the fact that the suspension deflection (z t − z g ) and the car-body velocity (ż s ) can be directly measured in practice. The dynamic equilibrium equation for the quarter-car system in state space form is where FIGURE 7 | Schematic configuration of a quarter-car system with a MS rubber semi-active suspension.
and w =ż g .
The parameters α and β in matrix E are weighting factors for the suspension and tire deflection, respectively. By tuning the weight factors, the trade-off among the car-body acceleration, suspension deflection, and the tire deflection can be controlled. The relation between the controlled force u and the measured output y is where K = [ k 1 k 2 ] is the static output feedback control gain to be designed. To reduce the object transfer functions, robust control algorithms, such as LQR (Zhang et al., 2008;Zhang and Zhuan, 2020), H-2 (Li et al., 2003;Shukla et al., 2016), and Hinfinity (Xie et al., 2004;Alma et al., 2011;Khot et al., 2017) controller, are often utilized. Compared with other controllers, the H-infinity norm measure the energy-to-energy gain between the external disturbance and the target outputs directly, which can be viewed as the worst-case gain in frequency domain. In this study, the H-infinity index is used for the controller design. According to Du et al. (2003), the system in Equation (19) is asymptotically stable with a disturbance attenuation γ > 0 if and only if there exist matrices G = G T > 0 and K such that where sym(·) represents symmetric term, ( * ) denotes the symmetry matrix block, superscript T denotes matrix transpose, and I is the unit matrix. Normally, the bi-linear matrix inequality problem in Equation (29) cannot be solved directly. According to Ebihara et al. (2015), a coordinate-descent algorithm can be used to transfer the bi-linear matrix inequality problem in Equation (29) to several iterative linear matrix inequality problems, which can be solved in MATLAB LMI package directly. The dimensionless weighting factor α and β in Equation (24) are set to be 4 and 10, respectively. After applying the coordinatedescent algorithm-based iterative method, the H-infinity norm and the static output feedback control gain are obtained with γ = 8.441 and K = [ 8.9459 N m −1 2.9522 N s · m −1 ] × 10 3 . The corresponding transfer functions from the road disturbance to the car acceleration (z s ), suspension deflection (z s − z t ), and the tire deflection (z t − z g ) are in Figure 8. Compared with the open loop case, after applying H-infinity control algorithm, the peak gain for the car acceleration, suspension deflection, and the tire deflection are all reduced.

Design of MS Rubber Suspension and Numerical Implementation Method
The whole flowchart corresponding to the semi-active control of the MS rubber-based vehicle suspension system is given in Figure 9. After the external disturbance excitation is applied, the measured outputs y in Equation (19) is transmitted to the Hinfinity controller as the input for the controller. Subsequently, a desired control force output is calculated by the H-infinity controller. However, the desired control force cannot be exerted directly. Therefore, the MS rubber inverse model in section 2.2 is utilized to determine the magnetic field corresponding to the desired control force. Finally, the actual force by the MS rubber suspension is applied to the quarter car system and new measured outputs are detected. In order to verify the vibration control performance of the MS rubber semi-active suspension, simulations under a bump and a random ground excitation are conducted. Before the simulation, the dimension of the MS rubber-based suspension is designed. To reduce the static deformation of the MS rubber suspension caused by the preloading of the car body mass, a linear spring with a stiffness of 2,100 N m −1 is installed in parallel with the MS rubber-based suspension. Sun et al. (2018) applied a multi-layer strategy for the MS rubber suspension in order to improve the deformation ability of the whole suspension, meanwhile ensuring the strain of each layer of MS rubber is small. Considering the simplest case, a multi-layer sandwich typed structure made of steel plate and MS rubber working in a shear mode are utilized. The total layer number of MS rubber is 10. For each layer MS rubber, the thickness is 0.04 m and the area is A = (0.10 m) 2 . Therefore, the approximate initial stiffness of the suspension by using the shear modulus of MS rubber at 5% strain amplitude under 1 Hz is k s ≈ 2 0.10 2 × 8 × 10 4 5 × 0.04 + 2100 N m −1 = 10100 N m −1 .
Frontiers in Materials | www.frontiersin.org To illustrate the effectiveness of the control strategy, MS rubber suspension with the maximum stiffness (passive-on), MS rubber suspension with the minimum stiffness (passive-off), sky-hook on-off control strategy, and the active control strategy are simulated as well. The corresponding magnetic field control strategy for the sky-hook on-off control strategy is In order to fully consider the magnetic field-, frequency-, and the strain amplitude-dependent dynamic behavior of MS rubber and obtain the structure response, the Newmark-beta method (Géradin and Rixen, 2014) based on implicit time integration is applied. The dynamic equilibrium equation of the whole system can be written as where and First, an initial accelerationq trial is obtained by the known response at the previous time step, q n−1 , where symbols with subscript n − 1 and n represent the response at the previous and the current time step, respectively. Subsequently, the predicted displacement q, velocityq, and accelerationq at the current time step are obtained by the average acceleration method Following these, the residual force R is obtained A threshold for the residual force κ = 1 × 10 −3 N is set as the criteria for convergence. If R is larger than κ, a correction is applied by   where q = −R/S, with S as the Jacobian matrix for the system After q is updated, the new residual force R is compared with the threshold κ. If the converge criteria is not satisfied, new rounds of iterations are applied until the convergence criteria is met. Finally, the structure response at the current time step can be obtained. Regarding ∂F ∂q, a finite difference numerical method can be applied to obtain its value .

Comparison on Bump Response
According to Sun et al. (2020), a road bump shock can be represented by where H = 0.04 m, L = 5 m, and v= 25 m s −1 are the bump height, bump length, and the passing car velocity, respectively.
After applying the numerical simulation analysis by the numerical methods in sections 2.2, 3.1, and 4.1 , the comparison of the passive-on, passive-off, H-infinity semi-active, sky-hook on-off, and the active control strategies of the car-body acceleration, suspension deflection, and the tire deflection is given in Figure 10. Compared with passive-on and passive-off cases, a better vibration control effect can be achieved by the H-infinity semi-active, sky-hook on-off, and the active control strategies. The root mean square (RMS) values for the target responses (z s , z s −z t , and z t −z g ) under different control strategies for the bump excitation are in Table 1. It can be found that the vibration control effect is the best for the active control case, followed by the H-infinity semi-active control and the sky-hook on-off control strategies. In Figure 11, the desired control force for the active control strategy and the actual applied force by the MS rubber semi-active suspension vs. suspension deflection are compared. The first and third quadrants in the figure correspond to the positive stiffness area, while the second and forth quadrants correspond to the negative stiffness area. Due to the positive shear modulus and the limited MR effect of MS rubber, the advantage of the active control strategy cannot be fully utilized by the MS rubber suspension. This conclusion can be verified by the partially overlapping between the desired control force-deflection curve and the actual applied force-deflection curve.

Comparison on a Random Response
After conducting the simulation under the bump response, a new round of simulation analysis under a random response is conducted. According to Sun et al. (2020), the road roughness can be treated as a random process with a given displacement power spectral density under frequency domain G q (f ) by where G q (n 0 ) is the road roughness coefficient, n 0 = 0.1 m −1 and v is the car velocity. In the simulation, C road roughness is applied, where G q (n 0 ) = 64 × 10 −6 m 3 . The car velocity is set to be v = 60 km h −1 . The corresponding road vibration signal  Sky-hook on-off 3.5575 × 10 −2 7.3842 × 10 −3 6.0066 × 10 −3 generated by a fast Fourier transformation method (Beaulieu and Tan, 1997) with a sampling frequency of 1,000 Hz and a time length of 5 s is applied to the quarter car system. The comparison of the passive on, passive off, H-infinity semiactive, sky-hook on-off, and the active control strategies of the car-body acceleration, suspension deflection, and the tire deflection is given in Figure 12. The magnetic field vs. time under the H-infinity semi-active strategy and the sky-hook on-off control strategy is given in Figure 13. The result of the RMS values of the targeted responses are in Table 2. It can be found that compared with passive off and passive on cases, a better vibration reduction effect can be achieved by the H-infinity semi-active control strategy in controlling the car-body acceleration and the suspension deflection. Strangely, according to Table 2, the RMS value of the suspension and the tire deflection under the passive-on control strategy is the smallest. This can be attributed to the increase of the overall stiffness of the quarter-car system after implementing the maximum stiffness of the MS rubber suspension. Therefore, the deflection response under external excitation is smaller under the passive on case compared to other cases. However, the reduction of the car-body acceleration cannot be achieved directly by simply increasing the overall stiffness of the quartercar system. This conclusion can be verified by the comparison of the car-body acceleration results under different control strategies where the acceleration under the passive on case is the largest. Furthermore, it can be found that there exists a flutter phenomena in the car body acceleration for the sky-hook on-off control strategy compared with the H-infinity control strategy. Due to the continuous adjusting of the magnetic field, the vibration control effect of the H-infinity semi-active control strategy is superior as compared to the sky-hook on-off control strategy.

CONCLUSION
The fabrication, characterization, and the modeling of the mechanical performance of MS rubber along with the application of MS rubber in vehicle semi-active suspension are studied in this research. The matching between the simulated and measured shear modulus of MS rubber demonstrates that the applied constitutive model can depict the magnetic field-, frequency-, and strain amplitude-dependent shear modulus of MS rubber with accuracy. In order the desired control force can be implemented by the MS rubber semi-active device, an inverse model is developed to obtain the magnetic field corresponding to the desired control force. The matching between the target stress and the predicted stress generated by the calculated magnetic field through the inverse model indicates that the inverse model developed can track the desired control force with high accuracy. Furthermore, in order to explore the possible application of MS rubber in vehicle vibration control, a MS rubber semiactive suspension with the H-infinity control strategy applied to a quarter car model under a bump and a random ground excitation is studied. The results reveal that compared with the passive and the sky-hook on-off control cases, a better vibration reduction effect can be achieved by the MS rubber semi-active suspension with the H-infinity control strategy. However, in order to implement MS rubber in the vehicle vibration control area and achieve a good vibration reduction effect, issues such as further increase the MR effect of MS rubber, modeling the mechanical behavior of MS rubber with considering preload effect and the design optimization of the magnetic field generation circuit should be investigated. Besides, an increased vibration reduction effect is expected while adding an adjustable damping component to the current MS rubberbased suspension since the former adapts its damping while the latter adapts its stiffness. Furthermore, as the deflection of the suspension increases, the MR effect of MS rubber will weaken. Although the constitutive model can reflect the strain dependence of the MS rubber modulus, in practical applications, it is necessary to optimize the design of the MS rubber suspension to ensure that the MR effect of MS rubber can be utilized to the ultimate extent and achieve a good vibration control effect.

AUTHOR CONTRIBUTIONS
BW developed the theoretical formalism and performed the numerical simulations. TH helped to carry out the fabrication and characterization of MS rubber. XG and LK supervised the project and provided revisions to the scientific content of the manuscript. LS, JL, and ZX provided funding. All authors contributed to manuscript revision, read, and approved the submitted version.