Open AccessArticle

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

Jong-Seok Oh

^1,*

Chang Won Shul

²,

Tae Hyeong Kim

²,

Tae-Hoon Lee

³,

Sung-Wan Son

⁴ and

Seung-Bok Choi

^3,*

Division of Mechanical and Automotive Engineering, Kongju National University, Chungnam 31080, Korea

Agency for Defense Development, Daejeon 305-600, Korea

Department of Mechanical Engineering, Smart Structures and Systems Laboratory, Inha University, Incheon 22212, Korea

⁴

RMS Technology, Pungse Industrial Complex A17-6, Boseong-ri, Chungcheongnam-do 31217, Korea

Authors to whom correspondence should be addressed.

Int. J. Mol. Sci. 2020, 21(3), 1149; https://doi.org/10.3390/ijms21031149

Submission received: 30 November 2019 / Revised: 25 January 2020 / Accepted: 8 February 2020 / Published: 9 February 2020

(This article belongs to the Special Issue Magnetic-Responsive Molecular Particles Based Smart Materials: Model, Characterization and Applications)

Download

Browse Figures

Figure 1
MR fluid: (a) Phenomenological characteristics; (b) Yield stress results; (c) Shape of particle. "> Figure 2
Shape of the double pulse waveform profile. "> Figure 3
Schematic diagram of the horizontal hydraulic double pulse shock test setup. "> Figure 4
Photograph of the double pulse shock test setup: (a) Manufacture MR damper; (b) Experimental apparatus. "> Figure 4 Cont.
Photograph of the double pulse shock test setup: (a) Manufacture MR damper; (b) Experimental apparatus. "> Figure 5
Measured damping force in the time and frequency domains: (a) Measured damping force in the time domain; (b) Measured damping force in the frequency domain. "> Figure 6
MR damper experimental results at an impact velocity of 2.5 m/s (Test 1): (a) Damping force Piston velocity; (b) Piston Velocity; (c) F-V diagram of MR damper. "> Figure 7
MR damper experimental results at an impact velocity of 5.37 m/s (Test 2): (a) Damping force Piston velocity; (b) Piston Velocity; (c) F-V diagram of MR damper. "> Figure 8
MR damper experimental results at an impact velocity of 7.75 m/s (Test 3): (a) Damping force Piston velocity; (b) Piston Velocity; (c) F-V diagram of MR damper. "> Figure 9
MR damper experimental results at an impact velocity of 9.12 m/s (Test 4): (a) Damping force; (b) Piston velocity; (c) force-velocity (F-V) diagram of MR damper. "> Figure 10
F-V curve (Test. 4). "> Figure 11
Estimation results of F-V model for MR damper: (a) Comparison results between experimental and simulator results; (b) Simulation F-V curve according to continuous current inputs. "> Figure 12
GUI of the waveform prediction simulator program: (a)Flow chart; (b) GUI program. "> Figure 12 Cont.
GUI of the waveform prediction simulator program: (a)Flow chart; (b) GUI program. "> Figure 13
Simulation result of the waveform prediction program. "> Figure 14
Experimental result of the waveform prediction program. ">

Review Reports Versions Notes

Abstract

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.

Keywords:

Magnetorheological (MR) Fluid; MR Damper; Double Waveform-generating test System; Shock-Wave Profile

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

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.

2. Force-Velocity Diagram Modeling of the High Capacity Shock MR Damper

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

2.2. 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 Figure 6, Figure 7, Figure 8 and Figure 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 Figure 6, Figure 7, Figure 8 and Figure 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.

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

τ = η \dot{γ} + τ_{y} (B)

(1)

where

τ

and

η

are the total shear stress and viscosity constant of the MR fluid, respectively.

\dot{γ}

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:

F (v) = c (i) v + d (i)

(2)

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:

i f v e l \geq 0, F (v) = c (i) (v - v e l .) + d (i)

(3)

i f v e l \leq 0, F (v) = - c (i) (- v - v e l .) - d (i)

(4)

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.

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

3.1. 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₃), the governing equation of the primary shock waveform can be derived as follows:

