Impact Behaviour Modelling of Magnetorheological Elastomer Using a Non-Parametric Polynomial Model Optimized with Gravitational Search Algorithm

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.


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 MREbased 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 neurofuzzy 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 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 forcedisplacement 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.

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].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].

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-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.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  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.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.

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.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): 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 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].

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.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).
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.
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.
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).
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).
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, 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.

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. -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].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.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.

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:

X
The interpolation model for MRE was executed in the Matlab-Simulink environment using input currents 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].

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.

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 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.

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.

Fig. 9
Fig. 9 shows the superimposed forcedisplacement 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

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 Aof 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 .
Fig. 10.The structure of the interpolation current prediction model 4.1 Validation of the Dual-Acting MRE Behaviours for The Middle Current

Fig. 9 .
Fig. 9. Comparison of force-displacement curve for simulation at varying current 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:

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

Table 1 .
Share of MRE samples

Table 2 .
Experimental parameters set for impact test

Table 3 .
Initial simulation parameters for GSA

Table 4 .
Optimized model parameter for each input current

Table 5 .
Force prediction error for input current 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

Table 6 .
Force prediction error for the interpolation current