Fractional Modeling and Analysis of Coupled MR Damping System

288 IEEE/CAA JOURNAL OF AUTOMATICA SINICA, VOL. 3, NO. 3, JULY 2016 Fractional Modeling and Analysis of Coupled MR Damping System Bingsan Chen, Chunyu Li, Benjamin Wilson, and Yijian Huang Abstract—The coupled magnetorheological (MR) damping system addressed in this paper contains rubber spring and magnetorheological damper. The device inherits the damping merits of both the rubber spring and the magnetorheological damper. Here a fractional-order constitutive equation is introduced to study the viscoelasticity of the combined damper. An introduction to the definitions of fractional calculus, and the transfer function representation of a fractional-order system are given. The fractional-order system model of a magnetorheological vibration platform is set up using fractional calculus, and the function of displacement is presented. It is indicated that the fractional-order constitutive equation and the transfer function are feasible and effective means for investigating of magnetorheological vibration device. Index Terms—Fractional calculus, magnetorheological (MR) fluid, fractional-order constitutive equation, fractional-order system, system modeling. I. I NTRODUCTION M AGNETORHEOLOGICAL (MR) fluids are particulate suspensions whose rheological properties are dramatically altered by magnetic fields. In shear flow, an applied magnetic field can increase the apparent viscosity by several orders of magnitude. This phenomenon is currently being exploited in commercial applications. MR dampers are a new research development in the field of semi-active control. The mechanical model of an MR fluid is a key way to reach the ideal control effect of the device. In fact, the mechanical properties of MR fluids and their dampers are also influenced by many factors including the vibration displacement, the acceleration, the vibration frequency among other factors. The dynamics of an MR damper can be described through both theoretical and empirical relationships. Manuscript received September 21, 2015; accepted February 1, 2016. This work was supported by National Natural Science Foundation of China (51305079), Natural Science Foundation of Fijian Province (2015J01180), Outstanding Young Talent Support Program of Fijian Provincial Education Department (JA14208, JA14216), and the China Scholarship Council. Recommended by Associate Editor YangQuan Chen. Citation: Bingsan Chen, Chunyu Li, Benjamin Wilson, Yijian Huang. Fractional modeling and analysis of coupled MR damping system. IEEE/CAA Journal of Automatica Sinica, 2016, 3(3): 288−294 Bingsan Chen is with the School of Mechanical and Automotive Engineering, Fujian University of Technology, Fuzhou 350118, China and the Department of Chemical and Biological Engineering, University of Wisconsin, Madison, WI 53705, USA (e-mail: [email protected]). Chunyu Li and Yijian Huang are with the School of Mechanical and Automotive Engineering, Fujian University of Technology, Fuzhou 350118, China (e-mail: chunyuli [email protected]; [email protected]). Benjamin Wilson is with the Department of Chemical and Biological Engineering, University of Wisconsin, Madison, WI 53705, USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Stanway[1−2] established a rational mechanics model based on MR fluids viscosity. The Stanway model contains Coulomb friction and viscous damping, but the elastic characteristic of the MR fluids is not included; Zhou and Qu[3] modified the Bingham model based on a constitutive relation for MR fluids, the precise calculation of mechanical characteristics is given, but the model is inconvenient due to its many parameters; Gamota and Filisko[4] also proposed a similar viscoelasticplastic mechanical model. In this paper, the viscoelastic model of the MR damper is established by fractional calculus. As the physical meaning of fractional calculus is not clear, not achieving its genetic characteristics and infinite memory function, so its practical engineering application is latter than the integer order calculus, although they were present almost at the same time. Fractional calculus has been introduced into rheology by Slonimsky[5] and Friedrich[6] , et al., to study the nonlinear constitutive relation. Considerable progress has been made in using fractional calculus to study nonlinear viscoelasticity. Bagley and Torvik[7] used fractional calculus to study the three- dimensional constitutive relation as well as find limits of the model parameters caused by the thermodynamic effects. Paggi et al.[8] modeled the thermoviscoelastic rheological behavior of ethylene vinyl acetate (EVA) to assess the deformation and the stress state of photovoltaic (PV) modules and their durability; Jóźwiak et al.[9] studied the dynamic behavior of biopolymer materials with fractional Maxwell and Kelvin-Voigt rheological models. Fractional calculus has been a breakthrough in the theory and application of the constitutive equation, and emerged as a new principle and method for the constitutive equation of viscoelastic materials. Therefore, the constitutive equation applying fractional calculus theory of viscoelastic materials is always one key research field. In this paper, the fractional calculus is introduced to explore the viscoelastic properties of the composite MR-rubber damper, and the mechanical properties of the composite are also studied. The dynamic characteristics of the composite damper are verified by experiments, which provide the practical basis for verification of the theoretical results on MR shock absorber. II. M ODEL E STABLISHMENT A. Fractional Order Model The fractional order derivative rheological model is based on the spring, dashpot and friction element. CHEN et al.: FRACTIONAL MODELING AND ANALYSIS OF COUPLED MR DAMPING SYSTEM As shown in Fig. 1 (a), for a = 0, the model is a typical Hook theorem, given by (1). σ(t) = τ 0 EDt0 ε(t), (1) where σ(t) represents applied stress, E is the elastic modulus, Dt0 ε(t) is the 0 order time derivative with respect to t of the strain ε(t). 289 The advantage of the Caputo fractional calculus definition is that the physical meaning of the initial value is the same as integer order calculus. So for an arbitrary real number p, the definition of fractional calculus is given by dn (Dp−n f (t)), 0 < n − p < 1, dtn t Equation (6) can also be simplified to Dtp f (t) = Dtq f (t) = Fig. 1. Elastic coefficient and viscosity: (a) Hookinan spring, a = 0; (b) Newtonian dashpot, a = 1; (c) Abel sticky pot, 0 < a < 1. When a = 1 shown as Fig. 1 (b), the behavior obeys the laws of Newtonian fluid, and the constitutive equation is given by (2). σ(t) = τ 1 ηDt1 ε(t), (2) where τ 1 = η/E is the relaxation time for the dashpot, η is dynamic viscosity, and E represents the elasticity modulus of the dashpot, and Dt1 ε(t) is the first time derivative of strain with respect to time t. In practical application of some materials or devices, the fluid behaves viscoelastically and the mechanical properties exhibit both spring and dashpot characteristics. We can use Fig. 1 (c) to describe the Abel sticky pot. σ(t) = τ α EDtα ε(t), (3) Dtα ε(t) where is the fractional derivative of α order of the strain with respect to time with evidently 0 ≤ α ≤ 1. B. Definition of Fractional Derivative The most common definition of Riemann-Liouville (R-L) fractional integral is given by[10] Z t dn 1 q (t − ξ)(n−q)−1 f (ξ)dξ, a0 Dt f (t) = Γ(n − q) dxn a0 n − 1 ≤ q < n, (4) dq . dtq (6) (7) C. Fractional Order Model of MR Damper In Fig. 2, the shock absorber is composed of an MR damper and a rubber damper, which possesses the advantages of the rubber and the MR damper. The damping force can be adjusted rapidly with little control energy requirement. In view of the structural characteristics of the coupled shock absorber, a typical standard linear solid model is presented, also called the Zener model[10] , as shown in Fig. 2 (b). The fractional order Zener model can be obtained by replacing the traditional Newton dashpot with the Abel dashpot. The constitutive relation can be written as[10] : σ + τ α Dα σ(t) = E2 τ α Dα ε(t) + E1 ε(t), 0 ≤ α < 1, (8) where E1 is the relaxed modulus, E2 is the unrelaxed modulus shown in Fig. 2. When a sinusoidal pressure is applied, the storage modulus (E ′ ) and loss modulus (E ′′ ) can be got from (8) ³ π´ + E1 E2 (ωτ )2α + (ωτ )α (E1 + E2 ) cos α 2 E′ = h ³ π ´i2 , (9) ³ π ´i2 h + (ωτ )α sin α 1 + (ωτ )α cos α 2 2 ³ π´ (E2 − E1 )(ωτ )α sin α ′′ 2 E =h ³ π ´i2 , (10) ³ π ´i2 h α α + (ωτ ) sin α 1 + (ωτ ) cos α 2 2 where in both (9) and (10), ω is angular frequency (rad/s), ω = 2πf , and f is frequency (Hz). where Γ(·) is gamma function, q is a non-integer order, a0 is the iterative initial value. In addition, the Caputo definition is often adopted in engineering applications, given by the following equation: Z t 1 C q (t − ξ)(n−q)−1 f (n) (ξ)dξ, a Dt f (t) = Γ(n − q) a n − 1 < q ≤ n, (5) In order to distinguish Caputo definition from R-L fractional calculus definition, we decorate it with the additional apex C. The fractional calculus definitions given by R-L and Caputo are all defined in time domain as a function f (t). The Laplace transformation of the R-L definitions is related to the initial value of the fractional differential and fractional calculus. Although the solutions can be found, a reasonable physical interpretation to these solutions is difficult to understand[6] . Fig. 2. The principle of the damper and the simplified model: (a) The schematic diagram of the shock absorber; (b) Simplified model; (c) The shock absorber. 290 IEEE/CAA JOURNAL OF AUTOMATICA SINICA, VOL. 3, NO. 3, JULY 2016 Substituting the structure parameters obtained from the coupled shock absorber into the (9) and (10), the storage modulus (E ′ ) and loss modulus (E ′′ ) can be calculated, shown in Figs. 3 and 4. The storage modulus E ′ increases as a function of the system frequency, while the loss modulus E ′′ is nonlinear. Using verified parameters and changing the order of the Abel dashpot from 0.2 to 1, E ′ and E ′′ exhibit different characteristics. When the frequency is less than 10 Hz, E ′ decreases with the increase of the order α. When the frequency is larger than 10 Hz, the law is opposite, showing that the smaller α produces larger elastic properties of the shock absorber. When 0 < α < 1, the E ′′ increases with the increase of the order α, showing an increase in the viscous behavior. When α = 1, the E ′′ has a fast drop when the frequency is larger than 10 Hz, when the frequency increases beyond a certain value, the loss modulus is smaller than a small α, as the Fig. 4 showing, the loss modulus in α = 1 is smaller than α = 0.8 when the frequency is larger than 32 Hz. Fig. 3. storage modulus is plotted as a function of frequency for two different values of α. In Fig. 5 (a), as I is increased, the storage modulus asymptote increases. Also as I is increased, the storage modulus approaches the asymptote more rapidly. In Fig. 5 (b) the asymptotic value of each current is larger than the corresponding current in Fig. 5 (a). Similar to Fig. 6 (a), the storage modulus gradually approaches the asymptote as frequency is increased. As I is increased, the storage modulus approaches the asymptote much more rapidly. The storage modulus of Zener model E ′ . Fig. 5. The storage modulus E ′ of the Zener model with different currents. Fig. 4. The loss modulus of Zener model E ′′ . When the applied magnetic field is manipulated according to the damping part of the coupled MR damper, the viscoelastic properties of the entire shock absorber can be changed dramatically. The magnetic field can be manipulated by changing the current, I, of the system. In Fig. 5, the In Fig. 6, the loss modulus is plotted as a function of frequency. In Fig. 6 (a), for I = 0, the loss modulus reaches a maximum for f = 8. As frequency is increased, the loss modulus gradually decreases. For I > 0, the modulus rapidly increases and reaches a maximum value for f = 4. The maximums for all values of I are all very similar in magnitude. However, as frequency is increased, larger currents possess a smaller loss modulus. In Fig. 6 (b), for I = 0, the maximum in loss modulus occurs at f = 8. The loss modulus then gradually decreases. For I > 0, the maximum occurs for f = 4. Furthermore, the decrease in loss modulus is much more rapid than what we observed in Fig. 6 (a) for α = 0.6. In addition, the maximum for all values of I is larger than the maximums observed in Fig. 6 (a). The viscous characteristics of the coupled MR damper are reflected in the low working CHEN et al.: FRACTIONAL MODELING AND ANALYSIS OF COUPLED MR DAMPING SYSTEM frequency. At large frequency the viscous performance of the damper is decreased, which is directly related to the working magnetic field. 291 p √ = k/m, critical damping coefficient cc = 2 km, and the damping factor µ = c/cc . So (11) can also be written as F sin ωt , (12) m and the two order vibration system in fractional order form can be given as: ẍ(t) + 2µωn2 ẋ(t) + ωn2 x(t) = D2 x(t) + 2µωn Dβ x(t) + ωn2 x(t) = P (t), 0 < β ≤ 1. (13) Fig. 7. The simplified model and the real experimental platform: (a) The simplified model of the experimental platform; (b) The real experimental platform. In order to simplify (13), here A1 = µωn , A2 = wn2 , so (13) can be written as follows: D2 x(t) + A1 Dβ x(t) + A2 x(t) = F (t). (14) Laplace transform was applied on the fractional differential (14) to get: s2 X(s) + A1 sβ X(s) + A2 X(s) = F (s). Fig. 6. The loss modulus E ′′ of the Zener model with different currents. III. E XPERIMENTAL P LATFORM From the above analysis, it can be seen that the fractional order model can accurately describe the viscoelasticity of the shock absorber. In order to further the study and analyze the dynamic performance of the damper, equivalent viscous damping is introduced[11] . The equivalent viscous damping is used to replace the complex damping machine. (15) The Caputo fractional derivative operator can also be used with initial values x(0+ ) = c0 , ẋ(0+ ) = c1 , this is called the composite fractional vibration equation. The transfer function for the fractional order system can be obtained by using the Laplace transform[12−13] : 1 G(s) = , 0 < β < 2. (16) s2 + µβ ωnβ sβ + ωn2 For the differential equation (11), the Grünwald-Letnikov (G-L) definition is used to solve the differential equation. The Grünwald-Letnikov method is the direct numerical method for solving fractional calculus. The G-L definition of fractional calculus is as follows: t−a A. Experiment Platform and Model Analysis From (8), the resistance, f (t), is provided, and the direction is opposite to the speed of the mass, m. The force applied to the system is F sin ωt, as shown in Fig. 7 (a), the two-order mode for a single degree of freedom dynamic system is defined as: mẍ(t) + cẋ(t) + kx(t) = F sin ωt, (11) where k is the stiffness of the damper, c is the damping coefficient. Here, some characteristic parameters of the vibration system are introduced: natural frequency of the system ωn h 1 X (β ) βi w i xt−jh a Dt x(t) = βi h j=0 j   t−a h X 1 (β ) = βi xt + wj i xt−jh  . h j=1 (17) In (17), a is the initial value for the numerical calculation, and (β ) to meet 0 < a < 1, h is the calculation step size, wj i is the coefficient of a polynomial (1 + z)βj , which can be derived from the following recursive formula[14] : µ ¶ βi + 1 (β ) (βi ) w0βi = 1, wj i = 1 − wj−1 , j = 1, 2, . . . . (18) j 292 IEEE/CAA JOURNAL OF AUTOMATICA SINICA, VOL. 3, NO. 3, JULY 2016 Equation (18) is substituted into (17), the numerical solution of (17) can be directly derived from differential equation:   t−a n h X X 1 a i (β ) P (t) − wj i xt−jh  , (19) xt = n βi X ai h j=1 i=0 i=0 hβi where xt is the sampled displacement data. P (t) is the controllable input external force, and ai denotes the iterative value during the process of calculation. The numerical solution and the analytical solution for the trinomial model (i.e., α2 = 2, α1 = α, α0 = 0) excited by unit step function are shown in Fig. 8. The solutions are calculated based on (15) and by the method of Adomian decomposition, respectively. In Fig. 8, the displacement is plotted as a function of time. From Fig. 8, both the numerical and analytical solutions exhibit similar displacement for all values of time considered. Therefore, the numerical solution of G-L can be applied to engineering analysis. In this paper, the fractional order model and the integer order model are analyzed using the numerical solution. instrument system, data acquisition terminal, software system using LabVIEW7.0 version, programmable current source, etc. IV. E XPERIMENT A NALYSIS The dynamic characteristics in the MR fluids are considered with changing mass percent of carbonyl iron. In Figs. 9 and 10, the mass percentages of carbonyl considered are 74 % and 78 %. The eccentricity is large for vibration frequencies of 10 Hz and 11 Hz. The dynamic parameters of the system model are described in Tables I and II. TABLE I f = 10 Hz, THE PARAMETERS A1 , A2 , β, D(e) OF THE MODEL IN DIFFERENT WORKING FLUIDS MRF (%) 74 I (A) 1 78 1 74 3 78 3 P A1 (s−2 ) A2 (s−1 ) β 28.625 28.723 0.6800 20.122 22.980 1 42.389 42.436 0.6890 36.234 38.452 39.015 1 195.236 58.963 58.967 0.8400 32.149 51.273 51.519 1 128.265 64.398 64.386 0.8410 34.572 59.581 59.815 1 168.426 D(e) 35.624 189.356 TABLE II f = 11 Hz, THE PARAMETERS A1 , A2 , β, D(e) OF THE MODEL UNDER DIFFERENT WORKING FLUIDS MRF (%) 74 Fig. 8. Adomian decomposition method solution vs G-L definition numerical solution. B. Measurement and Control Device The measurement and control system of MR damper is shown in Fig. 7. In the MR damper, the sensor collects the signals of vibration, displacement and acceleration. LabVIEW is used to process, analyze, and display the collected data. Then, based on the specific vibration control requirements and other related parameters (system structure, magnetorheological material characteristics, etc.), the required control current is calculated using the GBIP mode in LabVIEW. Through LabVIEW, the vibrational damping force can be controlled. The changes to the vibrational parameters of the experimental platform can be observed, and the output current can be adjusted to achieve a more desirable vibrational damping effect. The response of the shock absorber is of the order of several tens of milliseconds. A high signal sampling rate is required in order to meet the required vibrational reduction. The system consists of temperature, acceleration, and displacement sensors, as well as a data acquisition card of virtual I (A) 1 78 1 74 3 78 3 P A1 (s−2 ) A2 (s−1 ) β 30.058 29.264 0.6800 23.612 24.532 1 62.424 61.426 0.6950 40.096 100.000 100.000 1 258.812 63.912 62.117 0.8400 35.012 100.000 100.000 1 244.707 69.046 66.084 0.8430 37.155 100.000 100.000 1 234.314 D(e) 38.605 249.437 1) The order β of fractional order model is related to the vibration damping performance of MR fluids. With the same control current and the nonmagnetic saturation situation, the damping capacity and the model order β increase with the increase of the mass fraction of the carbonyl iron powder. As shown in Fig. 9, at f = 10 Hz, I = 1 A, the fractional order β increases from 0.68 to 0.689 as mass fraction M increased from 74 % to 78 % accordingly. Similar situation can be seen from Fig. 10, at f = 11 Hz, I = 1 A, the fractional order β increases from 0.68 to 0.695 as mass fraction M adjusted from 74 % to 78 %. From the results we can find that the order number is changed with the different working fluids. 2) The viscoelastic characteristics of the system with higher iron content is stronger: such as, f = 10 Hz, I = 1 A, when M = 74 %, 78 % respectively, the viscosity coefficients of A1 are 28.625, 42.389, which shows a significant increase; viscoelastic ratio ξ is respectively 0.9965 and 0.9988, which is also increased weakly, so can also be viewed unchanged. CHEN et al.: FRACTIONAL MODELING AND ANALYSIS OF COUPLED MR DAMPING SYSTEM 293 1.148, indicating that different MR liquids of the system have certain influence on the system. TABLE III I = 1 A, THE VARIANCE OF THE SAMPLED DATA SEGMENTS WITH DIFFERENT CURRENTS Fig. 9. f = 10 Hz, the displacements and the fractional order in different working fluids. 74 %, 78 %, σ12 (mm2 ) σ22 (mm2 ) 6000-6500 0.0313 0.0270 1.159 6500-7000 0.0316 0.0267 1.184 3 7000-7500 0.0313 0.0268 1.168 4 7500-8000 0.0313 0.0274 1.142 5 8000-8500 0.0308 0.0277 1.112 6 8500-9000 0.0306 0.0270 1.133 7 9000-9500 0.0307 0.0269 No. Sampled data 1 2 Average value σ12 /σ22 1.141 8.039/7 = 1.148 V. C ONCLUSIONS Fig. 10. f = 11 Hz, the vibration displacements and the fractional order in different working fluids. 3) ∆x1 and ∆x2 represent the variation displacement of the two working fluids in 1 A, 3 A respectively, it can be seen that ∆x1 > ∆x2 , and under the working current of 3 A, the changes of MR fluids damping characteristics are reducing. 4) The displacement of the theoretical fractional model and the integer order model are given in Figs. 9 and 10. It can be seen that the fitting curve of fractional order model is more close to the sampled displacement curve than that of the integer order model, and the results are in agreement with the computed results of Tables I and II. BasedPon the same sampled signal, the residual sum of squares D(e) obtained by fitting the fractional order models is less than that of the integer order models obtained by fitting the integer model, indicating that the fractional order system model is more accurate than the integer order system model. The effect of working fluids on the vibrational energy of the system is analyzed quantitatively by using the variance analysis. As shown in Table III, taking I = 1 A in Fig. 9 for example, σ12 is the variance at M = 74 %, and σ22 is variance at M = 78 %, σ12 /σ22 denotes the energy coefficient, we can find that the replacement of the working fluids has great influence on the dynamic energy coefficient, whose average value is The above analysis shows the mechanical properties of the coupled MR damper using viscous and elastic characteristics, presenting the properties of an elastic solid and a viscous fluid, and through the experiment, we have shown that: 1) the constitutive equation with fractional derivative method is derived from a strict formula, which has definite physical meaning; 2) the viscoelastic constitutive equation with the fractional derivative can be used to describe the mechanical vibration performance of the coupled MR damper with great accuracy than the integer order model; 3) the dynamic characteristics of the system are related to the order number of the fractional order model: under the same operating frequency, with the increase of the control current, the order of the fractional model is increased, and the viscoelastic properties of the shock absorber are enhanced. R EFERENCES [1] Stanway R, Sproston J L, Stevens N G. Non-linear modelling of an electro-rheological vibration damper. Journal of Electrostatics, 1987, 20(2): 167−184 [2] Spencer B F Jr, Dyke S J, Sain M K, Carlson J D. Phenomenological model for magnetorheological damper. Journal of Engineering Mechanics, 1997, 123(3): 230−238 [3] Zhou Qiang, Qu Wei-Lian. Two mechanic models for magnetorheological damper and corresponding test verification. Earthquake Engineering and Engineering Vibration, 2002, 22(4): 144−150 (in Chinese) [4] Gamota D R, Filisko F E. Dynamic mechanical studies of electrorheological materials: moderate frequencies. Journal of Rheology, 1991, 35(3): 399−425 [5] Slonimsky G L. Laws of mechanical relaxation processes in polymers. Journal of Polymer Science Part C: Polymer Symposia , 1967, 16(3): 1667−1672 [6] Friedrich C. Relaxation and retardation functions of the maxwell model with fractional derivatives. Rheologica Acta, 1991, 30(2): 151−158 [7] Bagley R L, Torvik P J. On the fractional calculus model of viscoelastic behavior. Journal of Rheology, 1986, 30(1): 133−155 294 IEEE/CAA JOURNAL OF AUTOMATICA SINICA, VOL. 3, NO. 3, JULY 2016 [8] Paggi M, Sapora A. An accurate thermoviscoelastic rheological model for ethylene vinyl acetate based on fractional calculus. International Journal of Photoenergy, 2015, 2015: Article ID 252740 [9] Jóźwiak B, Orczykowska M, Dziubiński M. Fractional generalizations of maxwell and kelvin-voigt models for biopolymer characterization. PLoS One, 2015, 10(11): e0143090 Chunyu Li Engineer at Fujian University of Technology, Fuzhou, China. She received the B. Sc. degree in mechanical engineering from Beihua University, Jilin, China, in 2004. Her research interests include smart materials and laser cladding. [10] Alcoutlabi M, Martinez-Vega J J. Application of fractional calculus to viscoelastic behaviour modelling and to the physical ageing phenomenon in glassy amorphous polymers. Polymer, 1998, 39(25): 6269−6277 [11] Li Zhuo, Xu Bing-Ye. Equivalent viscous damping system for viscoelastic fractional derivative model. Journal of Tsinghua University (Science and Technology), 2000, 40(11): 27−29 (in Chinese) [12] Wang Zhen-Bin, Cao Guang-Yi, Zhu Xin-Jian. Application of fractional calculus in system modeling. Journal of Shanghai Jiaotong University, 2004, 38(5): 802−805 (in Chinese) [13] Wang Zhen-Bin, Cao Guang-Yi. Two system modeling methods using fractional calculus. Journal of System Simulation, 2004, 16(4): 810−812 (in Chinese) Benjamin Wilson Research assistant at University of Wisconsin-Madison. He received the B. S. degree from Purdue University in Chemical Engineering, USA. He received the Ph. D. degree in January 2016, also from Purdue University in Chemical Engineering, USA. His research interests include magnetorheological (MR) fluids. [14] Xue D Y, Chen Y Q. Solving Applied Mathematical Problems with Matlab. Boca Raton: CRC Press, 2008. Bingsan Chen Associate professor at Fujian University of Technology, Fuzhou, China. He received the B. Sc. degree in mechanical engineering from Beihua University, Jilin, China, in 2004. He received the Ph. D. degree in mechanical manufacture from Huaqiao University, Xiamen, China, in 2009. His research interests include smart materials and signal processing of the mechanical system. Corresponding author of this paper. Yijian Huang Professor at Huaqiao University, Xiamen, China. He received the B. Sc. degree in physics from Xiamen University, Xiamen, China, in 1968. He received the M. Sc. degree in fluid drive and control from Zhejiang University, Hangzhou, China, in 1981. His research interests include smart materials and signal processing of the mechanical system.