Dynamic Analysis of Sphere-Like Iron Particles Based Magnetorheological Damper for Waveform-Generating Test System

In this study, a new double pulse waveform-generating test system with an integrated magnetorheological (MR) damper is proposed. Since the total shear stress of MR fluid can be varied according to the shape of particles, sphere-like iron particles-based MR fluid is filled into the MR damper. The test system consists of a velocity generator, three masses (impact, test, and dummy), a spring, and an MR damper. To tune the double pulse waveform profile, a damping force model is constructed to determine the fundamental parameters of the simulator. Then, the first and second shock waveform profiles are analyzed to solve the governing equation of motions representing the damping force and velocity. The mathematical model of the MR damper is formulated and applied to a simulator with a graphical user interface programmed using MATLAB. The effectiveness of the proposed simulator-featuring controllable MR damper is demonstrated by comparing the simulation and experimental results.


Introduction
The first discovery of a magnetorheological (MR) fluid was made by Jacob Rainbow in the late 1940s [1,2]. Typically, MR fluids are composed of paramagnetic or soft ferromagnetic particles (~0.03-10 µm) dispersed in a carrier fluid. The MR particles are made of pure iron, carbonyl iron, or cobalt powder, and the carrier fluid is a non-magnetic, organic, or aqueous liquid, usually a silicone or mineral oil [3]. Without a magnetic field, the MR particles are randomly distributed in the carrier fluid. However, in magnetic field, they undergo rheological phase changes according to the magnitude of the field as shown in Figure 1a. The MR particles attain a dipole moment aligned with the external magnetic field and form chains. This chain formation causes an additional yield stress in the fluid, which can be rapidly and continuously controlled by magnetic input. As shown in Figure 1b, it is generally known that maximum yield stress of MR fluid is approximately 60-80 kPa. The power requirement for magnetic inputs is up to 50 Watts. In addition, it has been reported that the response time of commercial MR fluid is below 20 ms. As a result, MR fluid-based mechanical systems have Int. J. Mol. Sci. 2020, 21, 1149 2 of 16 salient advantages such as continuously adjustable dynamic range and fast response. However, the performance of MR fluid can be varied by the particle shape of MR fluid. As shown in Figure 1c, the shapes of particles are mostly sphere shape type and plate-lie shape type. Some of the research has analyzed the effect of particle shape of MR fluid [4][5][6][7]. From the previous work, it is revealed that sphere-like iron particles-based MR fluid is proper for MR damper, generating high damping force [8]. So, sphere-like iron particles-based MR fluid (RMS Corp., T5-500) is used in this work.
shocks. G. Wang et al. [21] investigated a vertical heavyweight shock test machine. Dynamic modeling and a mechanism analysis verified that the proposed test machine can produce shock acceleration wave profiles consistent with the MIL-S-901D criteria (via computer simulation). C. W. Shul et al. [22][23][24] carried out a dynamic design analysis for heavy dual shock generation systems. However, until now, a complete determination of the necessary parameter values of testing machines has not been possible because of the nonlinearity of large-scale MR dampers.
In this study, the results of a performance evaluation for a single unit large-scale MR damper were used to establish a mathematical model that can express the nonlinearity and to determine the set values of each element needed to realize the desired waveform. To obtain the desired shock wave magnitude and duration for each waveform in the double pulse waveform, the MR damper was employed to modify the shock Owing to the properties of MR fluids, MR dampers have inherent advantages such as fast response time and continuously controllable damping force. MR dampers are used in application devices, including semi-active controllable devices, such as shock absorbers for vehicle suspension systems [9][10][11][12][13][14], landing gears [15,16], and seismic dampers [17,18]. At present, shock wave profile generators based on MR dampers are actively researched. Many studies have been undertaken to evaluate the durability of on-board ship equipment using random shock waves generated from underwater explosions [19]. W. Zhaodong et al. [20] suggested a dual-wave shock-testing machine to simulate a non-contact underwater explosion and test the survivability of shipboard equipment to shocks. G. Wang et al. [21] investigated a vertical heavyweight shock test machine. Dynamic modeling and a mechanism analysis verified that the proposed test machine can produce shock acceleration wave profiles consistent with the MIL-S-901D criteria (via computer simulation). C. W. Shul et al. [22][23][24] carried out a dynamic design analysis for heavy dual shock generation systems. However, until now, a complete determination of the necessary parameter values of testing machines has not been possible because of the nonlinearity of large-scale MR dampers.
In this study, the results of a performance evaluation for a single unit large-scale MR damper were used to establish a mathematical model that can express the nonlinearity and to determine the set values of each element needed to realize the desired waveform. To obtain the desired shock wave magnitude and duration for each waveform in the double pulse waveform, the MR damper was employed to modify the shock magnitude and waveform. Then, MATLAB was used to program a graphical user interface (GUI) for the waveform prediction simulator, to verify the validity of the proposed system. Finally, the effectiveness of the proposed simulator-featuring controllable MR damper is demonstrated by comparing the simulated and experimental results.