\begin{array}{l} m_{1} {\ddot{x}}_{1} + k_{1} x_{1} - k_{1} x_{2} = 0 \\ m_{2} {\ddot{x}}_{2} - k_{1} x_{1} + (k_{1} + k_{2}) x_{2} - k_{2} x_{3} = 0 \\ m_{3} {\ddot{x}}_{3} - k_{2} x_{2} + k_{2} x_{3} = 0 \end{array}

(5)

Because x₁, x₂, and x₃ are the sum of sinusoidal waves, Equation (5) can be transformed as follows:

[\begin{matrix} (k_{1} - m_{1} ω^{2}) & - k_{1} & 0 \\ - k_{1} & (k_{1} + k_{2} - m_{2} ω^{2}) & - k_{2} \\ 0 & - k_{2} & (k_{2} - m_{3} ω^{2}) \end{matrix}] {\begin{matrix} X_{1} \\ X_{2} \\ X_{3} \end{matrix}} = {\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}}

(6)

From the determinant of Equation (6), the natural frequencies are obtained as follows:

ω_{1} = 0

(7)

ω_{2} = \frac{k_{2} m_{1} m_{2} + k_{1} m_{1} m_{3} + k_{2} m_{1} m_{3} + k_{1} m_{2} m_{3} + \sqrt{4 m_{1} m_{2} m_{3} (- k_{1} k_{2} m_{1} - k_{1} k_{2} m_{2} - k_{1} k_{2} m_{3}) + {(k_{2} m_{1} m_{2} + k_{1} m_{1} m_{3} + k_{2} m_{1} m_{3} + k_{1} m_{2} m_{3})}^{2}}}{2 m_{1} m_{2} m_{3}}

(8)

ω_{3} = \frac{k_{2} m_{1} m_{2} + k_{1} m_{1} m_{3} + k_{2} m_{1} m_{3} + k_{1} m_{2} m_{3} - \sqrt{4 m_{1} m_{2} m_{3} (- k_{1} k_{2} m_{1} - k_{1} k_{2} m_{2} - k_{1} k_{2} m_{3}) + {(k_{2} m_{1} m_{2} + k_{1} m_{1} m_{3} + k_{2} m_{1} m_{3} + k_{1} m_{2} m_{3})}^{2}}}{2 m_{1} m_{2} m_{3}}

(9)

Then, the general solution of Equation (5) can be expressed as follows:

\begin{array}{l} x_{1} (t) = X_{1}^{(1)} \cos (ω_{1} t + φ_{1}) + X_{1}^{(2)} \cos (ω_{2} t + φ_{2}) + X_{1}^{(3)} \cos (ω_{3} t + φ_{3}) \\ x_{2} (t) = r_{12_1} X_{1}^{(1)} \cos (ω_{1} t + φ_{1}) + r_{12_2} X_{1}^{(2)} \cos (ω_{2} t + φ_{2}) + r_{12_3} X_{1}^{(3)} \cos (ω_{3} t + φ_{3}) \\ x_{3} (t) = r_{13_1} X_{1}^{(1)} \cos (ω_{1} t + φ_{1}) + r_{13_2} X_{1}^{(2)} \cos (ω_{2} t + φ_{2}) + r_{13_3} X_{1}^{(3)} \cos (ω_{3} t + φ_{3}) \end{array}

(10)

\begin{array}{l} w h e r e \\ r_{12_1} = \frac{(k_{1} - m_{1} ω_{1}^{2})}{k_{1}}, r_{23_1} = \frac{k_{2}}{(k_{2} - m_{3} ω_{1}^{2})}, r_{13_1} = r_{12_1} \times r_{23_1} = \frac{k_{2} (k_{1} - m_{1} ω_{1}^{2})}{k_{1} (k_{2} - m_{3} ω_{1}^{2})}, \\ r_{12_2} = \frac{(k_{1} - m_{1} ω_{2}^{2})}{k_{1}}, r_{23_2} = \frac{k_{2}}{(k_{2} - m_{3} ω_{2}^{2})}, r_{13_2} = r_{12_2} \times r_{23_2} = \frac{k_{2} (k_{1} - m_{1} ω_{2}^{2})}{k_{1} (k_{2} - m_{3} ω_{2}^{2})}, \\ r_{12_3} = \frac{(k_{1} - m_{1} ω_{3}^{2})}{k_{1}}, r_{23_3} = \frac{k_{2}}{(k_{2} - m_{3} ω_{3}^{2})}, r_{13_3} = r_{12_3} \times r_{23_3} = \frac{k_{2} (k_{1} - m_{1} ω_{3}^{2})}{k_{1} (k_{2} - m_{3} ω_{3}^{2})} \end{array}

