*Corr. Author’s Address: Department of Mechanical Engineering, National Defence University of Malaysia, Malaysia, k.hudha@upnm.edu.my 471
Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482 Received for review: 2023-04-14
© 2023 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2023-09-18
DOI:10.5545/sv-jme.2023.604 Original Scientific Paper Accepted for publication: 2023-09-22
Impact Behaviour Modelling of Magnetorheological Elastomer 
Using a Non-Parametric Polynomial Model Optimized  
with Gravitational Search Algorithm
Zubir, A.R.
 
– Hudha, K. – Abd. Kadir, Z. – Amer, N.H.
Amrina Rasyada Zubir
 
– Khisbullah Hudha
*
 – Zulkiffli Abd. Kadir – Noor Hafizah Amer
Department of Mechanical Engineering, National Defence University of Malaysia, Malaysia
This paper presents an approach to model the impact behaviour of a dual-acting magnetorheological elastomer (MRE) damper using 4
th
-
order polynomial functions optimized with a gravitational search algorithm. MRE is a type of smart material that can change its mechanical 
properties in response to an injected current, making it well-suited for a wide range of applications such as vibration absorption, noise 
cancellation, and shock mitigation. The proposed model uses a combination of polynomial functions designed to predict the nonlinearity of 
MRE during compression and extension stages. The model is tuned and validated using experimental data from impact tests conducted on 
the MRE damper under various currents. The results indicate that the developed model can accurately track the impact behaviour of MRE 
with minimum error. Additionally, an interpolation model is proposed to estimate the appropriate forces for median currents. The interpolation 
model predicts the force between the upper and lower currents, demonstrating the model's ability to predict MRE behaviour accurately. The 
main contribution of this study is proposing a non-parametric model of MRE that is able to identify the hysteretic behaviour of the MRE based 
on specific current applied. In addition, an interpolation model is introduced in this study to cover not only the input current starting from 0 A 
to 2 A but also the intermediate current such as 0.3 A, 0.7 A, 1.3 A and 1.7 A.
Keywords: magnetorheological damper, polynomial model, gravitational search algorithm, force-displacement characteristic, 
interpolated model
Highlights
• 	 The use of magnetorheological elastomer in the automotive field specific for impact mitigation is one of the precautionary 
steps to reduce frontal collision towards the vehicle body and passenger. 
• 	 This study aims to investigate the behaviour of the magnetorheological elastomer damper using a polynomial model optimized 
with a gravitational search algorithm and compare it with experimental results from the drop impact test. 
• 	 The polynomic model that can perform at interpolation current has been developed, and the simulation results are well in 
agreement with the experimental results.
0  INTRODUCTION
The consumption of magnetorheological elastomer 
(MRE) has grown widely in the automotive field 
due to its potential to cope with variable stiffness 
and damping subjected to magnetic fields [1] and [2]. 
MREs are used in absorbers, isolator devices, and 
suspension systems to reduce unwanted vibration 
transmitted to the occupants [3] to [5]. However, 
during medium and hard collisions in which the 
vehicle’s velocity exceeds 10 mph, MREs cannot 
protect the vehicle structure and the passenger. To 
address this problem, researchers have initiated 
modelling the impact behaviour of the vehicle using 
non-parametric approaches and have developed a 
dual-acting MRE damper for the study of impact 
absorption characteristics.
Previous works on the design of MRE dampers 
have proposed adaptive variable stiffness absorbers 
made of two MRE parts attached to both piston 
and outer cylinder, a multi-layered MRE isolator 
device, and a semi-active system that incorporates 
both passive and active elements in the design [6] 
and [7]. Aside from creating a practical design of 
MRE-based devices, it is also critical to model and 
simulate the MRE’s dynamic behaviour under various 
excitation circumstances. In the modelling of MRE-
based devices, there are two common approaches: 
parametric and non-parametric. Non-parametric 
models are more favourable since they do not need 
any model parameters and only require input and 
output to operate [8]. Non-parametric modelling 
uses intelligent paradigms such as adaptive neuro-
fuzzy inference systems (ANFIS), which are capable 
of modelling MRE damper behaviour under impact 
loading [9] and [10].
Yu et al. [11] introduced a versatile model 
capable of predicting dynamic behaviour under 
diverse excitation conditions, specifically tailored 
for applications in structural seismic mitigation. 
This model employs the Kelvin-Voigt element to 
elucidate the viscoelastic behaviour of the device 
and leverages the Bouc-Wen hysteresis element to 
delineate the strain-stiffening phenomenon in its 

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
472 Zubir, A.R. – Hudha, K. – Abd. Kadir, Z. – Amer, N.H.
response. The model’s parameters are determined 
through an improved particle swarm optimization 
(IPSO) approach aimed at minimizing deviations 
between observed values and model predictions. 
Subsequently, the field-dependent properties of the 
isolator, influenced by varying current levels, are 
validated using both random and earthquake-induced 
test data.
A new hybrid model employing a curve-fitting 
method has also been developed to accurately depict 
the highly nonlinear and hysteresis relationships 
between shear force and displacement responses in 
magnetorheological (MR) elastomer-based isolators 
[12]. The hybrid content in this model involved 
the combination of the hyperbolic sine function 
and Gaussian function to capture the hysteresis 
loops exhibited by the device responses effectively. 
Furthermore, an enhanced fruit fly optimization 
algorithm (FOA) has been introduced for optimizing 
the model parameters. This improved FOA 
incorporates a self-adaptive step size mechanism, 
rather than a fixed step, to strike a balance between 
the global and local optimum search capabilities 
of the algorithm. To verify the performance of this 
innovative hybrid model, both harmonic and random 
excitations have been employed in comprehensive 
testing.
This study contributes to the establishment of 
an effective modelling concept that can analyse 
the impact behaviour of the MRE damper based 
on a non-parametric approach using a 4
th
-order 
polynomial model. The model is designed to suit 
various currents ranging from 0 A to 2 A, as well as 
in-between currents using an interpolation method 
for the currents of 0.3 A, 0.7 A, 1.3 A and 1.7 A. The 
effectiveness of the developed model is validated 
through a model verification procedure by comparing 
the results between the developed model in Matlab-
Simulink software and the experimental results 
obtained through the drop impact test on the MRE 
damper. The level of agreement for both responses is 
observed through the force prediction error, validating 
the effectiveness of the developed model in predicting 
the hysteretic behaviour of the MRE [13].
This paper is organized as follows: Section 1 
discusses several previous works on the design of 
an MRE damper; Section 2 describes the design of 
the dual-acting MRE damper; Section 3 explains 
the process of fabrication for MRE and the force-
displacement characteristics; Section 4 describes the 
modelling method of MRE behaviour under impact 
loading using 4
th
-order polynomial model; Section 
5 shows the simulation results and the validation of 
the model with experimental data. The last section 
consists of some discussion and conclusion of this 
study.
1  DESIGN OF DUAL-ACTING MAGNETORHEOLOGICAL 
ELASTOMER DAMPER FOR IMPACT MITIGATION
In this study, a conceptual design for a dual-acting 
MRE damper is presented. The damper is composed 
of an upper impactor, a piston, MREs, a coil, and a 
polyoxymethylene (POM) bobbin, which serves as 
the isolator of the coil to the housing. To enhance 
the dual-acting absorption capability of the damper 
during compression and extension stages, three units 
of MREs are installed vertically along the housing for 
both upper and lower sections, as illustrated in Fig. 1. 
The MREs exhibit variable damping properties and 
stiffness under changeable applied magnetic fields, 
allowing the MRE damper to return to its original 
position after impact [14].
Fig. 1.  Design of dual-acting MRE damper
The damper consists of six main components: the 
impactor, piston, housing, cap, POM bobbin, and the 
MRE itself, all made of mild steel, as shown in Fig. 2. 
This damper is designed to adapt to impacts in both 
upward and downward directions. To induce current 
flow to the MRE, a 0.7 mm copper coil is wrapped 
around the POM bobbin, which acts as an isolator unit 
between the MRE housing and the coil. This current 
flow strengthens the bonding between the ferrous 
particles and the rubber matrix in the MRE [15].
2  EXPERIMENTAL SETUP OF DUAL-ACTING 
MAGNETORHEOLOGICAL ELASTOMER DAMPER
This section outlines the experimental procedure for 
the fabrication of MRE samples and the subsequent 
testing using a drop impact machine to obtain force-

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
473 Impact Behaviour Modelling of Magnetorheological Elastomer Using a Non-Parametric Polynomial Model Optimized with Gravitational Search Algorithm 
displacement characteristics. The fabrication process 
was conducted at the Automotive Laboratory of 
Universiti Pertahanan Nasional Malaysia (UPNM), 
with the mixing and curing processes carried out 
manually using an isotropic approach without the 
presence of a magnetic field [16]. Firstly, the materials 
needed in the fabrication process of MRE, namely 5 
µm of carbonyl iron particles (CIP), room temperature 
vulcanization (RTV) silicone rubber, carbon black and 
hardening agent, were measured using a digital scale. 
The composition of each material was set according to 
the optimal composition proposed by [17] and shown 
in Table 1.
Table 1.  Share of MRE samples
Materials Share of MRE samples [%] Weight [g]
RTV silicone rubber 30 9
Carbonyl iron powder 60 18.9
Additive 10 of CIP weight used 1.89
Hardener 3 of the RTV weight used 0.27
The materials mentioned above were mixed 
in a container and stirred until a homogeneous 
combination was obtained. The mixture was then 
poured into two different types of moulds, solid and 
ring moulds and left to cure for approximately 24 
hours. Before pouring the mixture into the moulds, a 
releasing agent spray was applied to ease the removal 
of the samples from the mould [18]. Finally, the MRE 
samples were ready for the drop impact test. To 
prepare the specimen for the test, three solid MREs 
were stacked at the bottom of the damper, and another 
three ring MREs were placed along the piston. The 
overall specimen preparation process is illustrated in 
Fig. 3.
Upon completing the fabrication of MRE 
samples, the experimental test continues using the 
Instron Drop Impact Machine to conduct a drop 
impact test. The primary objective of this test is to 
assess the impact force of the MRE with the additive 
material added, which is expected to improve its 
impact absorption capabilities. To start the test, three 
Fig. 3.  MRE fabrication process
Fig. 2.  The prototype of a dual-acting MRE damper

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
474 Zubir, A.R. – Hudha, K. – Abd. Kadir, Z. – Amer, N.H.
units of MRE are loaded into the dual-acting damper 
in parallel positions for both upper and lower sections. 
The impactor is then placed at the top of the damper, 
as shown in Fig. 4. The striker holder is released 
from a pre-determined position set according to the 
parameters in the CEAST software, as detailed in 
Table 2. The impact velocity is set to 2.24 m/s, which 
represents the impact energy for a light impact.
Fig. 4.  Experiment setup of MRE drop impact
Table 2.  Experimental parameters set for impact test
Parameter Unit Value
Impact energy [J] 13.80
Impact velocity [m/s] 2.24
Falling height [mm] 256.00
Impact point offset [mm] 36.00
Tup holder mass [kg] 4.30
Tup nominal mass [kg] 1.20
Total mass [kg] 5.50
Current [A] 0 to 2
The hydraulic system of the impact machine is 
equipped with an accelerometer sensor and load cell to 
record the actual acceleration and impact force during 
testing. The recorded behaviours of the damper during 
impact loading by the sensors are then sent to the data 
acquisition system (DAQ). The impact absorption 
is analysed by observing the force-displacement 
characteristics, which will be used as the benchmark 
for MRE modelling in the next section. 
3 MODELLING AND OPTIMIZATION OF 
MAGNETORHEOLOGICAL ELASTOMER UNDER IMPACT 
LOADING BASED ON POLYNOMIAL MODEL
A non-parametric polynomial model was developed to 
analyse the dynamic behaviour of the MRE element 
under impact loading due to varying currents ranging 
from 0 A to 2 A. Fig. 5 illustrates the typical F-d 
characteristic of MRE. The blue line indicates the 
compression state of the MRE, where the force is 
exerted increasingly until maximum with increasing 
displacement during impact. Meanwhile, the red 
line indicates the extension state of the MRE, where 
the force decreases as the displacement of the MRE 
decreases to the original state.
Fig. 5.  Force-displacement characteristic
The polynomial function can be generated by 
considering a fourth-order polynomial function. In this 
function, y
1
 represents the function of the ascending 