Double Pulse Waveform Shock Test Setup
Since the full-scale submarine or naval vessel testing for heavy-weight components requires a lot of time and cost, BV-43/85 (developed in Germany) rules that a double half-sine acceleration profile, which can represent the shock loading caused by an underwater explosion, is applied to shipboard or submarine components [25,26]. As shown in Figure 2, the acceleration profile consists of three parts. In BV-43/85, after first half-sine acceleration waveform (Part I) strikes the test mass, second negative half-sine acceleration waveform (Part II) and very small residual vibration should be applied to test mass (Part III). In particular, the components are mounted on test mass. In order to dissipate the impact energy and realize the required acceleration waveform profile, MR damper with large power dissipation capacity is utilized. After striking, the impact mass is moved in the opposite direction to the test mass due to the elastic energy of programmer and air spring. Also, when impact mass is decoupled with test mass, the impact mass ideally moves at constant velocity. The positive half-wave (Part I) shows the effect of the magnitude of the impact, and the negative half-wave (Part II) and Part III are related to the energy dissipation system. In Figure 2, the magnitude of the acceleration in Part II is lower than that in Part I, and the magnitude of the residual vibration (Part III) is zero. Figure 3 shows a schematic diagram of the shock test experimental setup used in this study. The set up comprises a hydraulic velocity generator, impact mass, shock spring programmer, test mass, and MR damper. The specimen is attached to the test mass. In the impact-testing machine, the hydraulic velocity generator accelerates the impact mass, which then strikes the programmer shock spring and test mass with a uniform motion. The high-impact load and velocity are then applied to the MR damper, which is located between the test mass and a fixed wall. The MR damper is used to reduce the transmitted reaction force (Part II) and residual vibration (Part III). In other words, the generation of a shock-wave profile, particularly for Parts II and III, is the primary function of the MR damper.   The horizontal hydraulic double pulse shock tester is composed of the velocity generator, impact mass, shock programmer, test table, and MR damper. The shock test process was as follows. First, the velocity generator accelerates the impact mass at a specified velocity. The impact mass, accelerated to the maximum velocity, travels at constant velocity and impacts the shock programmer installed in front of the test table. After the impact, the impact mass bounces off and the impact acceleration of the test table is altered by the MR damper, located between the test table and the fixed wall. The impact acceleration is measured by an accelerometer and load cell, and the internal pressure of the MR damper is measured using a pressure sensor. The impact acceleration of the test table in this structure typically has the following two waveform types: the primary waveform, with a short period and positive acceleration value due to the impact between the accelerated impact mass and the test table, and the secondary impact acceleration, which is altered by the MR damper after the primary waveform. Figure 4 shows a photo of the MR damper and test equipment.

