Experimental study and load control of concentric electromagnetic damper under intensive impact load

Abstract

Background:

Conventional electromagnetic dampers impose a bias torque on artillery systems, which adversely affects system stability and performance.

Objective:

This study aims to propose a novel concentric electromagnetic damper (C-EMD). Based on the principle of electromagnetic induction, it achieves efficient braking to effectively eliminate the aforementioned bias torque effect.

Methods:

First, a nonlinear dynamic model of the C-EMD was established through theoretical analysis, finite element simulations, and impact experiments. This model was used to analyze the mechanism by which magnetic flux density affects damping performance and the parameter influence patterns. Subsequently, parameter sensitivity analysis was conducted using Optimal Latin Hypercube Design (OLHD) and polynomial regression methods. Finally, a surrogate model based on elliptic basis neural networks was developed, and multi-objective optimization was performed by integrating multi-island genetic algorithms with nonlinear programming by quadratic lagrangian.

Results:

Modeling analysis revealed that the resultant resistance curve exhibits a saddle-shaped profile with a concave center and elevated ends. Parameter sensitivity analysis results showed that liner thickness and air gap thickness have the most significant impact on performance (contributions >26%), while the influence of outer tube thickness is minimal (only 0.83%). Post-optimization results demonstrated: a reduction in resultant resistance curve flatness to 1.42, an increase in fill ratio to 89.86%, a 4.63% decrease in maximum displacement, and near-complete elimination of the saddle-shaped characteristics.

Conclusions:

This research provides an innovative solution for enhancing electromagnetic damper performance, with significant practical implications.

Keywords

Concentric electromagnetic damper intensive impact load experimental study sensitivity analysis load control

1. Introduction

An electromagnetic damper is an energy-dissipation device designed based on electromagnetic induction principles, with permanent magnets and conductors constituting its core components. Compared to viscous dampers, it eliminates concerns about sealing leakage and fluid replenishment in working media while offering minimal frictional wear, leading to broad applications in vehicle vibration reduction^1,2 and building seismic resistance.^3–5 In recent years, the emergence of high-performance NdFeB materials has significantly enhanced its energy dissipation density. Furthermore, electromagnetic dampers demonstrate advantages including non-contact operation, structural simplicity, and high reliability, driving increasing research in impact load scenarios such as weapon launching,^6,7 alongside rapid developments in modeling methodologies.

To address different configurations and application scenarios of permanent magnets, researchers have developed various eddy-current modeling theories. For example, Park et al.⁸ established a two-dimensional analytical model for Halbach-array eddy-current brakes using magnetic vector potential theory, validating the accuracy of braking torque calculations through finite element simulations and rotational speed-torque experiments, while analyzing critical parameter effects on braking performance. Min et al.⁹ employed a vector potential polar coordinate method to theoretically resolve eddy-current losses in Halbach-array magnetic couplers, combined with finite element verification to investigate segmentation effects on loss characteristics. Tian et al.¹⁰ proposed an efficient magneto-thermal coupling model for eddy-current brakes, simplifying electromagnetic field computations into mapping relationships between braking torque and stator skin-layer temperature, thereby significantly reducing computational costs. Their study revealed that optimizing thermal properties of stator materials enhances braking torque.

Research on artillery applications of electromagnetic dampers primarily employs theoretical analysis, numerical simulation, and experimental validation. Li et al.¹¹ investigated a high-energy density electromagnetic buffer (EMB). Based on design requirements, magnetic circuit parameters were determined. The equivalent current model was employed to analyze eddy current characteristics under demagnetization effects, while a segmented conduit strategy was adopted to resolve magnetization processes. Experimental validation confirmed the finite element results, demonstrating that the proposed EMB effectively accomplishes impact buffering. Li et al.¹² developed a novel permanent magnet linear eddy-current damper, constructing an improved equivalent magnetic circuit model based on magnetic field modulation theory. This model notably improved prediction accuracy for air-gap magnetic flux density and damping forces by comprehensively considering magnetic saturation effects, skin depth, and end effects. Li et al.¹³ introduced an equivalent subdomain method combining subdomain analysis and equivalent magnetic circuit modeling to calculate performance parameters of permanent magnet eddy-current dampers, establishing a prototype test platform to validate method effectiveness through experiments and finite element analysis. Shen et al.¹⁴ investigated the influence of different magnetic circuit configurations on braking performance, proposing a concentric recoil system design integrating eddy-current dampers with gun barrels, and analyzing parameter effects on recoil control characteristics. Liu et al.¹⁵ utilized equivalent subdomain models to study magnetic field behavior in permanent magnet eddy-current dampers under intense impact loads, developed time-varying finite element models to simulate operational processes, and designed a scaled prototype measurement system for dynamic magnetic field characterization under intensive impact conditions.

The Jiles-Atherton (J-A) model is a classical framework for characterizing hysteresis properties of ferromagnetic materials, capable of simultaneously describing macroscopic magnetization behavior and microscopic energy loss mechanisms. Research based on this theory and its improved variants has become an important approach for analyzing the nonlinear hysteresis characteristics and loss mechanisms of ferromagnetic materials. To elucidate demagnetization and damping mechanisms of permanent magnets under intensive impact loads and elevated temperatures, Li et al.¹⁶ proposed a modified exponential demagnetization model to simulate the demagnetization process of NdFeB permanent magnets, coupled with an improved vector hysteresis model accounting for hysteresis effects in soft magnetic materials. Li et al.¹⁷ systematically investigated demagnetization models of NdFeB permanent magnets and their application in eddy-current dampers, obtaining J-A model demagnetization curves for NdFeB at varying temperatures through experimental data fitting. Peng et al.¹⁸ developed a torque-adjustable bilateral axial-flux permanent magnet hysteresis damper, innovatively employing a modified vector J-A hysteresis model to precisely simulate flux distribution and hysteresis torque. Prototype validation demonstrated the model's high computational accuracy in characterizing hysteresis effects. Zheng et al.¹⁹ introduced an improved micromagnetic simulation method incorporating spatiotemporal correction coefficients for large-grain magnetic materials. Using sintered NdFeB as a case study, equivalent calculations based on thermal demagnetization theory yielded simulated results in excellent agreement with experimental thermal demagnetization curves.