line known as compression force. Meanwhile, y
2
 
represents the function of the descending line known 
as the extension force, as illustrated in Fig. 5. The 
functions of both compression and extension stages 
are shown in Eqs. (1) and (2):
	𝑦
1
 = 𝑎 y1–1
𝑥 4
 + 𝑎 y1–2
𝑥 3
 + 𝑎 y1–3
𝑥 2
 + 𝑎 y1–4
𝑥 + 𝑎 y1–5 , 
(1) 
 𝑦 2
 = 𝑎 y2–1
𝑥 4
 + 𝑎 y2–2
𝑥 3
 + 𝑎 y2–3
𝑥 2
 + 𝑎 y1–2
𝑥 + 𝑎 2–5 
. (2)
Here, 𝑎 yi–j
 is a real number, and 𝑥 is the 
displacement in meters obtained from the 
experimental data. The values of the constants depend 
on the curve of the line and can be tuned to obtain 
better force-displacement characteristics. In this study, 
the modelling of a non-parametric polynomial model 
requires the displacement, velocity, impact energy, 
and current to produce the output force, as explained 

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
475 Impact Behaviour Modelling of Magnetorheological Elastomer Using a Non-Parametric Polynomial Model Optimized with Gravitational Search Algorithm 
in Fig. 6. The unknown values 𝑎 yi–j 
represent the 
magnitude, following the trend, and connection factor, 
respectively. The optimized values for the model 
parameters 𝑎 yi–j
 are required to produce an ideal shape 