Double Pulse Shock Test Results and Analysis
To investigate the dynamic behavior and impact acceleration variation properties of the MR damper at high velocity and high mass, impact experiments were conducted using the horizontal hydraulic double pulse shock test equipment. Depending on the objective of the experiment, the damping force tuning properties of MR damper were evaluated by varying the current applied to the MR damper while fixing the impact mass, impact velocity, shock programmer stiffness, and test table mass. The damping force produced by the MR damper was plotted in the time domain and analyzed in the frequency domain, as shown in Figure 5. When high-mass structures impact at high velocity, various structural vibrations can be plotted simultaneously with major curves, as shown in Figures 6−9. The aim of this experiment was to obtain the force-velocity (F-V) diagram of the MR damper. However, when structural vibrations are mixed together in this process, the accurate modeling becomes difficult because a clear graph cannot be obtained. Thus, the structural vibration signal in the high frequency range, beyond that of the major waveform, needs to be removed. To accomplish this, a hamming window with a cut-off frequency of 500 Hz and a sampling frequency of 5120 Hz, and a finite impulse-response (FIR) low-pass filter were used for signal processing. The order of FIR filer is 100. The issue of time delay from the higher order was resolved and reorganized based on the impact velocity, and the final results are shown in Figures 6−9 and Table 1. The impact and test table masses were 140 and 280 kg, respectively, and the stiffness of the shock programmer was 4.6 MN/m. The results of the damping force of the MR damper in the time domain were measured using a load cell placed between the test table and the MR damper piston rod end. The results of the piston rod velocity in the time domain were measured using an accelerometer installed on the table mass. F-V curve results were obtained by plotting the measured damping force in the y-axis and the velocity in the x-axis.
The following observations were made from the shock experiment results. As the impact velocity increased, the damping force due to the viscous term of the fluid, or that of the MR damper at 0 A, increased proportionally, but the ratios of the magnetic field to the controllable damping force due to the MR yield stress and that due to the viscous term decreased. As the velocity increased, the effect of the damping force due to the magnetic field decreased, because the duration of the fluid flow through the orifice within the MR damper was shorter than the time required to observe the MR effect. Overall, the damping force increased but the controllability decreased.

Double Pulse Shock Test Results and Analysis
To investigate the dynamic behavior and impact acceleration variation properties of the MR damper at high velocity and high mass, impact experiments were conducted using the horizontal hydraulic double pulse shock test equipment. Depending on the objective of the experiment, the damping force tuning properties of MR damper were evaluated by varying the current applied to the MR damper while fixing the impact mass, impact velocity, shock programmer stiffness, and test table mass. The damping force produced by the MR damper was plotted in the time domain and analyzed in the frequency domain, as shown in Figure 5. When high-mass structures impact at high velocity, various structural vibrations can be plotted simultaneously with major curves, as shown in Figures 6-9. The aim of this experiment was to obtain the force-velocity (F-V) diagram of the MR damper. However, when structural vibrations are mixed together in this process, the accurate modeling becomes difficult because a clear graph cannot be obtained. Thus, the structural vibration signal in the high frequency range, beyond that of the major waveform, needs to be removed. To accomplish this, a hamming window with a cut-off frequency of 500 Hz and a sampling frequency of 5120 Hz, and a finite impulse-response (FIR) low-pass filter were used for signal processing. The order of FIR filer is 100. The issue of time delay from the higher order was resolved and reorganized based on the impact velocity, and the final results are shown in Figures 6-9 and Table 1. The impact and test table masses were 140 and 280 kg, respectively, and the stiffness of the shock programmer was 4.6 MN/m. The results of the damping force of the MR damper in the time domain were measured using a load cell placed between the test table and the MR damper piston rod end. The results of the piston rod velocity in the time domain were measured using an accelerometer installed on the table mass. F-V curve results were obtained by plotting the measured damping force in the y-axis and the velocity in the x-axis.
The following observations were made from the shock experiment results. As the impact velocity increased, the damping force due to the viscous term of the fluid, or that of the MR damper at 0 A, increased proportionally, but the ratios of the magnetic field to the controllable damping force due to the MR yield stress and that due to the viscous term decreased. As the velocity increased, the effect of the damping force due to the magnetic field decreased, because the duration of the fluid flow through the orifice within the MR damper was shorter than the time required to observe the MR effect. Overall, the damping force increased but the controllability decreased.