Performance optimization of electromagnetic dampers in artillery systems is critical, with core objectives centered on precise control of barrel recoil and counter-recoil motion. Li et al.²⁰ proposed an analytical method for electromagnetic field calculation in tubular linear permanent magnet eddy current brake. By solving boundary value problems using Fourier series expansion and variable separation techniques, the approach was validated against finite element analysis. To enhance braking performance within spatial constraints, parametric analysis identified feasible design ranges, followed by multi-objective optimization to determine optimal parameters. The optimized design achieved a 27% increase in braking force while reducing permanent magnet volume by 7.7%. Li et al.²¹ presented a robust design methodology for electromagnetic dampers based on interval uncertainty optimization. Through sensitivity analysis using Optimal Latin Hypercube Design and polynomial regression, key parameters were screened to establish a resultant resistance force (RRF) considering robustness interval uncertainty optimization model. A radial basis function neural network-based interval nested optimization algorithm was developed to efficiently solve Pareto frontiers, with numerical simulations confirming the effectiveness of optimization objectives and strategies. Xu et al.²² addressed peak recoil resistance optimization by proposing a segmented optimization scheme considering parameter uncertainties via robust game theory. A novel robust Nash equilibrium solution method was developed through affine arithmetic and Chebyshev polynomial-based deep integration, mitigating electromagnetic damping force fluctuations caused by uncertainties. Optimized inner/outer tube thicknesses significantly reduced peak recoil resistance. Liu et al.²³ investigated dynamic characteristics and structural optimization of linear eddy-current dampers under artillery firing loads, elevating recoil resistance curve fill ratio from 84.93% to 92.1%.

In artillery systems, conventional electromagnetic dampers are typically installed in a non-coaxial configuration relative to the gun barrel. The resulting eccentric moment deteriorates the launch mechanical environment, increases muzzle disturbance, and consequently affects firing stability. The typical structural arrangement and positioning of such eccentric electromagnetic dampers within a full artillery system can be referred to in Fig. 2 of Reference⁷ and Fig. 22 of Reference.²⁴ Xie et al.²⁵ comparatively studied the influence of concentric versus eccentric damper configurations on muzzle vibration. Based on an established artillery launch dynamics model, their simulation analysis demonstrated that a concentric layout effectively reduces muzzle angular displacement and angular velocity, thereby helping to mitigate muzzle disturbance and improve firing accuracy. Building on this analysis, this study proposes a novel concentric electromagnetic damper. This concentric design not only improves space utilization, allowing for a more compact lateral layout, but also ensures that loads are transmitted along the axis. This effectively eliminates the eccentric moment caused by an offset configuration, thereby reducing muzzle disturbance during firing. Furthermore, this paper establishes a nonlinear dynamic model of the proposed damper. An intensive impact experimental platform was constructed for prototype validation, complemented by transient electromagnetic finite element simulations to analyze damping characteristics and parametric influence mechanisms. Combined optimization algorithms were implemented to control dynamic responses and enhance damping performance under intense impact loads, providing theoretical foundations and practical solutions for optimized electromagnetic damper design.

2. Nonlinear dynamics modeling for C-EMD

2.1. Principle and configuration of the C-EMD

The two-dimensional structure of the concentric electromagnetic damper (C-EMD) adopted in this study is illustrated in Fig. 1, featuring a compact design where the permanent magnet assembly is coaxially fixed to the moving rod. The device comprises two primary components: a mover and a stator. The mover consists of a moving rod, liner, permanent magnets, and iron poles, while the stator is formed by a composite structure with an inner conductive tube and outer ferromagnetic tube. The inner tube utilizes 6061 aluminum alloy for its superior machinability, where high electrical conductivity generates eddy currents that impose damping forces on the mover to achieve energy dissipation. The outer tube employs 45 steel, whose ferromagnetic properties enhance magnetic field utilization and improve damping performance. Detailed structural parameters are provided in Table 1.

Figure 1.

Schematic diagram of the C-EMD.

Table 1.

Material and structural parameters of the C-EMD.

Part	Material	Inner radius/mm	Outer radius/mm	Length/mm	Number
The moving rod	42CrMo	60.7	72	1400	1
Liner	6061 alloy	72	80	340	1
Permanent magnet	Nd-Fe-B 52	80	120	16	12
Iron pole	DT4C	80	122	6	13
Inner tube	6061 alloy	123	126	600	1
Outer tube	45 steel	126	150	600	1

The C-EMD operates based on Faraday's law of electromagnetic induction and Lenz's law: Relative motion between the mover and stator induces eddy currents within the inner tube. These eddy currents generate opposing magnetic fields that interact with the original field, producing damping forces counteracting the motion direction. This mechanism converts transient intensive impact loads into prolonged and smoothed loading profiles, stabilizing the recoil and counter-recoil process of the mover and enhancing the resistance plateau effect. Furthermore, as shown in Fig. 2, adjacent permanent magnets are arranged in an opposing-polarity configuration to enhance magnetic field utilization and damping performance. Their length is $τ_{m}$ and radius is R.

Figure 2.

Schematic diagram of cylindrical permanent magnets in homopolar opposing configuration.

2.2. Finite element model under intensive impact load

Electromagnetic field modeling of the C-EMD shown in Fig. 1 is performed based on the differential form of Maxwell's equations:

rot H = J + \frac{\partial D}{\partial t}

(1)

rot E = - \frac{\partial B}{\partial t}

(2)

d i v D = ρ

(3)

d i v B = 0

(4)

where

H

is the magnetic field intensity,

J

the current density,

D

the electric displacement vector,

E

the electric field intensity,

B

the magnetic flux density, and

ρ

the charge density.

Additionally, the constitutive relations governing the physical quantities in the medium are expressed as:

J = σ (E + v \times B)

(5)

D = ε E

(6)

B = μ H

(7)

where

σ

is the electrical conductivity,

v

the relative motion velocity,

ε

the permittivity, and

μ

the permeability.

The magnetic insulation boundary conditions are imposed as:

n \times A = 0

(8)

where

n

is the unit normal vector and

A

denotes the magnetic vector potential.

Based on the Lorentz force equation, the electromagnetic damping force generated by the inner and outer tubes cutting magnetic flux lines can be expressed as:

F_{n} = \int_{V_{i n}} (J \times B) d V_{i n} + \int_{V_{o u t}} (J \times B) d V_{o u t}

(9)

where

V_{i n}

and