of the hysteresis loops for the yielding element [19].
Fig. 6.  MRE non-parametric polynomial model
3.1  Optimization of MRE Model with Gravitational Search 
Algorithm
To obtain the optimized parameters, a gravitational 
search algorithm (GSA) was implemented in the 
simulation model as it can identify the parameters that 
minimize a specific goal or fitness function. The flow 
chart for the model identification procedure based on 
GSA is presented in Fig. 7. 
Fig. 7.  Model identification procedure based on GSA
In GSA coding; the system incorporates a set of 
agents called masses, gravitational law, and Newton’s 
law of motion to achieve the optimal result with a 
flexible implementation concept. GSA employs a 
similar approach to PSO, whereby data optimization 
is achieved through the exploration and exploitation 
abilities of the agents in the search space [20]. The 
agents in GSA are analogous to particles in the 
universe, and their locations can be represented 
using Eq. (3). The entire system of N agents can be 
expressed as shown in Eq. (4).
 X
i
 = (x
i
1
, x
i
2
, …, x
i
d
, …, x
i
n
), (3)
 X = (X
1
, X
2
, …, X
i
, …, X
N
). (4)
In the preceding equations, x
i
d
 represents the 
location of agent i in dimension d, n represents the 
dimension of the search space and N represents the 
number of individuals (or agents). Eq. (5) is used to 
calculate the applied force on agent i by agent j at a 
given time t. M
aj
 denotes agent j's active gravitational 