Damping F-V Curve Modeling
The magnetorheological fluid is the suspension of iron particles dispersed in a non-magnetic medium. The behavior of the MR fluid without magnetic field are similar with that of Newtonian fluid, but the unrestrained particles transform to a chain-like structuring under an external magnetic field. A relationship between magnetic input and shear stress of MR fluid has been researched and established as follows [27,28]: where τ and η are the total shear stress and viscosity constant of the MR fluid, respectively.
. γ is the relative shear rate. τ y (B) is the dynamic yield stress of the MR fluid, which is a function of the magnetic field intensity. To generate magnetic field, a DC current input is applied to coil. Accordingly, the total shear stress of MR fluid is affected by the current input and relative shear rate. Especially, some research reported that the viscosity can be tuned by the magnetic field intensity and magnetic field applying time for nanoparticle-added MR fluid [29,30]. Then, the damping force of MR damper is determined by the total shear stress of MR fluid and geometric configuration. However, when big shock is given to MR damper, the response of MR damper is different with that of MR damper with relative low piston velocity. Additionally, the shear rate of MR fluid is dependent on the piston velocity of MR damper. It is judged that the relationship between dynamic yield stress, shear rate and total shear stress is nonlinear. So, to predict the response of MR damper, an empirical formula is constructed in this work. The experimental results are divided into intervals with different velocity ranges, and a linear function modeling based on piston velocity and current input is conducted for each interval. The empirical formula can be utilized as the command law for generating varied double pulse waveform. Figure 10 shows the F-V curve obtained in Test 4, divided into intervals for modeling. It can be seen from Figure 10 that the F-V diagram obtained from the experiments are significantly nonlinear. Because it is difficult to model this with one equation, the method of modeling intervals was adopted. The velocity intervals were 0.0, 1.0, 1.8, 2.0, 2.5, 3.25, 3.5, and 4.5 m/s based on the extent of the F-V curve slope variation. Each interval was expressed with a linear function of the velocity. A simplified form of the function is: where, c(i) and d(i) are coefficients that change according to the input current, and F(v) is the damping force of the MR damper, which varies according to the velocity ν and the input current. Table 2 shows the c(i) and d(i) values obtained. Using the coefficients shown in Table 2, depending on whether the velocity is positive or negative, the damping force can be expressed as follows: Figure 11 shows the MR damper damping force curve modeling result obtained from the above equations. Figure 11a shows a comparison of the experiment results, divided into intervals, with the modeling results; the graph on the right shows the damping force results obtained using Equations (3) and (4) with input velocities between −4.5 m/s and 4.5 m/s. In Figure 11a, the solid black line is the curve obtained from the modeling, and the other color-lines are obtained from the experiments. As shown in Figure 11a, the simulation and experimental curves are very similar. Figure 11b shows a continuous curve after checking that the next modeling value and the experiment value matched. A linear function was used from −1 to 1 m/s because hysteresis characteristics in the low velocity range were not considered, as the goal was to express the high velocity range using the MR damper modeling. Moving on to the maximum velocity portion, it can be observed that, similar to the trend in the experimental result, the controllable range was approximately constant with increasing velocity, but the minimum damping force increased proportionally, indicating a lower control performance. and (4) with input velocities between −4.5 m/s and 4.5 m/s. In Figure 11 (a), the solid black line is the curve obtained from the modeling, and the other color-lines are obtained from the experiments. As shown in Figure 11 (a), the simulation and experimental curves are very similar. Figure 11 (b) shows a continuous curve after checking that the next modeling value and the experiment value matched. A linear function was used from -1 to 1 m/s because hysteresis characteristics in the low velocity range were not considered, as the goal was to express the high velocity range using the MR damper modeling. Moving on to the maximum velocity portion, it can be observed that, similar to the trend in the experimental result, the controllable range was approximately constant with increasing velocity, but the minimum damping force increased proportionally, indicating a lower control performance.