V_{o u t}

represent the volumes of the inner and outer tubes, respectively.

Jiles and Atherton^26,27 proposed the hysteresis Jiles-Atherton model for ferromagnetic materials based on energy variations induced by pinning effects:

\frac{d M}{d H} = (1 - c) \frac{M_{a n} - M_{i r r}}{k δ - α (M_{a n} - M_{i r r})} + c \frac{d M_{a n}}{d H}

(10)

Once the five key parameters—saturation magnetization $M_{s}$ , domain wall density a, pinning loss k, magnetic reversibility c, and interdomain coupling $α$ —are determined, the hysteresis loop of the ferromagnetic material can be obtained by solving Equation (10). The B-H characteristic satisfies:

B = μ_{0} (M + H)

(11)

where

μ_{0}

is the vacuum permeability, and the total magnetization M comprises reversible (

M_{r e v}

) and irreversible (

M_{i r r}

) components:

M = M_{r e v} + M_{i r r}

(12)

Under a given magnetic field intensity, the reversible component reduces the difference between the irreversible magnetization and the anhysteretic magnetization $M_{a n}$ , expressed as:

M_{r e v} = c (M_{a n} - M_{i r r})

(13)

Based on ferromagnetic material magnetization theory, the magnetization behavior can be described by the Langevin function:

M_{a n} (H_{e}) = M_{s} [c o t h (\frac{H_{e}}{a}) - \frac{a}{H_{e}}]

(14)

where

H_{e}

is the effective magnetic field intensity, composed of the external field H and the molecular field

H_{m}

H_{e} = H + H_{m} = H + α M

(15)

In the Jiles-Atherton hysteresis model, the reversible and irreversible magnetization components are governed by the following differential equations:

\frac{d M_{r e v}}{d H} = c (\frac{d M_{a n}}{d H} - \frac{d M_{i r r}}{d H})

(16)

\frac{d M_{i r r}}{d H} = \frac{M_{a n} - M_{i r r}}{k δ - α (M_{a n} - M_{i r r})}

(17)

where

δ

indicates the direction of magnetic field change:

δ = + 1

for

d H / d t > 0

, and

δ = - 1

for

d H / d t < 0

Based on the theoretical framework described above and neglecting the influence of temperature variations, mechanical vibrations, and other external factors, this study employs the finite element analysis software COMSOL Multiphysics^® to simulate the dynamic response of the C-EMD. Given the rotationally axisymmetric geometry of the damper, a two-dimensional (2D) axisymmetric electromagnetic finite element model (Fig. 3) was developed to balance computational accuracy and efficiency. Components with negligible impact on the magnetic field were omitted in the model. Material properties and structural parameters in the simulation strictly adhere to the configurations listed in Table 1, ensuring consistency with experimental conditions.

Figure 3.

Finite element model of the C-EMD.

In the model, material properties and parameters were assigned to each domain, followed by the addition of relevant physics interfaces. The material properties in the simulation model are shown in Table 2, with a relative permittivity of 1 for all materials. Magnetic insulation boundary conditions were imposed, and the permanent magnet characteristics were defined using the hysteresis Jiles-Atherton model with parameters referenced from literature.¹⁷ The temperature was 293.15 K and the initial magnetization was 1.15 × 10⁶A/m, as specified in Table 3. The homopolar opposing magnetic circuit design was implemented by specifying the positive/negative of the initial magnetization intensity. To accurately characterize the electromagnetic damping force, the magnetic saturation characteristics of the iron pole is described using the curve shown in Fig. 4. This approach, which neglects magnetic hysteresis while accounting for nonlinear characteristics, offers a favorable trade-off: it incurs a minimal loss in computational accuracy while significantly enhancing computational efficiency. Experimentally measured intensive impact loading profiles were imported as input variables, motion process variables were defined, and the global ordinary differential equations (ODEs) of motion were solved via the moving mesh method. Triangular meshing was adopted, with local mesh refinement applied to the air gap and inner tube regions to mitigate skin effect-induced errors. The total mesh count was 1.32 × 10⁵ elements with an average element quality of 0.91. Steady-state simulation results served as initial conditions for transient analysis. A time step of 5 × 10⁻⁵ s was set for the 0.035 s simulation duration. Solver configurations were optimized for transient analysis, and numerical computations yielded electromagnetic field and dynamic response results under intense impact loading.

Figure 4.

B–H curve of the iron pole.

Table 2.

Material parameters of 42CrMo, 6061 alloy, Nd-Fe-B 52, DT4C, and 45 steel.

Material properties	Symbol	42CrMo	6061 alloy	Nd-Fe-B 52	DT4C	45 steel
Conductivity	$σ$ /(S/m)	4.61e6	2.43e7	7.14e5	2e6	5.78e6
Density	$ρ$ /(kg/m³)	7850	2700	7400	7870	7850

Table 3.

Material parameters of N52.

J-A model	$M_{s}$ /(A/m)	$a$ /(A/m)	$k$ /(A/m)	$c$	$α$
Nd-Fe-B 52	1.3e6	1.02e5	1.31e6	2e-4	0.84

During the motion of the C-EMD, the system is subjected to electromagnetic damping forces, counter-recoil forces, and frictional forces. These three forces collectively constitute the resultant resistance $F_{R}$ , expressed as:

F_{R} = F_{n} + F_{f} + F_{T}

(18)

where

F_{f}

is the counter-recoil force and

F_{T}

the frictional force,

F_{T} = 80 N

The counter-recoil force, generated by a spring, acts as a restoring force to return the mover to its initial position after the damping stroke. It exhibits a linear relationship with displacement, as shown in Fig. 5, and is calculated by:

F_{f} = F_{f 0} + K x

(19)

where

F_{f 0}

is the initial counter-recoil force (ensuring the mover overcomes friction to return to the initial position even under minor displacements)

F_{f 0} = 750 N

, K the spring stiffness coefficient

K = 2.8 \times 10^{4} N / m

, and x the displacement.

Figure 5.

Counter-Recoil force.

The governing differential equation for the damper's motion is:

m_{h} \frac{d^{2} x}{d t^{2}} = m_{h} \frac{d v}{d t} = F_{p t} - F_{R}

(20)

where

m_{h}

is the mover mass, and

F_{p t}