mass, M
pi
 denotes agent i's passive gravitational mass, 
G(t) denotes the gravitational constant at time t, ɛ is a 
small constant, and R
ij
 denotes the Euclidian distance 
between the two agents.
 Ft Gt
Mt Mt
Rt
xt xt
ij
d aj pi
ij
j
d
i
d
() ()
() ()
()
(( )( )). 




 (5)
Gravitational fixed G is a time-dependent 
function that begins with the initial value G
0
 and 
decreases over time to control the search accuracy. 
Eqs. (6) and (7) are used to calculate the value of this 
function.
 G ( t ) = G ( G
0
, t ), (6)
 Gt Ge
t
T
() , 

0

 (7)
where α and G
0
 are constants and T indicates all 
iterations. Eqs. (8) to (10) are used to update inertia 
and gravitational masses in Eq (5).
 (M
ai
 = M
pi
 = M
ii
 = M
i
),    i = 1, 2, ..., N,  (8)
 mt
fit tw orst t
best tw orst t
i
i
 
  
  
, (9)
 Mt
mt
mt
i
i
j
N
j





1
, (10)
where fit
i
(t) shows the fitness value of agent i at time 
t, and worst(t) and best(t) are calculated using Eq. (11) 
and (12). 
 best tm in fit t
jN j
 
 12 ,, ...,
, (11)
 worst tm ax t
jN
 
 12 ,, ...,
. (12)
Here, best(t) is the smallest fitness value among 
all agents, while worst(t) is the biggest fitness value 
among all agents. These two important values are used 
in a minimisation problem where the fitness value, 

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
476 Zubir, A.R. – Hudha, K. – Abd. Kadir, Z. – Amer, N.H.
fit
j
(t), influences the mass value of i
th
 agents, which 
correlates to the particle’s position in a search space. 
The method then examines each agent's fitness, fit (t) 
in relation to the goal function.
3.2  Simulation Results of Optimized MRE Model
This study utilizes the GSA method to determine 
the optimal values for 𝑎 yi-j
 in both the extension and 
compression stages of a 4
th
-order polynomial model 
to create the hysteresis loop. The objective is to 
optimize the polynomial parameters by comparing the 
simulation results with experimental responses and 
adjusting the performance index accordingly. Table 3 
displays the initial simulation parameters. The model 
configuration parameters were set in Matlab-Simulink 
with a fixed-step size of 0.01 and discrete time (no 
continuous states). The model was run for 36 s.
Table 3.  Initial simulation parameters for GSA
Parameters Values
Number of dimensions, D 5
Number of agents, N 50
Number of iterations, T 60
Range of  
dimension
Compression 
case:
–1100 ≤ a
y1–1
 ≤ 4000
–9000 ≤ a
y1–2
 ≤ 20000
–25000 ≤ a
y1–3
 ≤ 3000
4000 ≤ a
y1–4
 ≤ 20000
–300 ≤ a
y1–5
 ≤ 400
Extension 
case:
–9000 ≤ a
y2–1
 ≤ 33000
–62000 ≤ a
y2–2
 ≤ 21000
–8000 ≤ a
y2–3
 ≤ 35000
500 ≤ a
y2–4
 ≤ 800
–70 ≤ a
y2–5
 ≤ 3