Waveform Prediction Simulator
A waveform prediction simulator is a program that calculates the input current considering the characteristics of the MR damper and the parameter values of the shock system when there is a desired waveform in a large-scale high-mass impact system. To construct such a program, the dynamics equations of the impact system have to be calculated and the MR damper model has to be applied. Finally, the system has to be programmed.

Dynamic Modeling of Dual Shock Waveform: Primary Shock Waveform
When the impact mass strikes the test table, a primary shock and vibration are generated. With air as the spring, corresponding to the dummy mass (m 3 ), the governing equation of the primary shock waveform can be derived as follows: Because x 1 , x 2 , and x 3 are the sum of sinusoidal waves, Equation (5) can be transformed as follows: From the determinant of Equation (6), the natural frequencies are obtained as follows: Then, the general solution of Equation (5) can be expressed as follows: 1 cos(ω 1 t + ϕ 1 ) + X 1 cos(ω 2 t + ϕ 2 ) + X By using the initial conditions, Equation (5) can be re-written as follows: From Equation (12), the following angle value is obtained: Finally, the particular solution of Equation (5) is obtained from Equations (10) and (13): x 2 (t) = −r 12_1 X (1) ..
x 2 (t) = r 12_1 ω 1 2 X (1) Because the impact mass is decoupled after striking the test mass and ω 2 is very low, the frequency of the primary shock wave form is determined approximately by ω 3 . Accordingly, the shock duration time and maximum acceleration of the primary shock waveform can be expressed as follows: Max. Acc. = ..

Dynamic Modeling of Dual Shock Waveform: Secondary Shock Waveform
A similar approach is applied to derive the dynamic model of the secondary shock waveform. The difference between the primary and secondary waveforms is the initial conditions. From Equation (16) with T = π/ω 3 , the initial conditions of the secondary waveform can be expressed as follows: where D 2init. , V 2init. , D 3init. , and V 3init. are the initial displacement and velocity of the second waveform for test mass (m 2 ) and dummy mass (m 3 ), respectively. The governing equation of the secondary shock waveform can be derived as follows: Equation (21) can be transformed as follows: Then, the general solution of Equation (21) can be expressed as follows: From the initial conditions, Equations (19) and (20), the following equations can be: From Equation (24), the expression for the angle is obtained: Finally, the particular solution of Equation (21) is obtained from Equations (23) and (25): .. ..
Because the first natural frequency of the secondary waveform is zero, the frequency of the secondary shock wave form is ω 2 . Accordingly, the shock duration time and maximum acceleration of the secondary shock waveform can be expressed as follows: Max.Acc. = ..