denotes the intense impact load (the primary driving force of the mover's motion), with its time-history data obtained experimentally.

3. Impact tests

3.1. Concentric electromagnetic damper impact tests

Experiments were conducted using the C-EMD test bench shown in Fig. 6, which primarily consists of a projectile, composite tube, permanent magnet assembly, moving rod, frame, and spring. The permanent magnet assembly, spring, and moving rod are coaxially arranged, with this concentric configuration eliminating bias torque effects on the artillery system. The composite tube and frame form the stator, fixed to the test bench, while the moving rod, spring, and permanent magnet assembly constitute the mover. Under intensive impact loading, the projectile is launched along the moving rod, inducing relative motion between the mover and stator. This generates eddy currents in the composite tube, producing damping forces to counteract the mover's motion and achieve braking. The spring compresses to store energy during recoil and restores its original length during counter-recoil, providing restoring forces to reset the mover and complete the recoil and counter-recoil cycle.

Figure 6.

Test bench platform for the C-EMD.

The dynamic test system for the C-EMD (Fig. 7) includes a displacement sensor, acceleration sensor, force sensor, charge amplifier, and data acquisition system. The displacement sensor is fixed directly behind the damper, while the acceleration sensor and force sensor are mounted on the breech. Upon ignition trigger initiation, the mover completes one recoil and counter-recoil cycle under impact loading. The displacement sensor outputs voltage signals directly to the data acquisition system, whereas charge signals from the acceleration sensor and force sensor are converted to voltage signals via the charge amplifier before being recorded, with a sampling frequency of 200 kHz. The data acquisition system processes these signals into displacement, acceleration, and pressure data, with velocity derived by differentiating displacement. Specifications of the sensors are listed in Table 4.

Figure 7.

Dynamic testing system of the C-EMD.

Table 4.

Type and performance parameters of the sensors.

Sensor name	Type	Sensitivity	Measuring range
Displacement sensor	OptoNCDT 1402	6820mm/V	200～800mm
Acceleration sensor	B&K type4371V	9.315pC/g	±20000g
Force sensor	Kistler 6213B	−1.205pC/bar	0～8000bar

Test results for three experimental trials are summarized in Table 5, revealing significant differences in mover dynamics under varying peak impact loads. With an increase in the charge mass, the maximum pressure rises correspondingly, which consequently and directly leads to synchronous increases in both the maximum displacement and maximum acceleration of the mover. The recoil and counter-recoil process for Experiment 1 (Fig. 8) demonstrates complete operational reliability: the mover attains a maximum displacement of 137.67 mm within 31.43 ms before fully resetting. Smoke dispersion is visibly observed during reset (Fig. 8(d)). Measured impact load curves for all three trials (Fig. 9) exhibit peak values of 97.32 kN, 127.09 kN, and 159.2 kN, respectively, characteristic of high-transient, intense impact loading.

Figure 8.

Mover recoil and counter-recoil process.

Figure 9.

Intensive impact loads.

Table 5.

Experimental test results.

	Maximum pressure/bar	Maximum displacement/mm	Maximum acceleration/g
Test 1	346.99	137.67	217.19
Test 2	443.58	168.54	275.87
Test 3	564.04	215.15	367.73

3.2. Verification of electromagnetic finite element simulation results

As shown in Fig. 10, the displacement, velocity, and acceleration curves of the concentric electromagnetic damper's nonlinear dynamic model under three distinct impact loads exhibit excellent agreement with experimental measurements. Key metrics—including maximum displacement, maximum velocity, maximum acceleration, recoil completion time, and velocity peak timing—align closely with finite element method (FEM) results, with minimal errors (Table 6). This validates the simulation model's capability to accurately characterize the dynamic response of the C-EMD. Furthermore, the model and simulation methodology can be employed to analyze the influence of varying impact loads on damping performance and optimize the structural design of the C-EMD.

Figure 10.

Simulation and experiment comparison curves (a) displacement - time, (b) velocity - time, (c) acceleration – time.

Table 6.

Comparative analysis of numerical simulation and impact test.

Comparison item		Impact 1	Impact 2	Impact 3
Maximum displacement/mm	FE value	138.82	170.85	212.39
	test value	137.67	168.54	215.15
	error	0.84%	1.35%	1.28%
Maximum velocity/(m·s^-1)	FE value	8.45	10.52	12.78
	test value	8.37	10.22	12.84
	error	0.96%	2.94%	0.47%
Maximum acceleration/g	FE value	208.98	269.98	352.54
	test value	217.19	275.87	367.73
	error	3.78%	2.13%	4.13%
Recoil completion time/ms	FE value	31.80	32.1	32.8
	test value	31.43	32.85	31.90
	error	1.18%	2.28%	2.82%
Velocity peak time/ms	FE value	7.95	7.90	7.85
	test value	8.14	7.06	8.33
	error	2.23%	11.90%	5.76%

3.3. Mechanism of magnetic flux density effects on damping performance

The damping coefficient, a critical parameter for evaluating the damping performance of electromagnetic dampers, quantifies the braking effectiveness and can be expressed as²⁸:

C = σ_{i n} \int_{V_{i n}} B_{r}^{2} d V_{i n}

(21)

where

σ_{i n}

denotes the electrical conductivity of the inner tube. Equation (21) indicates that the damping coefficient is directly governed by the electrical conductivity of the inner tube and the spatial magnetic flux density distribution. Since the inner tube's conductivity remains nearly constant under varying impact loads, the damping coefficient is predominantly determined by the magnetic flux density distribution.

Selecting the inner tube of the C-EMD as the integration domain, the calculated C values under three distinct impact loads are illustrated in Fig. 11. At $t = 0$ , the damper is stationary, and the normalized damping coefficient C reaches its maximum value in all three cases, initiating the braking process. Subsequently, $C$ decreases sharply compared to the static magnetic field condition at rest. As braking progresses and the mover velocity diminishes, the C curve gradually recovers. This behavior stems from the interaction between the induced reverse magnetic field (opposing the original field) and the primary field, collectively termed the demagnetization effect. The temporal evolution of C reflects the dynamic demagnetization effects during damper operation. Divergences among the three C curves indicate varying demagnetization intensities under different impact loads, with Impact Load Case 3 exhibiting the most pronounced eddy current-induced demagnetization. Larger impact loads result in higher relative velocities between the inner tube and magnetic field, generating stronger eddy currents and consequently amplifying the demagnetization effect. As shown in Fig. 12, the damping coefficient varies with velocity: When $v < 2.5 m / s$ , the damping coefficient remains approximately constant at $c = 4.05 kN \cdot (m \cdot s^{- 1})^{- 1}$ ; when $v > 2.5 m / s$ , it gradually decreases with increasing mover velocity. In the low-speed region ( $v < 2.5 m / s$ ), the demagnetization effect is weak, resulting in a constant damping coefficient. When entering the high-speed region ( $v > 2.5 m / s$ ), the enhanced demagnetization effect causes the damping coefficient to decline.

Figure 11.

Damping coefficient-time curve.

Figure 12.

Damping coefficient-velocity curve.

A region spanning four times the pole pitch at the median diameter position of the inner tube was selected to analyze magnetic field variations under different impact loads. Figure 13 illustrates the temporal evolution of radial magnetic flux density at this location. Influenced by the permanent magnet arrangement and the flux-conducting properties of the iron poles, a periodic traveling-wave magnetic field is formed, with maximum radial magnetic flux density occurring at the iron poles and opposite magnetic flux density directions observed between adjacent iron poles. At t = 0, the magnetic field stabilizes, and the magnetic flux density curve exhibits an approximately symmetric distribution about the center of the iron poles, peaking at 1.74 T. At $t = 7.9 ms$ , when braking velocities reach peak values under all three impact loads, edge effects induce symmetry breakdown in the magnetic flux density profile, deviating significantly from the initial state. Impact Load Case 3 yields the smallest magnetic flux density peak (1.31 T), indicating the most pronounced demagnetization effect. While all three impact loads induce comparable magnetic flux density offsets, their demagnetization intensities differ markedly, with higher loads amplifying eddy current-driven field weakening.

Figure 13.

Radial magnetic flux density at the median diameter of the inner tube.

4. Sensitivity analysis of the C-EMD system responses

4.1. Analysis of parameters affecting damping performance

The influence of structural parameters on the damping performance of the C-EMD system manifests in two primary aspects: the parametric dependence of electromagnetic damping forces and the spring stiffness effects on counter-recoil forces. Under constant permanent magnet parameters, a single-variable control method was employed to investigate seven independent parameters—liner thickness, axial length of iron pole, iron pole wall thickness, air gap thickness, inner tube thickness, outer tube thickness, and spring stiffness—on damping characteristics. Given the high velocity of the mover during recoil, these parameters significantly affect dynamic behaviors, making the resultant resistance force (RRF) during recoil the focal metric for analysis. To ensure controlled variables, experimentally validated impact load data from Test Case 1 (with minimal error) was adopted as the input.

During recoil, the C-EMD is subjected to combined electromagnetic damping, counter-recoil, and frictional forces, collectively forming the saddle-shaped resultant resistance force shown in Fig. 14. Analysis of Test Case 1 reveals two distinct resistance peaks ( $F_{e_m a x}$ and $F_{f_m a x}$ ) with an intermediate valley ( $F_{v_m i n}$ ). The resultant resistance force reaches its initial peak $F_{e_m a x}$ with increasing displacement, where eddy currents in the inner tube intensify, amplifying demagnetization effects and causing the resultant resistance force to drop to valley $F_{v_m i n}$ . The secondary peak $F_{f_m a x}$ emerges from progressive counter-recoil force accumulation with displacement, reaching maximum displacement $λ_{m a x}$ when velocity decays to zero.

Figure 14.

Schematic diagram of the RRF curve.

The influence of seven independent parameters on the resultant resistance force versus recoil displacement characteristics is illustrated in Fig. 15, with subplots (a)-(g) corresponding to parametric analyses of liner thickness ( $x_{1}$ ), axial length of iron pole ( $x_{2}$ ), iron pole wall thickness ( $x_{3}$ ), air gap thickness ( $x_{4}$ ), inner tube thickness ( $x_{5}$ ), outer tube thickness ( $x_{6}$ ), and spring stiffness ( $x_{7}$ ), respectively.

Figure 15.

Influence of different parameters on the resultant resistance force-recoil travel curve.

As a magnetic isolation layer between the permanent magnet assembly and moving rod, increased liner thickness reduces magnetic flux leakage through the rod, thereby enhancing magnetic field utilization. This results in elevated RRF peaks/valleys and reduced maximum displacement (Fig. 15(a)). The initial peak $F_{e_m a x}$ grows with increasing axial iron pole length ( $x_{2}$ ), albeit with diminishing returns due to structural and material saturation effects on flux-carrying capacity. Minimum maximum displacement occurs at $x_{2} = 8 mm$ (Fig. 15(b)). Thicker iron pole walls ( $x_{3}$ ) increase outer diameters, raising magnetic circuit reluctance and reducing air gap magnetic flux density. Consequently, RRF peaks/valleys decrease while maximum displacement increases (Fig. 15(c)). Air gap expansion ( $x_{4}$ from 0.5 mm to 2.0 mm) elevates magnetic circuit reluctance, weakening the primary field and lowering RRF peaks/valleys while extending maximum displacement (Fig. 15(d)). For inner tube thicknesses ( $x_{5}$ ) between 1.0∼1.5 mm, eddy currents remain unsaturated. Thicker inner tubes enhance $F_{e_m a x}$ and reduce maximum displacement. Beyond 1.5 mm, skin effects, demagnetization from eddy-induced fields, and increased magnetic resistance collectively diminish peaks/valleys while prolonging displacement (Fig. 15(e)). The outer tube ( $x_{6}$ ) structurally supports the damper and confines magnetic flux. Increased thickness enhances flux closure within the composite tube, marginally boosting $F_{e_m a x}$ but exerting limited influence on RRF due to peripheral positioning in the main magnetic circuit (Fig. 15(f)). Higher spring stiffness ( $x_{7}$ ) amplifies counter-recoil forces, disproportionately increasing the secondary peak $F_{f_m a x}$ compared to $F_{e_m a x}$ and $F_{v_m i n}$ , while reducing maximum displacement (Fig. 15 g).

4.2. The sensitivity analysis based on OLHD and polynomial regression analysis

The sensitivity analysis of the C-EMD system's dynamic response is performed using a combined experimental design (DOE) and polynomial regression analysis approach, with the following procedural steps:

Step 1: Parameter space sampling. Parameter vectors $X = (X_{1} X_{2} \dots X_{k} \dots X_{n})$ are generated based on the number of design variables, constrained within their respective ranges $[X_{i}^{\min}, X_{i}^{\max}]$ . To ensure uniformity of projection space in the design matrix, an Optimal Latin Hypercube Design (OLHD) method is adopted. OLHD improves upon traditional Latin Hypercube Design (LHD) by enhancing spatial uniformity, thereby achieving higher accuracy in factor-response fitting. A schematic of the OLHD sampling strategy is illustrated in Fig. 16.

Figure 16.

Sampling schematic of the OLHD method.

Step 2: Parametric modeling. A nonlinear dynamic parametric model of the C-EMD is established. The parameter design matrix $D_{m \times n}$ is linked to the parametric model via a CSV-formatted file.

Step 3: Finite element simulation. The mesh is updated and a new electromagnetic finite element model is reconstructed based on the parameter design matrix $D_{m \times n}$ . Automated electromagnetic field simulations are performed for all 100 samples.

Step 4: Sensitivity analysis. A linear polynomial regression model is employed for sensitivity analysis of the system response:

\begin{matrix} y_{i} = β_{0} + β_{1} χ_{1 i} + β_{2} χ_{2 i} + \dots + β_{k} χ_{k i} + ε, ε - N (0, σ^{2}) \end{matrix}

(22)

where

β_{k}

and

σ^{2}

are estimated via ordinary least squares (OLS), minimizing the residual sum of squares (RSS):

Q = \sum_{i = 1}^{n} (y_{i} - β_{0} - β_{1} χ_{1 i} - β_{2} χ_{2 i} - \dots - β_{k} χ_{k i})^{2}

(23)

The regression coefficients are obtained by solving partial derivatives of Equation (23) set to zero.

Step 5: Data normalization. Sample data are normalized to the [0, 1] interval interval using min-max scaling, eliminating dimensional and scale effects.

Step 6: Contribution ratio analysis. The influence of each parameter is quantified as a percentage contribution by scaling regression coefficients:

N_{i} = 100 β_{i} / \sum_{j}^{k} | β_{i} |

(24)

Step 7: Comprehensive sensitivity synthesis. Absolute values of sensitivity results for four key dynamic response metrics ( $F_{e_m a x}$ , $F_{v_m i n,} F_{f_m a x}$ , $λ_{m a x}$ ) are aggregated. The total contribution ratio of each parameter is calculated as:

χ_{i} = \frac{| N_{1 i} |}{n} + \frac{| N_{2 i} |}{n} + \dots + \frac{| N_{n i} |}{n}

(25)

Step 8: Dynamic response variance analysis. Analysis of variance (ANOVA) is conducted to evaluate statistical significance. Total variance comprises model-explained variance and residual error. The coefficient of determination ( $R^{2}$ ) quantifies goodness-of-fit:

R^{2} = 1 - \frac{\sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}}{\sum_{i = 1}^{n} {(y_{i} - {\bar{y}}_{i})}^{2}}