(11)

By using the initial conditions, Equation (5) can be re-written as follows:

\begin{array}{l} x_{1} (0) = 0 = X_{1}^{(1)} \cos φ_{1} + X_{1}^{(2)} \cos φ_{2} + X_{1}^{(3)} \cos φ_{3} \\ {\dot{x}}_{1} (0) = v_{i m p a c t} = - ω_{1} X_{1}^{(1)} \sin φ_{1} - ω_{2} X_{1}^{(2)} \sin φ_{2} - ω_{3} X_{1}^{(3)} \sin φ_{3} \\ x_{2} (0) = 0 = r_{12_1} X_{1}^{(1)} \cos φ_{1} + r_{12_2} X_{1}^{(2)} \cos φ_{2} + r_{12_3} X_{1}^{(3)} \cos φ_{3} \\ {\dot{x}}_{2} (0) = 0 = - ω_{1} r_{12_1} X_{1}^{(1)} \sin φ_{1} - ω_{2} r_{12_2} X_{1}^{(2)} \sin φ_{2} - ω_{3} r_{12_3} X_{1}^{(3)} \sin φ_{3} \\ x_{3} (0) = 0 = r_{13_1} X_{1}^{(1)} \cos φ_{1} + r_{13_2} X_{1}^{(2)} \cos φ_{2} + r_{13_3} X_{1}^{(3)} \cos φ_{3} \\ {\dot{x}}_{3} (0) = 0 = - ω_{1} r_{13_1} X_{1}^{(1)} \sin φ_{1} - ω_{2} r_{13_2} X_{1}^{(2)} \sin φ_{2} - ω_{3} r_{13_3} X_{1}^{(3)} \sin φ_{3} \end{array}

(12)

From Equation (12), the following angle value is obtained:

{\begin{matrix} φ_{1} \\ φ_{2} \\ φ_{3} \end{matrix}} = {\begin{matrix} \frac{π}{2} \\ \frac{π}{2} \\ \frac{π}{2} \end{matrix}}

(13)

Finally, the particular solution of Equation (5) is obtained from Equations (10) and (13):

\begin{array}{l} x_{1} (t) = - X_{1}^{(1)} \sin ω_{1} t - X_{1}^{(2)} \sin ω_{2} t - X_{1}^{(3)} \sin ω_{3} t \\ {\dot{x}}_{1} (t) = - ω_{1} X_{1}^{(1)} \cos ω_{1} t - ω_{2} X_{1}^{(2)} \cos ω_{2} t - ω_{3} X_{1}^{(3)} \cos ω_{3} t \\ {\ddot{x}}_{1} (t) = ω_{1}^{2} X_{1}^{(1)} \sin ω_{1} t + ω_{2}^{2} X_{1}^{(2)} \sin ω_{2} t + ω_{3}^{2} X_{1}^{(3)} \sin ω_{3} t \end{array}

(14)

\begin{array}{l} x_{2} (t) = - r_{12_1} X_{1}^{(1)} \sin ω_{1} t - r_{12_2} X_{1}^{(2)} \sin ω_{2} t - r_{12_3} X_{1}^{(3)} \sin ω_{3} t \\ {\dot{x}}_{2} (t) = - r_{12_1} ω_{1} X_{1}^{(1)} \cos ω_{1} t - r_{12_2} ω_{2} X_{1}^{(2)} \cos ω_{2} t - r_{12_3} ω_{3} X_{1}^{(3)} \cos ω_{3} t \\ {\ddot{x}}_{2} (t) = r_{12_1} ω_{1}^{2} X_{1}^{(1)} \sin ω_{1} t + r_{12_2} ω_{2}^{2} X_{1}^{(2)} \sin ω_{2} t + ω_{3}^{2} r_{12_3} X_{1}^{(3)} \sin ω_{3} t \end{array}

(15)