Waveform Prediction Simulator Program Structure
The MATLAB GUI program was used to produce the waveform prediction simulator. In the structure, the equations for the primary and secondary waveforms analyzed in Sections 3.1 and 3.2 were all input and solved to obtain the required values. The procedure is described in the following paragraph. The impact and test table masses were fixed and the material properties of the variables that could be obtained analytically were determined, including the impact velocity, dummy mass, shock programmer stiffness, and air spring stiffness. To integrate these into the MR damper model, the values were changed based on the obtained dummy mass. The waveform properties were determined by solving the response of the entire impact system, including the MR damper model, using the Runge Kutta method, based on the following equations: Multiple cases that yielded the desired waveform were identified and that with the most reasonable variables was selected. Figure 12 shows the proposed GUI program and flow chart. First, the waveform properties of the impact duration and maximum impact acceleration for the primary and secondary waveforms are input, followed by the setting of the impact, test, and dummy mass ranges for the system. In <1>, the variables that can realize first shock waveform are determined, except for the MR damper. Based on this, the numerical analysis, including the MR damper, is carried out; the various cases that approach the waveform properties are listed in process <2>. Among the listed cases, the proper one is manually selected and the predicted waveform is calculated. The secondary waveform time and acceleration, where the residual vibrations are attenuated, are expressed according to the damping force model of the MR damper (Equations (2)-(4)), in process <3>.   Figures 13 and 14 show the experimental and simulation results obtained with the waveform properties period and maximum acceleration set to 0.01 s and 590 m/s 2 (60G), respectively, for the primary waveform and 0.06 s for the period of secondary waveform. The parameter values obtained through the waveform prediction program were used to carry out the simulation. From the simulator, the masses of impact, test table, and dummy were determined to be 140, 200, and 100 kg, respectively. The stiffness of programmer shock spring and air spring were 4.6 MN/m and 0.39 MN/m, respectively. The magnitude of the impact velocity was 6.1 m/s. The simulation results were found to be in relatively good agreement with the experimental results. From the simulation experimental results, it can be inferred that the proposed dynamic modeling of the test system and experimental damping force model of MR damper are effective for predicting shockwave form. From Figures 13 and 14, the magnitude and duration time of the second shockwave are increased according to current input. This is why the reaction force of the MR damper is proportional to the damping force of the MR damper. Also, the dynamic energy should be dissipated during part III, owing to the MR damper. In particular, the dynamic energy of the testing system is proportional to the acceleration of test mass. From Figure 14, it can be known that the Part III acceleration results for 0.5 A is reduced from 84.15 m/s 2 to 43.87 m/s 2 in 0.149 s.
However, the measured magnitude of the first shockwave is different from that of the predicted shockwave. During the impact between impact and test mass, the high-impact velocity is unequally applied to test mass. This is the reason for the difference at the first shockwave. The maximum error value of the first shockwave is 173.14 m/s 2 at 0 A. The other tangible difference is the duration time of the second shockwave (t 2 ). The maximum difference value of the second shockwave is 0.01 s at 0 A. This is related with the nonlinearity of the elastic and damping force for sudden impact and slow response time of the MR damper. In addition, the structure of the proposed program involves determining the value by varying the dummy mass, which sometimes becomes too large, resulting in a lower MR damper efficiency. Thus, in future studies, the waveform should be found by varying other material properties while leaving the maximum mass unchanged.

Conclusions
In this study, we proposed an MR damper modeling and a waveform prediction program to simulate the double pulse from an underwater explosion. Impact experiments were carried out to obtain the damping F-V curve, and that with the highest velocity was modeled mathematically. For this, the curve with high nonlinearity was divided into intervals with different velocity ranges, with rapid variations in the damping force. A linear function modeling was conducted for each interval. Then, the simulation results were compared with the actual impact experiment results. From the comparison results between simulation and experimental results, the proposed dynamic modeling of the system and damping are quite self-explanatory, showing good agreement with experimental results. In addition, it has been confirmed that MR damper featuring sphere-like iron particle shows a high damping force up to 200 kN. Since the performance on damping force, durability, sedimentation, and dispersible stability of MR fluid can be varied according to the shape of particles, the performance of MR dampers with different particle shapes will be evaluated as a second phase of this work in the future.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
In this study, we proposed an MR damper modeling and a waveform prediction program to simulate the double pulse from an underwater explosion. Impact experiments were carried out to obtain the damping F-V curve, and that with the highest velocity was modeled mathematically. For this, the curve with high nonlinearity was divided into intervals with different velocity ranges, with rapid variations in the damping force. A linear function modeling was conducted for each interval. Then, the simulation results were compared with the actual impact experiment results. From the comparison results between simulation and experimental results, the proposed dynamic modeling of the system and damping are quite self-explanatory, showing good agreement with experimental results. In addition, it has been confirmed that MR damper featuring sphere-like iron particle shows a high damping force up to 200 kN. Since the performance on damping force, durability, sedimentation, and dispersible stability of MR fluid can be varied according to the shape of particles, the performance of MR dampers with different particle shapes will be evaluated as a second phase of this work in the future.