(26)

where

{\hat{y}}_{i}

is the predicted response for the i-th sample,

{\bar{y}}_{i}

the mean response, and n the sample size. Values closer to 1 indicate superior fit accuracy.

The sensitivity analysis results for each parameter's influence on $F_{e_m a x}$ , $F_{v_m i n,} F_{f_m a x}$ , $λ_{m a x}$ , and the total contribution ratio are illustrated in Fig. 17. The analysis reveals that liner thickness ( $x_{1}$ ) and air gap thickness ( $x_{4}$ ) exert the most significant impact on all four dynamic responses, while the outer tube thickness ( $x_{6}$ ) exhibits negligible influence. Statistical validation through Equation (26) yields $R^{2}$ -values of 0.992, 0.975, 0.978, and 0.981 for $F_{e_m a x}$ , $F_{v_m i n}$ , $F_{f_m a x}$ , $λ_{m a x}$ , respectively, indicating high goodness-of-fit and confirming the reliability of the sensitivity analysis. Liner thickness and air gap thickness dominate the total contribution ratios at 27.91% and 32.31%, respectively, whereas outer tube thickness contributes minimally (0.83%). Consequently, $x_{6}$ can be disregarded in subsequent optimization studies. The optimized parameters include $x_{1} \sim x_{5}$ and $x_{7}$ , with their initial values and bounds specified in Table 7.