\begin{array}{l} x_{3} (t) = - r_{13_1} X_{1}^{(1)} \sin ω_{1} t - r_{13_2} X_{1}^{(2)} \sin ω_{2} t - r_{13_3} X_{1}^{(3)} \sin ω_{3} t \\ {\dot{x}}_{3} (t) = - r_{13_1} ω_{1} X_{1}^{(1)} \cos ω_{1} t - r_{13_2} ω_{2} X_{1}^{(2)} \cos ω_{2} t - r_{13_3} ω_{3} X_{1}^{(3)} \cos ω_{3} t \\ {\ddot{x}}_{3} (t) = r_{13_1} ω_{1}^{2} X_{1}^{(1)} \sin ω_{1} t + r_{13_2} ω_{2}^{2} X_{1}^{(2)} \sin ω_{2} t + ω_{3}^{2} r_{13_3} X_{1}^{(3)} \sin ω_{3} t \end{array}

(16)

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:

S h o c k d u r a t i o n t i m e = \frac{π}{ω_{3}}

(17)

M a x . A c c . = {\ddot{x}}_{2} (T) = r_{12_2} ω_{2}^{2} X_{1}^{(2)} \sin ω_{2} T + ω_{3}^{2} r_{12_3} X_{1}^{(3)} \sin ω_{3} T (T = \frac{π}{2 ω_{3}})

(18)

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

\begin{array}{l} D_{2 i n i t .} = - r_{12_1} X_{1}^{(1)} \sin ω_{1} T - r_{12_2} X_{1}^{(2)} \sin ω_{2} T - r_{12_3} X_{1}^{(3)} \sin ω_{3} T \\ V_{2 i n i t .} = - r_{12_1} ω_{1} X_{1}^{(1)} \cos ω_{1} T - r_{12_2} ω_{2} X_{1}^{(2)} \cos ω_{2} T - r_{12_3} ω_{3} X_{1}^{(3)} \cos ω_{3} T \end{array}

(19)

\begin{array}{l} D_{3 i n i t .} = - r_{13_1} X_{1}^{(1)} \sin ω_{1} T - r_{13_2} X_{1}^{(2)} \sin ω_{2} T - r_{13_3} X_{1}^{(3)} \sin ω_{3} T \\ V_{3 i n i t .} = - r_{13_1} ω_{1} X_{1}^{(1)} \cos ω_{1} T - r_{13_2} ω_{2} X_{1}^{(2)} \cos ω_{2} T - r_{13_3} ω_{3} X_{1}^{(3)} \cos ω_{3} T \end{array}

(20)

where

D_{2 i n i t .}, V_{2 i n i t .}, D_{3 i n i t .},

and

V_{3 i n i t .}

are the initial displacement and velocity of the second waveform for test mass (m₂) and dummy mass (m₃), respectively. The governing equation of the secondary shock waveform can be derived as follows:

\begin{array}{l} m_{2} {\ddot{x}}_{2} + k_{2} x_{2} - k_{2} x_{3} = 0 \\ m_{3} {\ddot{x}}_{3} - k_{2} x_{2} + k_{2} x_{3} = 0 \end{array}

(21)

Equation (21) can be transformed as follows:

[\begin{matrix} (k_{2} - m_{2} ω^{2}) & - k_{2} \\ - k_{2} & (k_{2} - m_{3} ω^{2}) \end{matrix}] {\begin{matrix} X_{2} \\ X_{3} \end{matrix}} = {\begin{matrix} 0 \\ 0 \end{matrix}}; ω_{1} = 0, ω_{2} = \sqrt{\frac{k_{2} (m_{2} + m_{3})}{m_{2} m_{3}}}

(22)

Then, the general solution of Equation (21) can be expressed as follows:

\begin{array}{l} x_{2} (t) = X_{2}^{(1)} \cos (ω_{1} t + φ_{1}) + X_{2}^{(2)} \cos (ω_{2} t + φ_{2}) \\ x_{3} (t) = r_{1} X_{2}^{(1)} \cos (ω_{1} t + φ_{1}) + r_{2} X_{2}^{(2)} \cos (ω_{2} t + φ_{2}) \\ w h e r e \\ r_{1} = \frac{X_{3}^{(1)}}{X_{2}^{(1)}} = \frac{- m_{2} ω_{1}^{2} + k_{2}}{k_{2}} = 1, r_{2} = \frac{X_{3}^{(2)}}{X_{2}^{(2)}} = \frac{- m_{2} ω_{2}^{2} + k_{2}}{k_{2}} = - \frac{m_{2}}{m_{3}} \end{array}