During the initialization step, the algorithm 
generates a random number and assigns it to a 
variable. Next, a random number of agents, N, 
is assigned to each entity in the population. The 
simulation evaluates the fitness function for each 
entity and performs random value selection processes 
[21]. Each agent’s fitness is then evaluated for the next 
iteration, and the constraints are checked. This process 
is iterated over multiple generations to achieve the 
best results. Finally, the optimized polynomial model 
parameters for MRE are tabulated in Table 4.
The optimized parameter values for the simulation 
model are determined by fitting the model to the 
experimental data obtained from the drop impact test. 
The responses of MRE in terms of force-displacement 
characteristics are observed for both simulation and 
experiment, from 0 A to 2 A of input currents with 
0.5 A increments. As shown in Fig. 8, each impact 
energy produced by the model can be used to track the 
energy recorded in the experimental test by forming 
hysteresis loops for each applied current. To evaluate 
the accuracy of this model, the differences between 
the experimental and simulation are examined based 
on the highest error. The highest errors are marked in 
the red region in Fig. 8 for all cases. The percentage 
errors are then calculated based on the difference 
in force produced as tabulated in Table 5. From the 
overall results, the percentage error obtained is less 
than 13 % where within the acceptable range of error 
[22].
Table 4.  Optimized model parameter for each input current
Input 
current [A]
Model parameter
Compression Extension
0 A
a
y1–1
–1201.40 a
y2–1
–2198.30
a
y1–2
4792.40 a
y2–2
5862.70
a
y1–3
–8317.70 a
y2–3
–1052.60
a
y1–4
8945.10 a
y2–4
465.53
a
y1–5
–236.86 a
y2–5
2.19
0.5 A
a
y1–1
–10630.00 a
y2–1
–8708.80
a
y1–2
15621.00 a
y2–2
20027.00
a
y1–3
–3332.90 a
y2–3
–7723.00
a
y1–4
3968.60 a
y2–4
1209.50
a
y1–5
93.28 a
y2–5
–19.39
1 A
a
y1–1
3422.00 a
y2–1
–2193.60
a
y1–2
–8986.60 a
y2–2
8677.00
a
y1–3
2936.20 a
y2–3
–7363.90
a
y1–4
9330.20 a
y2–4
7180.70
a
y1–5
336.52 a
y2–5
–22.37
1.5 A
a
y1–1
–7302.70 a
y2–1
–1962.30
a
y1–2
8103.30 a
y2–2
7596.90
a
y1–3
–1859.90 a
y2–3
–2696.10
a
y1–4
8441.90 a
y2–4
4795.50
a
y1–5
37.99 a
y2–5
–65.56
2 A
a
y1–1
–6619.60 a
y2–1
32293.00
a
y1–2
19772.00 a
y2–2
–61171.00
a
y1–3
–24021.00 a
y2–3
35010.00
a
y1–4
19016.00 a
y2–4
3166.00
a
y1–5
–105.53 a
y2–5
–58.19
Table 5.  Force prediction error for input current
Current 
[A]
Experimental 
data
Simulation 
data
Percentage 
error [%]
0 2079.22 1821.72 12.38
0.5 5458.78 5810.15 6.44
1.0 6603.41 6097.63 7.66
1.5 7231.52 6501.59 10.09
2 6421.51 6035.60 6.01

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
477 Impact Behaviour Modelling of Magnetorheological Elastomer Using a Non-Parametric Polynomial Model Optimized with Gravitational Search Algorithm 
Fig. 9 shows the superimposed force-
displacement model responses for all currents. 
The trend indicates that the MRE force increases 
proportionally with the input current while the MRE 
displacement decreases proportionally with the input 
current. These simulation results demonstrate that 
the MRE model can provide variable stiffness for an 
active impact device. The hysteresis loop results also 
show a consistent increment of impact force, with a 
maximum value of 7600 N for both 1.5 A and 2 A. 
a)      b) 
c)     d) 
e) 
Fig. 8.  Force-displacement characteristics of the MRE damper at different currents: a) 0 A, b) 0.5 A, c) 1 A, d) 1.5 A, and e) 2 A

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
478 Zubir, A.R. – Hudha, K. – Abd. Kadir, Z. – Amer, N.H.
of 0.3 A, 0.7 A, 1.3 A, and 1.7 A, as shown in Fig. 10. 
The output force of the subsystems for each current 
was combined to generate a force-displacement curve 
for the developed interpolation model.
Fig. 10.  The structure of the interpolation current prediction 
model
4.1  Validation of the Dual-Acting MRE Behaviours for The 
Middle Current
The behaviours of a dual-acting magnetorheological 
(MRE) damper for middle currents have been tested 
by comparing simulation results obtained from 
Matlab with experimental results obtained from a 
drop impact machine. To begin, drop impact tests 
were carried out on the MRE damper by applying a 
set of input currents of 0.3 A, 0.7 A, 1.3 A, and 1.7 
A. The force-displacement characteristic responses 
obtained in simulation for each input current were 
then compared with the experimental data for 
validation purposes. Based on the results shown in 
Fig. 11, the model’s behaviours showed an increase in 
force with increasing input current despite a decrease 
in displacement. The trend of the input current at 0.3 
A and 0.7 A demonstrated the hysteresis behaviour of 
the MRE, which experienced compression and tension 
when an external force was applied. Similarly, for 
input currents of 1.3 A and 1.7 A, the model managed 
to follow the trend for the compression part, but the 
tension part for both currents diverged, and the trends 
were not as well-defined as those for the 0.3 A and 0.7 
A currents. 
Table 6.  Force prediction error for the interpolation current
Current 
[A]
Experimental 
data
Simulation data Percentage 
error [%]
0.3 5025.37 5187.61 3.23
0.7 6376.98 6253.14 1.94
1.3 6884.49 6270.71 8.92
1.7 7423.21 6240.26 15.94
To check the similarity between the developed 
model and the experimental results of the MRE, the 
percentage error is evaluated by undergoing a similar 
procedure as the input currents of 0 A to 2 A. The 
The insignificant difference in impact force for 2 A 
compared to 1.5 A might be due to the limitations of 
the MRE design in terms of its stroke, resulting in an 
almost similar magnitude of impact force. However, 
the reduced displacement characteristic of 2 A 
compared to 1.5 A is still within the acceptable range 
of impact force that can be used as an active impact 
device.
Fig. 9.  Comparison of force-displacement curve  
for simulation at varying current
4  INTERPOLATION MODEL FOR  
MEDIAN FORCES OF MRE DAMPER
The next contribution of this study is proposing a 
method for predicting median forces based on median 
applied currents using an interpolation approach. 
Interpolation is a curve-fitting method that uses linear 
polynomials to create new data points within a given 
range. For instance, if the range of the current set is 
from 0.5 A to 1 A, with known forces of f
0.5
 and f
1 
obtained in the previous section, the interpolation 
approach can predict the force produced by the MRE 
for a median current of 0.7 A. The predicted force for 
f
0.7
 can be calculated using the following equation:
 
CC
CC
ff
ff
07 05
10 5
07 05
10 5
..
.
..
.
.





 X
By inserting the values of C
0.5 
and C
1
 as well as 
data for both f
0.5
 and f
1
, the final equation to predict 
the hysteresis loop for f
0.7
 is written as follows:
 f
ff f
07
10 50 5
07 05
10 5
.
..
..
.
. 
    
 
 X
The interpolation model for MRE was executed in 
the Matlab-Simulink environment using input currents 

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
479 Impact Behaviour Modelling of Magnetorheological Elastomer Using a Non-Parametric Polynomial Model Optimized with Gravitational Search Algorithm 
difference between the simulation and experimental 
forces are marked as presented in Fig. 11, and their 
values are tabulated in Table 6. From the percentage 
of error, it can be seen that the error is less than 20 %, 
which is within an acceptable range of error [22].
4.2  Verification of the MRE Middle Current with Upper and 
Lower Responses
To check the acceptable range of interpolation models, 
the force-displacement characteristics were compared 
with upper and lower currents for each case. For 
instance, the interpolation model with a 0.3 A current 
it is compared with 0 A and 0.5 A. Based on the 
comparison results, it can be observed that the force 
for the compression stage is within the range of forces 
produced for 0 A to 0.5 A, which indicates that it falls 
within an acceptable range. Similarly, for the extension 
stage, the force produces a similar range between 0 A 
to 0.5 A currents, which is also within the acceptable 
range. Fig. 12 illustrates all the interpolation models 
for 0.7 A, 1.3 A, and 1.7 A currents, and the model 
successfully generates the in-between currents and 
validates the proposed interpolation model. 
4.3  Effects of Varying GSA Size and the Number of 
Iterations
The number of iterations plays a crucial role in GSA, 
as it sets the duration for which particles can search 
for the optimal solution. Generally, increasing the 
number of iterations can enhance the algorithm's 
capacity to discover the global optimum, as it provides 
the agents with more time to explore the search space 
and converge on the best solution. Nevertheless, 
after a certain point, additional iterations may not 
contribute to any further enhancement in the solution, 
as the agents may have already converged to a local 
a)      b) 
c)      d) 
Fig. 11.  Force-displacement curve of the interpolation current; a) 0.3 A, b) 0.7 A, c) 1.3 A, and d) 1.7 A

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
480 Zubir, A.R. – Hudha, K. – Abd. Kadir, Z. – Amer, N.H.
optimum. Conversely, decreasing the number of 
iterations can lead to faster execution times but at 
the risk of potentially missing out on better solutions. 
Thus, it is vital to strike a balance between the number 
of iterations and the desired level of performance index 
and execution time. Based on Fig. 13, 100 iterations 
were chosen as they have a quick convergence rate 
and can achieve a superior performance index.
Enlarging the agent size can enhance the 
exploration of the solution space on a global level 
since a greater number of agents are involved in the 
search. Nonetheless, it may escalate the computational 
expenses and impede the convergence rate because 
of an increased need for communication and updates 
among the agents. Based on the results depicted 
in Fig. 14, this study opted for a size of 100, which 
exhibited a swift convergence rate.
In order to check the validity of the proposed 
model, the simulation results for currents 0 A, 1 A and 
2 A are compared with the Multi-Layers Exponential 
based Preisach Model, the MRE modelling carried out 
by [23] as presented in Fig. 15.
15.
Fig. 13.  Effect of varying the number of iterations on the 
performance of GSA
a)      b) 
c)      d) 
Fig. 12.  Comparison of the force-displacement curve for simulation at varying currents;  
a) 0 A to 0.5 A, b) 0.5 A to 1 A, c) 1 A to 1.5 A, and d) 1.5 A to 2 A

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
481 Impact Behaviour Modelling of Magnetorheological Elastomer Using a Non-Parametric Polynomial Model Optimized with Gravitational Search Algorithm 
Fig. 14.  Effect of varying the number of particles  
on the performance of GSA
Fig. 15.  Comparison between Exponential model (EM)  
and Polynomial model (PM)
5  CONCLUSIONS
This paper presents a comprehensive study on the 
hysteresis behaviour modelling of magnetorheological 
elastomers under impact loadings. The study 
introduces a 4
th
-order polynomial model that is 
optimized with the gravitational search algorithm, 
resulting in a reliable and accurate framework that 
can capture the complex hysteresis behaviour of 
the material. The developed model is shown to 
perform exceptionally well in capturing the dynamic 
response of magnetorheological elastomers under 
various impact-loading scenarios. The model’s 
response closely matches the experimental data, 
with a maximum prediction error of less than 13 %. 
Additionally, the interpolated model’s response agrees 
well with the experimental data, with a maximum 
percentage error of less than 15.94 %. The study also 
investigates the impact of changing the number of 
iterations and the number of agents on the performance 
of the gravitational search algorithm. Overall, the 
results suggest that the proposed model is a promising 
approach for accurately predicting the hysteresis 
behaviour of magnetorheological elastomers under 
impact loads. For future study of this research, the 
developed MRE model will be used as an actuator 
for the active bumper system to absorb impact force 
during frontal collision. It is expected that the MRE 
actuator can reduce the impact force and thus reduce 
vehicle acceleration and jerks. 
6  ACKNOWLEDGEMENTS
The research has been carried out under the 
Fundamental Research Grant Scheme project 
FRGS/1/2021/TK02/UPNM/02/1 provided by the 
Ministry of Education of Malaysia.
7  REFERENCES
[1] Fu, J., Liu, J., Lai, J., Zhong, C., Dai, Z., Yu, M. (2023). 
Robustness analysis of magnetorheological elastomer-based 
vibration isolation system with optimal fuzzy controller. Smart 
Materials and Structures, vol. 32, no. 3, DOI:10.1088/1361-
665X/acb577.
[2] Ahamed, R., Rashid, M.M., Ferdaus, M. M., Yusof, H.M. (2016). 
Advancement in energy harvesting magneto-rheological fluid 
damper. Korea Australia Rheology Journal, vol. 28, p. 67-74, 
DOI:10.1007/s13367-016-0007-6.
[3] Masa’id, A., Lenggana, B.W., Ubaidillah, U., Susilo, D.D., 
Choi, S.-B. (2023). A review on vibration control strategies 
using magnetorheological materials actuators: Application 
perspective. Actuators, vol. 12, no. 3, art.ID. 113, 
DOI:10.3390/act12030113.
[4] Yu, Y., Li, J., Li, Y., Li, S., Li, H., Wang, W. (2019). Comparative 
investigation of phenomenological modeling for hysteresis 
responses of magnetorheological elastomer devices. 
International Journal of Molecular Sciences, vol. 20, no. 13, 
art. ID 3216, DOI:10.3390/ijms20133216.
[5] Choi, Y.T., Wereley, N.M. (2022). Adaptively tunable 
magnetorheological elastomer-based vibration absorber for 
a propeller aircraft seat. AIP Advances, vol. 12, no. 3, art. ID 
035332, DOI:10.1063/9.0000323.
[6] Zhao, L., Yu, M., Fu, J., Zhu, M., Li, B (2017). A miniature MRE 
isolator for lateral vibration suppression of bridge monitoring 
equipment: design and verification. Smart Materials and 
Structures, vol. 26, no. 4, art. ID 047001, DOI:10.1088/1361-
665X/aa5d97.
[7] Gu, X., Li, J., Li, Y., Askari, M. (2015). Frequency control of 
smart base isolation system employing a novel adaptive 
magneto-rheological elastomer base isolator. Journal of 