Figure 17.

Sensitivity analysis results of the C-EMD response under intensive impact load.

Table 7.

Optimization parameters.

Parameters	Initial value	Lower bound	Upper bound
x₁/mm	4	2	6
x₂/mm	6	4	14
x₃/mm	21	20	22.5
x₄/mm	0.5	0.5	2
x₅/mm	1.5	1	3
x₇/(N/mm)	28	15	35

5. Load control

5.1. Surrogate model construction

During the optimization of the resultant resistance in the C-EMD, the high computational cost of repeatedly invoking the coupled simulation model—owing to its nonlinearity and complexity—poses significant challenges. Direct optimization using the original model would substantially increase computational complexity, making it impractical to achieve optimal results within acceptable resource constraints. To address this, surrogate models with higher computational efficiency are widely adopted as effective substitutes for complex coupled simulations. This study evaluated multiple surrogate models, including third- and fourth-order response surface models (RSMs), radial basis function (RBF) and elliptic basis function (EBF) neural networks, Chebyshev orthogonal polynomial models, and Kriging models. Using a design of experiments approach, 150 experimental samples were generated, randomized, and normalized, with 120 samples allocated to the training set and 30 to the testing set.

The root mean square error (RMSE) was selected as the metric to evaluate the predictive performance of the surrogate models. An RMSE closer to 0 indicates smaller deviations between predicted and actual values, reflecting higher model accuracy. The RMSE is calculated as:

R M S E = \sqrt{\frac{1}{n} \sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}}

(27)

As shown in Table 8, the RMSE values for the six models demonstrate that the elliptic basis function (EBF) neural network achieves the highest accuracy and optimal performance. Consequently, the EBF model was selected for subsequent multi-objective optimization.

Table 8.

Error analysis of surrogate model.

Surrogate model	$F_{e_m a x}$	$F_{v_m i n}$	$F_{f_m a x}$	$λ_{m a x}$
RSM3	0.026	0.036	0.036	0.028
RSM4	0.027	0.033	0.037	0.025
RBF	0.027	0.04	0.036	0.032
EBF	0.021	0.033	0.024	0.023
COPM	0.025	0.038	0.028	0.032
KM	0.076	0.073	0.063	0.082