(23)

From the initial conditions, Equations (19) and (20), the following equations can be:

\begin{array}{l} x_{2} (0) = X_{2}^{(1)} \cos ϕ_{1} + X_{2}^{(2)} \cos ϕ_{2} = D_{2 i n i t .} \\ {\dot{x}}_{2} (0) = - ω_{1} X_{2}^{(1)} \sin ϕ_{1} - ω_{2} X_{2}^{(2)} \sin ϕ_{2} = V_{2 i n i t .} \\ x_{3} (0) = r_{1} X_{2}^{(1)} \cos ϕ_{1} + r_{2} X_{2}^{(2)} \cos ϕ_{2} = D_{3 i n i t .} \\ {\dot{x}}_{3} (0) = - ω_{1} r_{1} X_{2}^{(1)} \sin ϕ_{1} - ω_{2} r_{2} X_{2}^{(2)} \sin ϕ_{2} = V_{3 i n i t .} \end{array}

(24)

From Equation (24), the expression for the angle is obtained:

ϕ_{1} = \tan^{- 1} (\frac{- r_{2} V_{2 i n i t .} + V_{3 i n i t .}}{ω_{1} (r_{2} D_{2 i n i t .} - D_{3 i n i t .})}) = \frac{π}{2}, ϕ_{2} = \tan^{- 1} (\frac{r_{1} V_{2 i n i t .} - V_{3 i n i t .}}{ω_{2} (- r_{1} D_{2 i n i t .} + D_{3 i n i t .})})

(25)

Finally, the particular solution of Equation (21) is obtained from Equations (23) and (25):

\begin{array}{l} x_{2} (t) = X_{2}^{(2)} \cos (ω_{2} t + ϕ_{2}) \\ {\dot{x}}_{2} (t) = - ω_{2} X_{2}^{(2)} \sin (ω_{2} t + ϕ_{2}) \\ {\ddot{x}}_{2} (t) = - ω_{2}^{2} X_{2}^{(2)} \cos (ω_{2} t + ϕ_{2}) \end{array}

(26)

\begin{array}{l} x_{3} (t) = r_{2} X_{2}^{(2)} \cos (ω_{2} t + ϕ_{2}) \\ {\dot{x}}_{3} (t) = - ω_{2} X_{2}^{(2)} \sin (ω_{2} t + ϕ_{2}) \\ {\ddot{x}}_{3} (t) = - ω_{2}^{2} X_{2}^{(2)} \cos (ω_{2} t + ϕ_{2}) \end{array}

(27)

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:

S h o c k d u r a t i o n t i m e = \frac{π}{ω_{2}}

(28)

M a x . A c c . = {\ddot{x}}_{2} (T) = - ω_{2}^{2} X_{2}^{(2)} \cos (ω_{2} T + ϕ_{2}) (T = \frac{π}{2 ω_{2}})

(29)

3.3. 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 Section 3.1 and Section 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:

\begin{array}{l} m_{2} {\ddot{x}}_{2} + k_{2} x_{2} - k_{2} x_{3} = 0 \\ m_{3} {\ddot{x}}_{3} - k_{2} x_{2} + k_{2} x_{3} - {(F (v))}_{MR} = 0 \end{array}

(30)

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