Strojniški vestnik - Journal of Mechanical Engineering 69(2023)11-12, 471-482
482 Zubir, A.R. – Hudha, K. – Abd. Kadir, Z. – Amer, N.H.
Intelligent Material Systems and Structures, vol. 27, no. 7, p. 
849-858, DOI:10.1177/1045389X15595291.
[8] Rahmat, M.S., Hudha, K., Kadir, Z.A., Amer, N.H., Nor, N.M., 
Choi, S.B. (2021). A hybrid skyhook active force control for 
impact mitigation using magneto-rheological elastomer 
isolator. Smart Materials and Structures, vol. 30, no. 2, art. ID 
25043, DOI:10.1088/1361-665X/abd911.
[9] Yang, J., Sun, S., Tian, T.F., Li, W., Du, H., Alici, G., Nakano, 
M. (2016). Development of a novel multi-layer MRE isolator 
for suppression of building vibrations under seismic events. 
Mechanical Systems and Signal Processing, vol. 70-71, p. 
811-820, DOI:10.1016/j.ymssp.2015.08.022.
[10] Norouzi, M., Sajjadi Alehashem, S.M., Vatandoost, H., Ni, Y.Q., 
Shahmardan, M.M. (2016). A new approach for modeling of 
magnetorheological elastomers. Journal of Intelligent Material 
Systems and Structures, vol. 27, no. 8, p. 1121-1135, 
DOI:10.1177/1045389X15615966.
[11] Yu, Y., Hoshyar, A.N., Li, H., Zhang, G., Wang, W. (2021). 
Nonlinear characterization of magnetorheological elastomer-
based smart device for structural seismic mitigation. 
International Journal of Smart and Nano Materials, vol. 12, 
no. 4, p. 390-428, DOI:10.1080/19475411.2021.1981477.
[12] Yu, Yang, Yousefi, A.M., Yi, K., Li, J., Wang, W., Zhou, X. (2021). 
A new hybrid model for MR elastomer device and parameter 
identification based on improved FOA. Smart Structures 
and Systems, vol. 28, no. 5, p. 617-629, DOI:10.12989/
sss.2021.28.5.617.
[13] Rossi, A., Orsini, F., Scorza, A., Botta, F., Belfiore, N.P., Sciuto, 
S.A. (2018). A review on parametric dynamic models of 
magnetorheological dampers and their characterization 
methods. Actuators, vol. 7, no. 2, art. ID 16, DOI:10.3390/
act7020016.
[14] Kwon, S.H., Lee, J.H., Choi, H.J. (2018). Magnetic particle filled 
elastomeric hybrid composites and their magnetorheological 
response. Materials, vol. 11, no. 6, art. ID 1040, DOI:10.3390/
ma11061040.
[15] Alias, N.F., Muthalif, A.G.A., Arpan, K.A.M., Nordin, N.H.D. 
(2018). Experimental investigation of static properties of 
magnetorheological elastomer. Iranian Journal of Science and 
Technology - Transactions of Mechanical Engineering, vol. 42, 
p. 185-197, DOI:10.1007/s40997-017-0081-5.
[16] Kang, S.S., Choi, K., Nam, J.-D., Choi, H.J. (2020). 
Magnetorheological elastomers: Fabrication, characteristics, 
and applications. Materials, vol. 13, no. 20, art. ID 4597, 
DOI:10.3390/ma13204597.
[17] Sobri, N.S., Lin, P.G., Haniffah, N.A., Kadir, Z.A., Hudha, K., 
Rahmat, M.S. (2021). Impact testing of magnetorheological 
elastomer based matrix material and magnetic particle. IEEE 
9th Conference on System, Process and Control, p. 70-74, 
DOI:10.1109/ICSPC53359.2021.9689137.
[18] Brancati, R., Di Massa, G., Di Vaio, M., Pagano, S., Santini, S. 
(2019). Experimental investigation on magneto-rheological 
elastomers. Mechanisms and Machine Science. Springer, 
Cham, vol. 68, p. 441-448, DOI:10.1007/978-3-030-03320-
0_48.
[19] Siddique, N., Adeli, H. (2016). Gravitational search algorithm 
and its variants. International Journal of Pattern Recognition 
and Artificial Intelligence, vol. 30, no. 8, art. ID 1639001, 
DOI:10.1142/S0218001416390018.
[20] Rashedi, E., Rashedi, E., Nezamabadi-pour, H. (2018). A 
comprehensive survey on gravitational search algorithm. 
Swarm and Evolutionary Computation, vol. 41, p. 141-158. 
DOI:10.1016/j.swevo.2018.02.018.
[21] Wang, Y., Gao, S., Zhou, M., Yu, Y. (2020). A multi-layered 
gravitational search algorithm for function optimization and 
real-world problems. IEEE/CAA Journal of Automatica Sinica, 
vol. 8, no. 1, p. 94-109, DOI:10.1109/JAS.2020.1003462.
[22] Subari, M.A., Hudha, K., Kadir, Z.A., DArdin, S.M.F.S.M., 
Amer, H.H. (2022). Path following control of tracked vehicle 
using modified sup controller optimized with particle swarm 
optimization (PSO). International Journal of Dynamics and 
Control, vol. 10, p. 1461-1470, DOI:10.1007/s40435-021-
00900-6.
[23] Mohd. Alawi, A.H., Hudha, K., Kadir, Z.A., Amer, N.H. (2023). 
Hysteresis behavior modeling of magnetorheological 
elastomers under impact loading using a multilayer 
exponential-based preisach model enhanced with particle 
swarm optimization. Polymers, vol. 15, no. 9, art. ID 2145, 
DOI:10.3390/polym15092145.