The elliptic basis function (EBF) neural network is a three-layer feedforward architecture comprising input, hidden, and output layers. The hidden layer performs nonlinear transformations using the Mahalanobis distance between test points and sample centroids as inputs, expressed as:

g_{j} (x) = g (‖ x - x_{j}^{c} ‖_{m}) (j = 1, 2, \dots, N)

(28)

‖ x - x_{j} ‖_{m} = (x - x_{j})^{T} Q^{- 1} (x - x_{j})

(29)

Q = \frac{1}{N} \sum_{j = 1}^{N} (x_{j} - μ) (x_{j} - μ)^{T}

(30)

where x is the input vector, Q the covariance matrix, N the number of samples,

μ

the centroid of the j-th cluster, and

c

the basis function (

0.2 < c < 3

The mapping from the hidden layer to the output layer is linear, with response values calculated via weighted summation:

y_{i} = \sum_{j = 1}^{N} ω_{j} g_{j} (x_{i}) + ω_{N + 1}

(31)

\sum_{j = 1}^{N} ω_{j} = 0

(32)

where

ω

denotes the connection weights. The topology of the EBF model is illustrated in Fig. 18.

Figure 18.

Topology diagram of the EBF model.

5.2. Multi-variable optimization design based on combinatorial algorithm

Optimization parameters (Table 7) include electromagnetic variables (e.g., axial iron pole length) and mechanical properties (e.g., spring stiffness). The objectives are to minimize the initial peak ( $F_{e_m a x}$ ), secondary peak ( $F_{f_m a x}$ ), and maximum displacement ( $λ_{m a x}$ ) of the resultant resistance force, while maximizing the valley resistance ( $F_{v_m i n}$ ), with the corresponding objective functions being $f_{1} (x)$ , $f_{3} (x)$ , $f_{4} (x)$ , and $f_{2} (x)$ . A multi-parameter, multi-objective optimization model is formulated as:

{\begin{matrix} \min [f_{1} (x), f_{3} (x), f_{4} (x)] \\ \max f_{2} (x) \\ s . t . f_{4} (x) \leq 138 \\ X = (x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{7}) \\ x_{i}^{L} \leq x_{i} \leq x_{i}^{U}, i = 1, 2, 3, 4, 5, 7 \end{matrix}

(33)

where

x_{i}^{L}

and

x_{i}^{U}

define parameter bounds.

To address high-dimensional search challenges, a hybrid algorithm combining the Multi-Island Genetic Algorithm (MIGA) and Nonlinear Programming by Quadratic Lagrangian (NLPQL) was implemented. MIGA partitions populations into subpopulations (“islands”) for parallel evolution,²⁹ while NLPQL refines solutions locally using sequential quadratic programming. The optimization workflow is summarized in Fig. 19.

Figure 19.

Optimization flowchart for RRF.

5.3. Optimization results analysis

Based on Pareto multi-objective theory, an optimal parameter set (Table 9) was selected from the Pareto front, balancing engineering constraints and objective weights. Post-optimization dynamics (Fig. 20) demonstrate reduced saddle-shaped characteristics and eliminated peak protrusions, confirming enhanced damping stability. Key metrics (Table 10) show decreased initial/secondary peaks ( $F_{e_m a x}$ , $F_{f_m a x}$ ), increased valley resistance ( $F_{v_m i n}$ ), and reduced maximum displacement ( $λ_{m a x}$ ).

Figure 20.

Comparison of curves before and after optimization.

Table 9.

Optimization results of the combined algorithm.

Parameters	x ₁	x ₂	x ₃	x ₄	x ₅	x ₇
Results	6	4.47	20	0.5	2.65	23.36

Table 10.

Target values before and after optimization.

Optimization objectives	$F_{e_m a x}$	$F_{v_m i n}$	$F_{f_m a x}$	$λ_{m a x}$
Initial value a	29.40	28.31	29.38	138.82
Simulation-optimized value b	28.97	28.67	29.05	132.39
Algorithm-optimized value c	28.87	28.08	28.67	131.16
[(b-a)/a]/%	−1.46	1.27	−1.12	−4.63
[(c-b)/b]/%	−0.35	−2.06	−1.31	−0.93

To quantify curve smoothness, fill ratio ( $α$ ) and flatness ( $β$ ) are defined as:

α = S_{1} / S_{2} \times 100 %

(34)

β = \frac{1}{q} \sum_{ω = 1}^{q} k_{ω}, ω = 1, 2, \dots q

(35)

where

S_{1}

is the area under the resistance curve,

S_{2}

is the rectangular area spanned by maximum displacement and peak resultant resistance, q is the number of equidistant continuous sampling points over the entire recoil stroke, assigned as

q = 10000

, and

k_{ω}

is the derivative of the resultant resistance curve at the

ω

-th sampling point.

Under ideal conditions, the resultant resistance should maintain a constant value. Consequently, a plateau filling ratio ( $α$ ) approaching 1 and a diminished flatness ( $β$ ) indicate superior platform characteristics of the resultant resistance curve. Post-optimization results achieve $α = 89.86 %$ (vs. 86.98% pre-optimization) and $β = 1.42$ (vs. 1.46), with a 4.63% reduction in $λ_{m a x}$ , demonstrating improved plateau stability.

6. Conclusion

Based on the C-EMD structure, this study developed a nonlinear dynamic model and validated its accuracy through experiments and finite element simulations. The influence of load on damping performance was analyzed using three sets of impact experiments with different peak values, and the effects of magnetic flux density and structural parameters on damping performance were investigated. Key parameters were screened via sensitivity analysis. Subsequently, a multi-objective optimization was achieved by employing a surrogate model and a combinatorial optimization algorithm, which effectively controlled the dynamic response. The specific conclusions are as follows:

The principle and configuration of the C-EMD were introduced. Its nonlinear dynamic and electromagnetic FEM were established. An experimental platform was built to measure the mover's motion characteristics under different impact loads, obtaining data on displacement, velocity, and acceleration, thereby validating the accuracy of the FEM. The finite element analysis indicated that as the impact load increases, the relative speed between the inner tube and the original magnetic field accelerates, generating stronger eddy currents and a more pronounced demagnetization effect, thereby influencing the damping performance.

The influence patterns of structural parameters on damping performance were analyzed. Sensitivity analysis of the damper's structural parameters was conducted using OLHD and polynomial regression methods. The results showed that liner thickness and air gap thickness have the most significant impact on performance, contributing 27.91% and 32.21%, while the influence of outer tube thickness is minimal, at only 0.83%. By comparing the accuracy of six surrogate models, the EBF neural network was selected to construct the dynamic response model. After optimization using the MIGA combined with the NLPQL algorithm, the RRF curve characteristics improved significantly: the flatness decreased to 1.42, the fill ratio increased to 89.86%, the braking stability was enhanced, and the maximum displacement was reduced by 4.63%.

This study established an impact experimental platform for the C-EMD under intensive impact loads. The combinatorial optimization algorithm was successfully applied to control its dynamic response effectively, optimizing the braking performance. This work provides a new technical pathway for further enhancing the performance of EMD and holds significant value for promoting their technological development and engineering applications.

Footnotes

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The work was supported by the Postgraduate Research & Practice Innovation Program of Jiangsu Province under Grant KYCX24_0716.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Ebrahimi

Bolandhemmat

Khamesee

, et al. A hybrid electromagnetic shock absorber for active vehicle suspension systems. Veh Syst Dyn 2011; 49: 311–332.

Cao

Liu

, et al. Multi-Field coupling analysis and demagnetization experiment of permanent magnet retarder for heavy vehicles (MAY 2018). IEEE Access 2019; 7: 50734–50745.

Zhang

Chen

Hua

, et al. Design and dynamic characterization of a large-scale eddy current damper with enhanced performance for vibration control. Mech Syst Signal Process 2020; 145: 106879.

Feng

Zhang

Jing

, et al. Train-induced longitudinal motion in suspension bridges: synergetic control with eddy current damping and frictional damping. Structures 2025; 71: 108002.

Wang

Zhou

. Bi-directional adaptive eddy current pendulum tuned mass damper for structural nonlinear seismic response control. Journal of Building Engineering 2025; 99: 111455.

Wang

Yang

. Magnetomechanical degradation of sintered NdFeB induced by impact in eddy current recoil brake. J Magn Magn Mater 2022; 564: 170114.

Yang

Wang

. Dynamic mechanical characteristics of NdFeB in electromagnetic brake. Defence Technology 2023; 19: 111–125.

Min

Choi

Kim

, et al. Eddy-current loss analysis of noncontact magnetic device with permanent magnets based on analytical field calculations. IEEE Trans Magn 2015; 51: 1–4.

Park

Choi

Shin

, et al. Torque analysis and measurements of a permanent magnet type eddy current brake with a halbach magnet array based on analytical magnetic field calculations. J Appl Phys 2014; 115: 17E707.

10.

Tian

Ning

, et al. A timesaving transient magneto-thermal coupling model for the eddy current brake. IEEE Trans Veh Technol 2020; 69: 10832–10841.

11.

Yang

Liu

. Eddy current loss and its magnetization effect of electromagnetic buffer under intensive impact load. Int J Appl Electromagnet Mech 2021; 65: 281–300.

12.

, et al. Evaluation and analysis of flux-regulated permanent magnet linear eddy current brakes. IEEE Trans Ind Appl 2023; 59: 712–725.

13.

Yang

. Equivalent subdomain method for performance prediction of permanent magnet eddy current brakes. IET Electr Power Appl 2021; 15: 1174–1186.

14.

Shen

Liu

Sun

, et al. Analysis method of permanent magnet eddy current brake under impact load. Int J Appl Electromagnet Mech 2024; 74: 101–121.

15.

Liu

Peng

, et al. Theoretical and experimental study on magnetic field characteristics of electromagnetic buffer under intensive impact load. IEEE Trans Magn 2025; 61: 1–11.

16.

Yang

Fan

, et al. Irreversible demagnetization mechanism of permanent magnets during electromagnetic buffering. Defence Technology 2021; 17: 763–774.

17.

Yang

Sun

. A hysteresis model considering temperature effects for sintered NdFeB and its application in electromagnetic buffer. J Magn Magn Mater 2021; 537: 168158.

18.

Peng

Chen

Wang

. Analysis and design of adjustable-gap double-side axial flux permanent magnet hysteresis damper. IEEE Trans Ind Electron 2024; 71: 3514–3524.

19.

Zheng

Wang

, et al. Micromagnetic simulation of NdFeB considering spatial–temporal correction. Comput Mater Sci 2024; 239: 112957.

20.

Yang

Fan

. Modeling and optimization of tubular linear permanent magnet eddy current brake with halbach array. Int J Appl Electromagnet Mech 2021; 66: 91–113.

21.

Yang

, et al. Interval uncertain optimization for damping fluctuation of a segmented electromagnetic buffer under intensive impact load. Defence Technology 2021; 17: 884–897.

22.

Yang

Wang

, et al. A robust game optimization for electromagnetic buffer under parameters uncertainty. Eng Comput 2023; 39: 1791–1806.

23.

Liu

Wang

Shen

, et al. Performance study of a permanent-magnet eddy current damper for guns recoil control. J Phys: Conf Ser 2023; 2478: 092021.

24.

Yang

Sun

, et al. Nonlinear analytical model for performance prediction of eddy current recoil brake. IEEE Trans Energy Convers 2022; 37: 1739–1751.

25.

Xie

Liu

Huang

, et al. Comparative study on effect of eddy current damper arrangement on muzzle vibration. Journal of Ordnance Equipment Engineering 2022; 43: 204–209.

26.

Jiles

Atherton

. Theory of ferromagnetic hysteresis. J Magn Magn Mater 1986; 61: 48–60.

27.

Jiles

Thoelke

Devine

. Numerical determination of hysteresis parameters for the modeling of magnetic properties using the theory of ferromagnetic hysteresis. IEEE Trans Magn 1992; 28: 27–35.

28.

Ebrahimi

Khamesee

Golnaraghi

. Permanent magnet configuration in design of an eddy current damper. Microsyst Technol 2010; 16: 19–24.

29.

Wang

Zhang

, et al. Aerodynamic optimal design for a glider with the supersonic airfoil based on the hybrid MIGA-SA method. Aerosp Sci Technol 2019; 92: 224–231.