Figure 13 and Figure 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² (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 Figure 13 and Figure 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² to 43.87 m/s² 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² at 0 A. The other tangible difference is the duration time of the second shockwave (t₂). 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.

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

Author Contributions

Conceptualization, C.W. Shul; methodology, T.H.K.; data curation, T.H.K.; validation, S.-W.S.; investigation, T.-H.L.; writing—original draft preparation, J.-S.O.; writing—review and editing, J.-S.O. and S.-B.C.; supervision, S.-B.C.; project administration, C.W.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

This work was supported by the Agency for Defense Development (UG130119JD).

Conflicts of Interest

The authors declare no conflict of interest.

References

Rabinow, J. The magnetic fluid clutch. Electr. Eng. 1948, 67, 1167. [Google Scholar] [CrossRef]
Rabinow, J. Magnetic Fluid Torque and Force Transmitting Device. U.S. Patent 2,575,360, 20 November 1951. [Google Scholar]
Choi, S.B.; Han, Y.M. Magnetorheological Fluid Technology: Applications in Vehicle Systems; CRC Press: Boca Raton, FL, USA, 2012. [Google Scholar]
De Vicente, J.; Klingenberg, D.J.; Hidalgo-Alvarez, R. Magnetorheological fluids: A review. Soft Matter 2011, 7, 3701–3710. [Google Scholar] [CrossRef]
Kor, Y.K.; See, H. The electrorheological response of elongated particles. Rheol. Acta 2010, 49, 741–756. [Google Scholar] [CrossRef]
De Vicente, J.; Vereda, F.; Segovia-Gutiérrez, J.P.; del Puerto Morales, M.; Hidalgo-Álvarez, R. Effect of particle shape in magnetorheology. J. Rheol. 2010, 54, 1337–1362. [Google Scholar] [CrossRef]
Shah, K.; Oh, J.S.; Choi, S.B.; Upadhyay, R.V. Plate-like iron particles based bidisperse magnetorheological fluid. J. Appl. Phys. 2013, 114, 213904. [Google Scholar] [CrossRef]
Jiang, Z.; Christenson, R. A comparison of 200 kN magneto-rheological damper models for use in real-time hybrid simulation pretesting. Smart Mater. Struct. 2011, 20, 065011. [Google Scholar] [CrossRef]
Choi, S.B.; Sung, K.G. Vibration control of magnetorheological damper system subjected to parameter variations. Int. J. Veh. Des. 2008, 46, 94–110. [Google Scholar] [CrossRef]
Lu, S.B.; Li, Y.N.; Choi, S.B.; Zheng, L.; Seong, M.S. Integrated control on MR vehicle suspension system associated with braking and steering control. Veh. Syst. Dyn. 2011, 49, 361–380. [Google Scholar] [CrossRef]
Choi, S.-B.; Han, Y.-M. MR seat suspension for vibration control of a commercial vehicle. Int. J. Veh. Des. 2003, 31, 202–215. [Google Scholar] [CrossRef]
Nguyen, Q.H.; Choi, S.B. Optimal design of MR shock absorber and application to vehicle suspension. Smart Mater. Struct. 2009, 18, 035012. [Google Scholar] [CrossRef]
Nguyen, Q.H.; Choi, S.B. Optimal design of a vehicle magnetorheological damper considering the damping force and dynamic range. Smart Mater. Struct. 2008, 18, 015013. [Google Scholar] [CrossRef]
Choi, S.B.; Li, W.; Yu, M.; Du, H.; Fu, J.; Do, P.X. State of the art of control schemes for smart systems featuring magneto-rheological materials. Smart Mater. Struct. 2016, 25, 043001. [Google Scholar] [CrossRef] [Green Version]
Choi, Y.T.; Wereley, N.M. Vibration control of a landing gear system featuring electrorheological/magnetorheological fluids. J. Aircr. 2003, 40, 432–439. [Google Scholar] [CrossRef]
Han, C.; Kim, B.G.; Choi, S.B. Design of a New Magnetorheological Damper Based on Passive Oleo-Pneumatic Landing Gear. J. Aircr. 2018, 55, 2510–2520. [Google Scholar] [CrossRef]
Seong, M.S.; Choi, S.B.; Kim, C.H. Design and performance evaluation of MR damper for integrated isolation mount. J. Intell. Mater. Syst. Struct. 2011, 22, 1729–1738. [Google Scholar] [CrossRef]
Han, C.; Choi, S.B.; Lee, Y.S.; Kim, H.T.; Kim, C.H. A new hybrid mount actuator consisting of air spring and magneto-rheological damper for vibration control of a heavy precision stage. Sens. Actuators A Phys. 2018, 284, 42–51. [Google Scholar] [CrossRef]
Kim, H.C.; Oh, J.S.; Choi, S.B. The field-dependent shock profiles of a magnetorhelogical damper due to high impact: An experimental investigation. Smart Mater. Struct. 2014, 24, 025008. [Google Scholar] [CrossRef]
Zhaodong, W.; Yu, W.; Lei, Z.; Jianye, D. Modeling and dynamic simulation of novel dual-wave shock test machine. In Proceedings of the 2011 International Conference on Fluid Power and Mechatronics, Beijing, China, 17–20 August 2011; IEEE: Piscataway, NJ, USA, 2011; pp. 898–901. [Google Scholar]
Wang, G.; Xiong, Y.; Tang, W. A novel heavy-weight shock test machine for simulating underwater explosive shock environment: Mathematical modeling and mechanism analysis. Int. J. Mech. Sci. 2013, 77, 239–248. [Google Scholar] [CrossRef]
Kim, T.H.; Shul, C.W.; Yang, M.S.; Lee, G.S. An analytic investigation on the implementation method for required signal of heavy weight shock test machine. In Proceedings of the Korean Society of Mechanical enginees, Annual Spring/Fall Conference, Gwangju, Korea, 11–14 November 2014; KSME: Seoul, Korea, 2014; pp. 400–405. [Google Scholar]
Shul, C.W.; Kim, T.H.; Yang, M.S.; Lee, G.S. Development of large-scale heavy weight shock testing system. In Proceedings of the Korean Society of Mechanical enginees, Annual Spring/Fall Conference, Gwangju, Korea, 11–14 November 2014; KSME: Seoul, Korea, 2014; pp. 377–382. [Google Scholar]
Kang, M.S.; Shul, C.W.; Kim, T.H.; Yang, M.S.; Song, W.K.; Lee, G.S. A design of large and heavy dual shock generation system. In Proceedings of the Korean Society of Mechanical enginees, Annual Spring/Fall Conference, Gwangju, Korea, 11–14 November 2014; KSME: Seoul, Korea, 2014; pp. 383–388. [Google Scholar]
Chnstiansen, M. Under Water Explosive Shock testing (UNDEX) of a subsea mateable electrical connector, the CM2000. In MTS/IEEE Oceans 2001. An Ocean Odyssey. Conference Proceedings (IEEE Cat. No. 01CH37295); IEEE: Piscataway, NJ, USA, 2001; pp. 661–666. [Google Scholar]
Demir, M.E.; Çalışkan, M. Shock Analysis of an Antenna Structure Subjected to Underwater Explosions. 2015. Available online: http://etd.lib.metu.edu.tr/upload/12619172/index.pdf (accessed on 9 February 2020).
Kwon, S.H.; Na, S.M.; Flatau, A.B.; Choi, H.J. Fe–Ga alloy based magnetorheological fluid and its viscoelastic characteristics. J. Ind. Eng. Chem. 2020, 82, 433–438. [Google Scholar] [CrossRef]
Wu, J.; Hu, H.; Li, Q.; Wang, S.; Liang, J. Simulation and experimental investigation of a multi-pole multi-layer magnetorheological brake with superimposed magnetic fields. Mechatronics 2020, 65, 102314. [Google Scholar] [CrossRef]
Bica, I.; Anitas, E.M.; Chirigiu, L.; Daniela, C.; Chirigiu, L.M.E. Hybrid magnetorheological suspension: Effects of magnetic field on the relative dielectric permittivity and viscosity. Colloid Polym. Sci. 2018, 296, 1373–1378. [Google Scholar] [CrossRef]
Bica, I.; Anitas, E.M.; Lu, Q.; Choi, H.J. Effect of magnetic field intensity and γ-Fe2O3 nanoparticle additive on electrical conductivity and viscosity of magnetorheological carbonyl iron suspension-based membranes. Smart Mater. Struct. 2018, 27, 095021. [Google Scholar] [CrossRef]

Figure 1. MR fluid: (a) Phenomenological characteristics; (b) Yield stress results; (c) Shape of particle.

Figure 2. Shape of the double pulse waveform profile.

Figure 3. Schematic diagram of the horizontal hydraulic double pulse shock test setup.

Figure 4. Photograph of the double pulse shock test setup: (a) Manufacture MR damper; (b) Experimental apparatus.

Figure 5. Measured damping force in the time and frequency domains: (a) Measured damping force in the time domain; (b) Measured damping force in the frequency domain.

Figure 6. MR damper experimental results at an impact velocity of 2.5 m/s (Test 1): (a) Damping force Piston velocity; (b) Piston Velocity; (c) F-V diagram of MR damper.

Figure 7. MR damper experimental results at an impact velocity of 5.37 m/s (Test 2): (a) Damping force Piston velocity; (b) Piston Velocity; (c) F-V diagram of MR damper.

Figure 8. MR damper experimental results at an impact velocity of 7.75 m/s (Test 3): (a) Damping force Piston velocity; (b) Piston Velocity; (c) F-V diagram of MR damper.

Figure 9. MR damper experimental results at an impact velocity of 9.12 m/s (Test 4): (a) Damping force; (b) Piston velocity; (c) force-velocity (F-V) diagram of MR damper.

Figure 10. F-V curve (Test. 4).

Figure 11. Estimation results of F-V model for MR damper: (a) Comparison results between experimental and simulator results; (b) Simulation F-V curve according to continuous current inputs.

Figure 12. GUI of the waveform prediction simulator program: (a)Flow chart; (b) GUI program.

Figure 13. Simulation result of the waveform prediction program.

Figure 14. Experimental result of the waveform prediction program.

Table 1. Measured experimental results.

Test No.	Impact Force at Test Table	Damping Force		Max. Velocity of Piston	Dynamic Range of Damping Force
Test No.	Impact Force at Test Table	Max.	Min.	Max. Velocity of Piston	Dynamic Range of Damping Force
1	108.208 kN	84.0 kN	48.9 kN	1.5 m/s	1.72
2	183.708 kN	144.0 kN	97.9 kN	2.75 m/s	1.47
3	246.678 kN	177.7 kN	132.5 kN	4.2 m/s	1.341
4	301.194 kN	201.3 kN	150.7 kN	4.4 m/s	1.336

Table 2. Values of the coefficients for the estimation of the F-V model with different inputs.

Vel. [m/s]	$c (i)$	$d (i)$
~0.0–1.0	−0.4534 i ² + 3.7204 i + 1.9272	0
~1.0–1.8	0.0773 i ² – 0.6573 i + 3.3945	−0.4534 i ² + 3.7204 i + 1.9272
~1.8–2.0	−0.0045 i ² + 0.1255 i + 3.3809	−0.3916 i ² + 3.1946 i + 4.6428
~2.0–2.5	−0.0586 i ³ + 0.5758 i2 – 1.4175 i + 1.9272	−0.3925 i ² + 3.1675 i + 5.319
~2.5–3.25	−0.187 i ² + 1.3076 i + 0.9671	−0.2724 i ² + 2.6697 i +6.9998
~3.25–3.5	−0.4325 i ² + 1.813 i + 2.8962	−0.4127 i ² + 3.6504 i + 7.7252
~3.5–4.5	0.1386 i ² – 1.3686 i + 6.9173	−0.5208 i ² + 4.1037 i + 8.4492
>4.5	2	−0.3822 i ² + 2.7351 i + 15.366

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Oh, J.-S.; Shul, C.W.; Kim, T.H.; Lee, T.-H.; Son, S.-W.; Choi, S.-B. Dynamic Analysis of Sphere-Like Iron Particles Based Magnetorheological Damper for Waveform-Generating Test System. Int. J. Mol. Sci. 2020, 21, 1149. https://doi.org/10.3390/ijms21031149

AMA Style

Oh J-S, Shul CW, Kim TH, Lee T-H, Son S-W, Choi S-B. Dynamic Analysis of Sphere-Like Iron Particles Based Magnetorheological Damper for Waveform-Generating Test System. International Journal of Molecular Sciences. 2020; 21(3):1149. https://doi.org/10.3390/ijms21031149

Chicago/Turabian Style

Oh, Jong-Seok, Chang Won Shul, Tae Hyeong Kim, Tae-Hoon Lee, Sung-Wan Son, and Seung-Bok Choi. 2020. "Dynamic Analysis of Sphere-Like Iron Particles Based Magnetorheological Damper for Waveform-Generating Test System" International Journal of Molecular Sciences 21, no. 3: 1149. https://doi.org/10.3390/ijms21031149

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu