Damage of permanent magnets in the electromagnetic buffer under intensive impact load

Abstract

This study focused on the damage of permanent magnets(NdFeB) in the electromagnetic buffer under the intensive impact load. According to the constitutive characteristics of the permanent magnets, a constitutive model suitable for NdFeB called damage-modified ZWT was determined. Furthermore, the specific parameters in the constitutive model were determined by fitting the constitutive curve. Then the incremental expression was established for the constitutive model and the model was embedded into ABAQUS via the VUMAT subroutine interface. A dynamic model of the electromagnetic buffer was established and calculated. After that the evolution and distribution law of damage in permanent magnets under the intensive impact load were clarified by analyzing the finite element results. Finally the damage of the device is reduced by optimizing the geometric parameters. The results show that the permanent magnets near the inner ring are damaged more from the radial distribution of the damage and the permanent magnets on both sides are damaged more from the lateral distribution of the damage, which provides a reference for the design of electromagnetic buffers.

Keywords

Intensive impact load the damage in permanent magnets damage-modified ZWT constitutive model VUMAT subroutine

1. Introduction

According to electromagnetic induction theorem, eddy currents are caused due to relative motion of the source magnetic field and the conductor. Eddy currents create magnetic fields and cause a repulsive force proportional to the relative velocity of the field and conductor. The electromagnetic buffer is a new type of brake device designed based on the above principle. Electromagnetic buffers are divided into electromagnetic brakes, permanent magnet brakes and hybrid brakes depending on the excitation source. Compared with electromagnetic brakes and hybrid brakes, permanent magnet brakes do not require external excitation power supply and excitation winding. They save electricity and copper usage, reduce the size of the device, and simplify the structure of the device. On the other hand, their reliability is higher because there is no danger of brake failure during power failure [1].

According to the braking principle of the permanent magnet electromagnetic buffer, it can be seen that permanent magnets are the excitation source of the whole device. Whether permanent magnets can provide sufficient source magnetic field and whether their mechanical properties and magnetic properties are stable are the key factors of the performance of the brake. Due to the high remanence, high coercivity and high magnetic energy product of sintered NdFeB, the use of sintered NdFeB in electromagnetic buffers has great advantages in terms of magnetic field. Many scholars at home and abroad have done a lot of research on the magnetic field distribution and braking force characteristics of electromagnetic buffers during operation [2–5]. However, with the continuous expansion of the application field of electromagnetic buffers, their working conditions are more and more severe, and their load strength is also increasing. Structurally, due to the complex crystal structure and the low slip system, NdFeB have poor plasticity. At the same time, Li et al. [6] pointed out that NdFeB are mostly sintered magnets and there are a certain amount of pores and defects in it, which also reduce the material strength and toughness of NdFeB. Therefore, the permanent magnet electromagnetic buffer can not only consider the magnetic field during the working process under the intensive impact load. The mechanical properties of NdFeB are also critical to the operation of the brake. Wang et al. [7] obtained the engineering stress-strain curves at different strain rates by uniaxial compression and dynamic fracture tests of sintered NdFeB. They considered that The tensile hoop stress, axial stress and shear stress in the expanding specimen cause the surface cracks. Shan Tao et al. [8] studied the effect of different annealing temperatures on the compression curve of the magnet. It is found that the higher the annealing temperature, the lower the compressive stress value. At present, the mechanical properties of NdFeB are mostly studied itself. However there is almost no study of the mechanical properties when it is applied in a specific device.

In this study, the dynamic mechanical properties of NdFeB were described by one-dimensional elastic brittle damage- modified constitutive model and damage-modified Zhu–Wang–Tang (ZWT) constitutive model. By comparison, it was found that the damage-modified ZWT constitutive model can better describe the stress-strain characteristics of NdFeB at different strain rates. The incremental form of damage-modified ZWT constitutive model was established and the material subroutine was written. After that the subroutine was embedded into the dynamic model of electromagnetic buffer to analyze the damage of NdFeB caused by the intensive impact load. Finally the damage of the device is reduced by optimizing the geometric parameters.

2. Damage constitutive model

There are a certain number of pores and fine cracks inside the sintered body. These defects undergo the connection of pores and the growth of micro-cracks to form a dynamic damage process under the intensive impact load [9]. Damage constitutive model is to introduce damage factor D into the existing constitutive model [10]. Material damage behavior is the process of material damage accumulation. When the damage reaches the limit, the material is damaged. The process from no damage to gradual damage can be measured by damage factor D. The Seeger and Jhonson-Cook models, as constitutive models that often describe metal materials, are not accurate enough in describing the mechanical properties of NdFeB. In this paper, a one-dimensional elastic-brittle damage-modified constitutive model and a damage-modified ZWT constitutive model are introduced to fit dynamic mechanical properties of NdFeB and compare.

Through a lot of experimental research on different materials (including metal materials, polymer materials and concrete materials, etc.), it has been consistently shown that: The high-speed deformation process of materials is often accompanied by different forms of internal micro-damage evolution processes, which ultimately lead to the destruction of materials. The observed internal defects / damages of the material are mainly represented by microcracks, microvoids and shear bands on the meso-scale, while dislocations and luan crystals appear on the micro-scale. The various forms of dynamic damage are not actually a simple transient response, but a dynamic process that includes different forms of damage evolution.

Regarding the damage in the form of microcracks, it can be seen from many experimental data that the evolution of damage depends on both strain and strain rate, which is similar to the thermal activation mechanism of dislocation motion. The literature [11] proposes a damage evolution model based on the thermal activation mechanism.

As we all know, the strain rate-dependent constitutive relationship of metal materials is commonly explained by the thermally-activated motion of line defect dislocations inside materials on the microscopic mechanism. $\begin{eqnarray}\displaystyle \dot{{\varepsilon}}=\dot{{\varepsilon}}_{0}\exp \left(-\frac{U_{s}({\sigma})}{kT}\right) & & \displaystyle\end{eqnarray}$ (1) where $\dot{{\varepsilon}}$ is the material strain rate, $\dot{{\varepsilon}}_{0}$ is the frequency factor, U _s(σ) is the thermal activation energy related to the applied stress σ, k is Boltzman constant, and T is the absolute temperature.

Based on the concept of Kachnov’s continuity factor, Rabotonv first defined the damage variable as $\begin{eqnarray}\displaystyle D=(A_{0}-\overline{A})/A_{0}=1-\overline{A}/A_{0} & & \displaystyle\end{eqnarray}$ (2) where D is the damage variable, A ₀ is the initial cross-sectional area of the non-destructive material, and the A is effective load area after the material is damaged. When A = A ₀, D = 0, it means no damage; when A = 0, D = 1, it means that the material is completely damaged.

On this basis, Lemartire proposed in 1971 that the strain response of a damaged element under stress σ is the same as the strain response of a non-destructive element under a defined effective stress σ. The constitutive relationship of the damaged material under external force can be in the form of non-destructive, as long as the Cauchy stress is simply replaced by effective stress. Taking the one-dimensional linear elasticity problem as an example, $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}={\sigma}/\overline{E}=\overline{{\sigma}}/E\end{eqnarray}$ (3) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}=\frac{\overline{{\sigma}}}{E}=\frac{{\sigma}}{E(1-D)}=\frac{{\sigma}}{\overline{E}}\end{eqnarray}$ (4) where σ is the nominal stress on the cross section; σ is the effective stress; E is the elastic modulus of the non-destructive material, and E is the elastic modulus of the damaged material.

The relationship between effective stress σ and nominal stress σ is as follows $\begin{eqnarray}\displaystyle \overline{{\sigma}}={\sigma}/(1-D). & & \displaystyle\end{eqnarray}$ (5)

The constitutive equation of the damaged material in the one-dimensional problem based on the strain equivalence hypothesis is: $\begin{eqnarray}\displaystyle {\sigma}=(1-D)E{\varepsilon}. & & \displaystyle\end{eqnarray}$ (6)

The evolution process of micro-damage is regarded as a stress-promoted thermal activation process, which is similar to the thermal activation motion equation of metal materials. $\begin{eqnarray}\displaystyle {\dot{D}}=\frac{\partial D}{\partial t}={\dot{D}}_{0}\exp \left(-\frac{U_{D}({\sigma})}{kT}\right) & & \displaystyle\end{eqnarray}$ (7) where ${\dot{D}}$ is the evolution rate of microdamage, ${\dot{D}}_{0}$ is the frequency factor of damage evolution, and U _D is the thermal activation energy of microdamage evolution. As an approximation, for the same material, there is a proportional relationship between U _D and U _s, that is $\begin{eqnarray}\displaystyle U_{D}=aU_{S} & & \displaystyle\end{eqnarray}$ (8) where a is the material parameter. So from the formula ((1)) and formula ((7)), we can obtain $\begin{eqnarray}\displaystyle {\dot{D}}/{\dot{D}}_{0}=(\dot{{\varepsilon}}/\dot{{\varepsilon}}_{0})^{a},\quad \text{or }{\dot{D}}=K_{D}\dot{{\varepsilon}}^{a}. & & \displaystyle\end{eqnarray}$ (9) Integrate the above formula to get $\begin{eqnarray}\displaystyle D=K_{D}{\Psi}_{D}(\dot{{\varepsilon}}) & & \displaystyle\end{eqnarray}$ (10) where $K_{D}={\dot{D}}_{0}/\dot{{\varepsilon}}_{0}^{a},{\Psi}_{D}(\dot{{\varepsilon}})=\int _{0}^{t}\dot{{\varepsilon}}^{a}\text{d}t$ .

Assuming that there is a threshold ϵ_th for damage evolution, we can get $\begin{eqnarray}\displaystyle D=K_{D}\dot{{\varepsilon}}^{a-1}({\varepsilon}-{\varepsilon}_{th}),\quad {\varepsilon}>{\varepsilon}_{th}. & & \displaystyle\end{eqnarray}$ (11)

In more general cases, there may be a non-linear relationship with D and ϵ, generalize the above formula to a more general form $\begin{eqnarray}\displaystyle D=K_{D}\dot{{\varepsilon}}^{a-1}({\varepsilon}-{\varepsilon}_{th})^{n},\quad {\varepsilon}>{\varepsilon}_{th}. & & \displaystyle\end{eqnarray}$ (12) In this paper, the damage critical strain ϵ_th is considered to be 0 and the above formula can be shortened as $\begin{eqnarray}\displaystyle D=m{\varepsilon}^{n} & & \displaystyle\end{eqnarray}$ (13) where m is a parameter related to strain rate and n reflects the nonlinear relationship between D and ϵ.

2.1. One-dimensional elastic-brittle damage-modified constitutive model

NdFeB and Al₂O₃ ceramics are both elastic-brittle materials. The internal structure of sintered materials is difficult to produce dislocation and slip. At room temperature, almost no plastic deformation occurs and fracture failure occurs within the elastic range. According to the one-dimensional elastic-brittle damage-modified constitutive model of Al₂O₃ ceramics established by Zhang et al. [12], the one-dimensional elastic-brittle damage-modified constitutive model of NdFeB can be constructed as $\begin{eqnarray}\displaystyle {\sigma}=(1-D)E{\varepsilon}. & & \displaystyle\end{eqnarray}$ (14)

2.2. Damage-modified ZWT constitutive model

The dynamic and static mechanical properties of the sinter are quite different. The damage-modified ZWT constitutive model can describe the mechanical behavior under different strain rates well. It has been widely applied in materials such as concrete, ceramics and plexiglass [13–15].

Fig. 1.

Physical model of ZWT constitutive.

The damage-modified ZWT constitutive model consists of a nonlinear spring and a linear Maxwell element connected in parallel. Its physical model is shown in Fig. 1 and the equation of the constitutive model is expressed as below $\begin{eqnarray}\displaystyle \overline{{\sigma}}=E_{0}{\varepsilon}+{\alpha}{\varepsilon}^{2}+{\beta}{\varepsilon}^{3}+E_{1}\int _{0}^{t}\dot{{\varepsilon}}\exp \left(-\frac{t-{\tau}}{{\theta}_{1}}\right)\text{d}{\tau} & & \displaystyle\end{eqnarray}$ (15) where E ₀, 𝛼 and 𝛽 are elastic constants of the nonlinear spring. E ₁ and θ₁ are elastic coefficient and stress relax time of the Maxwell element respectively, and $\dot{{\varepsilon}}$ is the strain rate. NdFeB is almost coincident and linear in the initial stage under different strain rates. Therefore, the nonlinear spring only considers the primary term E ₀ϵ, ignoring 𝛼ϵ² and 𝛽ϵ³. Then the damage-modified ZWT constitutive model can be written as $\begin{eqnarray}\displaystyle {\sigma}=(1-D)\times \left(E_{0}{\varepsilon}+E_{1}\int _{0}^{t}\dot{{\varepsilon}}\exp \left(-\frac{t-{\tau}}{{\theta}_{1}}\right)\text{d}{\tau}\right). & & \displaystyle\end{eqnarray}$ (16)

2.3. Comparison of fitting curves of two constitutive models

According to the NdFeB dynamic compression experimental curve obtained by Lei et al. [16], the two constitutive models are fitted with the stress-strain curve in the literature by the least squares principle. The fitting comparison are shown in Fig. 2.

When the strain rate are 1456 s⁻¹, 1774 s⁻¹, 2764 s⁻¹, the correlation index R ² of the curve fitting with the one-dimensional elastic-brittle damage-modified constitutive model are 0.832, 0.845 and 0.744, respectively. The correlation index R ² of the curve fitting with the damage-modified ZWT constitutive model are 0.984, 0.986 and 0.983, respectively. On the whole, the damage-modified ZWT constitutive model has a higher degree of fitting with literature values and it can better reflect the stress-strain relationship of NdFeB under different strain rates.

After fitting, the NdFeB damage-modified ZWT constitutive equation is expressed as $\begin{eqnarray}\displaystyle {\sigma}=(1-m{\varepsilon}^{1.00001})\times \left(40129{\varepsilon}+801806\int _{0}^{t}\dot{{\varepsilon}}\exp \left(-\frac{t-{\tau}}{5.7e^{-7}}\right)\text{d}{\tau}\right) & & \displaystyle\end{eqnarray}$ (17) Where $\begin{eqnarray}\displaystyle m=0.0064\dot{{\varepsilon}}+10.49. & & \displaystyle\end{eqnarray}$ (18)

The equation should first be rewritten into an incremental form before embedded into the finite element program. The process of establishing the incremental form will be described in detail below.

Fig. 2.

The fitting and comparison of two constitutive models at different strain rates.

3. Subroutine compilation of Damage-modified ZWT constitutive model

3.1. Incremental form of damage-modified ZWT constitutive model [17]

Considering the finite deformation, the second Kirchhoff stress and Green strain are adopted to expand the ZWT model to the three-dimensional form. $\begin{eqnarray}\displaystyle S_{ij}(E_{ij})=E_{0}AE_{kl}+\int _{0}^{t}E_{1}A\frac{\partial E_{kl}}{\partial {\tau}}\exp \left(-\frac{t-{\tau}}{{\theta}}\right)\text{d}{\tau} & & \displaystyle\end{eqnarray}$ (19) where S _ij is the second Kirchhoff stress tensor; E _ij is the Green strain tensor.

In the case of small deformations, the Green strain is approximately equal to the Cauchy strain. $\begin{eqnarray}\displaystyle E_{ij}={\varepsilon}_{ij}. & & \displaystyle\end{eqnarray}$ (20)

Therefore the above formula (19) can be rewritten as $\begin{eqnarray}\displaystyle S_{ij}({\varepsilon}_{ij})=E_{0}{\varepsilon}_{ij}+\int _{0}^{t}E_{1}\frac{\partial {\varepsilon}_{ij}}{\partial {\tau}}\exp \left(-\frac{t-{\tau}}{{\theta}}\right)\text{d}{\tau}=E_{0}A{\varepsilon}_{kl}+\int _{0}^{t}E_{1}A\frac{\partial {\varepsilon}_{kl}}{\partial {\tau}}\exp \left(-\frac{t-{\tau}}{{\theta}}\right)\text{d}{\tau} & & \displaystyle\end{eqnarray}$ (21) Where $\begin{eqnarray}\displaystyle A=\frac{1}{(1+{\mu})(1+2{\mu})}\times \left[\begin{array}{@{}llllll@{}}1-{\mu} & {\mu} & {\mu} & 0 & 0 & 0\\ {\mu} & 1-{\mu} & {\mu} & 0 & 0 & 0\\ {\mu} & {\mu} & 1-{\mu} & 0 & 0 & 0\\ 0 & 0 & 0 & {\displaystyle \frac{(1-2{\mu})}{2}} & 0 & 0\\ 0 & 0 & 0 & 0 & {\displaystyle \frac{(1-2{\mu})}{2}} & 0\\ 0 & 0 & 0 & 0 & 0 & {\displaystyle \frac{(1-2{\mu})}{2}}\end{array}\right] & & \displaystyle\end{eqnarray}$ (22) where μ is Poisson’s ratio.

The process of establishing the incremental form of the second term in equation (21) is as follows. $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle S_{ij,(t)}^{v}=E_{1}\int _{0}^{t}A\frac{\partial {\varepsilon}_{kl}}{\partial {\tau}}\exp \left(-\frac{t-{\tau}}{{\theta}}\right)\text{d}{\tau}\end{eqnarray}$ (23) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\rm\Delta}S_{ij,(t+{\rm\Delta}t)}^{v}=S_{ij,(t+{\rm\Delta}t)}^{v}-S_{ij,(t)}^{v}\end{eqnarray}$ (24) $\begin{eqnarray}\displaystyle {\rm\Delta}S_{ij,(t+{\rm\Delta}t)}^{v} & = & \displaystyle \left[1-\exp \left(\frac{-{\rm\Delta}t}{{\theta}}\right)\right]\times \frac{E_{1}{\theta}}{{\rm\Delta}t}\times A{\rm\Delta}{\varepsilon}_{kl}^{(t+{\rm\Delta}t)}\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle -\left[1-\exp \left(\frac{-{\rm\Delta}t}{{\theta}}\right)\right]\times \left[E_{1}\int _{0}^{t}A\frac{\partial {\varepsilon}_{kl}^{(t+{\rm\Delta}t)}}{\partial {\tau}}\exp \left(-\frac{t-{\tau}}{{\theta}}\right)\text{d}{\tau}\right].\end{eqnarray}$ (25)

Converting the integral term in the equation (25), we obtain equation (26). $\begin{eqnarray}\displaystyle E_{1}\int _{0}^{t}A\frac{\partial {\varepsilon}_{kl}^{(t+{\rm\Delta}t)}}{\partial {\tau}}\exp \left(-\frac{t-{\tau}}{{\theta}}\right)\text{d}{\tau}=E_{1}{\theta}\left[1-\exp \left(\frac{-t}{{\theta}}\right)\right]\times A\frac{{\rm\Delta}{\varepsilon}_{kl}^{(t+{\rm\Delta}t)}}{{\rm\Delta}t}. & & \displaystyle\end{eqnarray}$ (26)

In this way, the incremental form of the ZWT constitutive model is obtained. $\begin{eqnarray}\displaystyle {\rm\Delta}S_{ij}^{(t+{\rm\Delta}t)}=E_{0}A{\rm\Delta}{\varepsilon}_{kl}^{(t+{\rm\Delta}t)}+E_{1}{\theta}\times A\frac{{\rm\Delta}{\varepsilon}_{kl}^{(t+{\rm\Delta}t)}}{{\rm\Delta}t}\times \left[1-\exp \left(\frac{-{\rm\Delta}t}{{\theta}}\right)\right]\times \exp \left(\frac{-t}{{\theta}}\right). & & \displaystyle\end{eqnarray}$ (27)

Kirchhoff stress is also called pseudo-stress and Cauchy stress is real stress. The relationship between Kirchhoff stress and Cauchy stress is $\begin{eqnarray}\displaystyle {\sigma}_{ij}^{(t+{\rm\Delta}t)}=({\rho}/{\rho}_{0})F_{ik}^{(t+{\rm\Delta}t)}S_{kl}^{(t+{\rm\Delta}t)}F_{lj}^{(t+{\rm\Delta}t)^{\text{T}}} & & \displaystyle\end{eqnarray}$ (28) where σ_ij is the Cauchy stress tensor; F is the deformation gradient tensor; the density ratio 𝜌∕𝜌₀ is equal to 1 when the material is incompressible.

The expression of the damage factor is $\begin{eqnarray}\displaystyle D=m{\varepsilon}^{n}. & & \displaystyle\end{eqnarray}$ (29)

Since the damage is unrecoverable, it’s necessary to determine whether the strain is increasing. $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}{\rm\Delta}D^{(t+{\rm\Delta}t)}=\left\{\begin{array}{@{}ll@{}}m[|{\varepsilon}^{(t+{\rm\Delta}t)}|^{n}-|{\varepsilon}^{(t)}|^{n}] & {\rm\Delta}{\varepsilon}^{(t+{\rm\Delta}t)}\geq 0\\ 0 & {\rm\Delta}{\varepsilon}^{(t+{\rm\Delta}t)}<0,\end{array}\right.\\[15.0pt] {\rm\Delta}{\varepsilon}^{(t+{\rm\Delta}t)}=|{\varepsilon}^{(t+{\rm\Delta}t)}|-|{\varepsilon}^{(t)}|.\end{array} & & \displaystyle\end{eqnarray}$ (30)

The incremental expression for the damage factor is $\begin{eqnarray}\displaystyle D^{(t+{\rm\Delta}t)}=D^{(t)}+{\rm\Delta}D^{(t+{\rm\Delta}t)}. & & \displaystyle\end{eqnarray}$ (31)

The expression of the stress tensor containing damage is $\begin{eqnarray}\displaystyle S_{ij,d}^{(t)}=(1-D^{(t)})S_{ij}^{(t)}. & & \displaystyle\end{eqnarray}$ (32)

The incremental form of the stress tensor containing damage is $\begin{eqnarray}\displaystyle {\rm\Delta}S_{ij,d}^{(t+{\rm\Delta}t)} & = & \displaystyle S_{ij,d}^{(t+{\rm\Delta}t)}-S_{ij,d}^{(t)}\nonumber\\ \displaystyle & = & \displaystyle (1-D^{(t+{\rm\Delta}t)})S_{ij}^{(t+{\rm\Delta}t)}-(1-D^{(t)})S_{ij}^{(t)}\nonumber\\ \displaystyle & = & \displaystyle (1-D^{(t+{\rm\Delta}t)})(S_{ij}^{(t)}+{\rm\Delta}S_{ij}^{(t+{\rm\Delta}t)})-(1-D^{(t)})S_{ij}^{(t)}\nonumber\\ \displaystyle & = & \displaystyle (1-D^{(t+{\rm\Delta}t)}){\rm\Delta}S_{ij}^{(t+{\rm\Delta}t)}-(D^{(t+{\rm\Delta}t)}-D^{(t)})S_{ij}^{(t)}\nonumber\\ \displaystyle & = & \displaystyle (1-D^{(t+{\rm\Delta}t)}){\rm\Delta}S_{ij}^{(t+{\rm\Delta}t)}-{\rm\Delta}D^{(t+{\rm\Delta}t)}S_{ij}^{(t)}\end{eqnarray}$ (33)

The Kirchhoff stress is converted to Cauchy stress according to equation (16): $\begin{eqnarray}\displaystyle {\sigma}_{ij,d}^{(t+{\rm\Delta}t)} & = & \displaystyle F_{ik}^{(t+{\rm\Delta}t)}S_{kl,d}^{(t+{\rm\Delta}t)}F_{lj}^{(t+{\rm\Delta}t)^{\text{T}}}\nonumber\\ \displaystyle & = & \displaystyle F_{ik}^{(t+{\rm\Delta}t)}(S_{kl,d}^{(t)}+{\rm\Delta}S_{kl,d}^{(t+{\rm\Delta}t)})F_{lj}^{(t+{\rm\Delta}t)^{\text{T}}}\nonumber\\ \displaystyle & = & \displaystyle F_{ik}^{(t+{\rm\Delta}t)}S_{kl,d}^{(t)}F_{lj}^{(t+{\rm\Delta}t)^{\text{T}}}+(1-D^{(t+{\rm\Delta}t)})\times (F_{ik}^{(t+{\rm\Delta}t)}{\rm\Delta}S_{kl,d}^{(t)}F_{lj}^{(t+{\rm\Delta}t)^{\text{T}}})\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle -∼{\rm\Delta}D^{(t+{\rm\Delta}t)}\times F_{ik}^{(t+{\rm\Delta}t)}S_{kl}^{(t)}F_{lj}^{(t+{\rm\Delta}t)^{\text{T}}}.\end{eqnarray}$ (34)

Fig. 3.

ABAQUS calculation process of embedded material subroutine.

3.2. Writing and verification of ABAQUS VUMAT subroutine

Although ABAQUS has a powerful library of material constitutive models, it can not cover all situations. When the user needs to define the constitutive model of the material autonomously, the material subroutine interface VUMAT of ABAQUS is needed. The running process of the subroutine is shown in Fig. 3. At the beginning of the incremental step, the main program passes the initial values to the corresponding variables via the subroutine interface. The subroutine VUMAT processes and updates these variables and then returns the processed variables to the main program at the end of the subroutine [18,19]. For the damage-modified ZWT constitutive model mentioned in this paper, the calculation flow of ABAQUS calling subroutine is as follows, where the specific process of updating the value of each variable is given by the above formula.

The damage-modified ZWT constitutive equation was implemented in Fortran language to set up its incremental form. As shown in Fig. 4, the numerical simulation model of a single NdFeB was established and the subroutine was embedded into the ABAQUS solver for finite element analysis. It can be seen from Fig. 5 that the fitted values are in good agreement with the simulated values at each strain rate. This verifies the correctness of the written subroutine. Therefore, the following part will establish the dynamic model of the electromagnetic buffer and associate it with the VUMAT subroutine of NdFeB. And the dynamic response of electromagnetic buffer under the intensive impact load and the damage distribution of permanent magnets will be analyzed.

Fig. 4.

Uniaxial compression model of a single permanent magnet.

Fig. 5.

Comparison of fitted values of damage-modified ZWT constitutive model and simulation values under different strain rates.

4. Dynamic model of the electromagnetic buffer

As stated in the introduction, electromagnetic buffers utilize the relative motion between the conductor and the source magnetic field. The two relative moving parts of brakes are usually referred to as the primary and secondary. As shown in Fig. 6, the electromagnetic buffer includes primary: magnetic shoes, permanent magnets, the moving rod, the guide ring and the end nut; secondary: the conductor cylinder and the end cover. The material properties of each part are shown in Table 1.

Fig. 6.

Structure diagram of the electromagnetic buffer. 1 - End cover; 2 - The moving rod; 3 - Conductor cylinder; 4 - Guide ring; 5 - Magnetic shoes; 6 - Permanent magnets; 7 - End nut.

Table 1

Material properties sheet for each part

Part name	Material	Density/(g/cm³)	Elastic modulus/MPa	Poisson’s ratio
Permanent magnets	N38	7.5	Defined by material subroutine	0.24
Magnetic shoes	DT4c	7.8	80650	0.29
Guide ring	Bronze	8.9	110000	0.34
End nut	45 steel	7.8	210000	0.30
The moving rod	Non-magnetic steel	7.8	206000	0.25

The dynamic model of the electromagnetic buffer is shown in Fig. 7. When the electromagnetic buffer is in operation, the intensive impact load is applied to the left end of the brake, so that the primary accelerate to the left. When the magnetic shoes and the permanent magnets move relative to the conductor cylinder, it generates the damping force opposite to the movement direction. At the same time, the displacement-related resistance is applied to the right side of the moving rod to simulate the mechanical buffering force during the actual working process. The curves of the intensive impact load and the mechanical buffering force are shown in Fig. 8 and Fig. 9, respectively. The contact collision relations are defined where the moving rod, the guide ring, magnetic shoes, permanent magnets and the end nut are in contact.

Fig. 7.

Dynamic model of the electromagnetic buffer.

Fig. 8.

Curve of the intensive impact load.

Fig. 9.

Curve of the mechanical buffering force.

5. Result analysis

The stress distribution of the electromagnetic buffer is determined by the external load. In the initial stage, the intensive impact load rises rapidly in a short time and reaches its peak near 5 ms. During this time, the intensive impact load is transmitted to the entire brake through the moving rod, and the overall stress of the brake is also rising. At 5.40 ms, the peak stress appears at the collision between the right end nut and the right guide ring, and the peak stress is 590.4 MPa. After that, the intensive impact load decreases rapidly, which results in the stress decrease of the electromagnetic buffer.

Fig. 10.

Stress distribution in the cross-section of the electromagnetic buffer at different times.

Fig. 11.

Damage factor distribution of the rightmost permanent magnet at different times.

The rightmost permanent magnet is taken to analyze the damage over time. The calculation results in Fig. 11 show that the damage factor D is increasing throughout the whole operation of the electromagnetic buffer. At the same time, it can be clearly seen that the value of the damage factor gradually decreases from the inner ring to the outer ring. This is due to the fact that the intensive impact load first acts on the moving rod. During the process of transferring the load to the permanent magnets, the moving rod and the permanent magnets generate a large friction force, resulting in a large strain of the permanent magnet inner ring. As the radial distance increases, the other positions are not in direct contact with the moving rod. So the strain at other positions is small, which results in the gradual decrease of the damage factor from the inside to the outside.

Fig. 12.

The position of three points (A,B,C) on the permanent magnet.

Fig. 13.

The curve of damage factor with time.

In order to better analyze the distribution law of the damage factor along the radial direction of the permanent magnets, three points ABC are taken out as shown in Fig.12. Through the analysis according to Fig. 13, it can be found that the two points of BC do not directly contact with the moving rod and the damage factor decreases with the increase of the distance from the inner ring. The damage factors of BC at the end of the operation are 0.119 and 0.061, respectively. However, due to the constant contact and collision between the inner ring and the moving rod during the operation, the damage factor at point A is relatively large and the damage factor value at the end of the operation is 0.189.

As stated above, we define the damage value according to the load capacity under the same strain. As the formula (35) shown below, there is a difference between the stress for the same strain when the material is damaged and not damaged. This is because the damage of internal microstructure leads to the degradation of the macroscopic elastic modulus E of the material. The damage is characterized by the degradation of elastic modulus and the relation between the equivalent elastic modulus after damage and the elastic modulus of undamaged material is E = (1 − D)E. When the damage value is 0.189, it means that the degradation degree of elastic modulus of the material under the damage is 0.189 and equivalent elastic modulus E = (1 − 0. 189)E. $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}{\varepsilon}={\sigma}/\overline{E}=\overline{{\sigma}}/E\\[7.0pt] \displaystyle {\varepsilon}=\frac{\overline{{\sigma}}}{E}=\frac{{\sigma}}{E(1-D)}=\frac{{\sigma}}{\overline{E}}.\end{array} & & \displaystyle\end{eqnarray}$ (35)

As can be seen in Fig. 14, with the increase of the intensive impact load, the extreme value of strain almost appeared at the same time as the stress peak. Later, due to the combined effect of the intensive impact load, eddy current buffering force and friction, the strain value fluctuates continuously and gradually decreases to 0 in the fluctuation with the advance of time. According to the definition of damage in the ZWT constitutive model, damage occurs when the absolute value of strain increases in two adjacent analysis steps as shown in the formula (36). And because the damage value D is irreversible, although the overall strain is decreasing, the damage factor D is still increasing due to the existence of fluctuations. $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}{\rm\Delta}D^{(t+{\rm\Delta}t)}=\left\{\begin{array}{@{}ll@{}}m[|{\varepsilon}^{(t+{\rm\Delta}t)}|^{n}-|{\varepsilon}^{(t)}|^{n}] & {\rm\Delta}{\varepsilon}^{(t+{\rm\Delta}t)}\geq 0\\ 0 & {\rm\Delta}{\varepsilon}^{(t+{\rm\Delta}t)}\lt 0,\end{array}\right.\\[10.0pt] {\rm\Delta}{\varepsilon}^{(t+{\rm\Delta}t)}=|{\varepsilon}^{(t+{\rm\Delta}t)}|-|{\varepsilon}^{(t)}|.\end{array} & & \displaystyle\end{eqnarray}$ (36)

Figure 15 shows the cross-section damage factor distribution of permanent magnets at different times. When the intensive impact load increases rapidly, the moving rod generates a large acceleration along the load direction and the right side of the permanent magnets have a large contact collision with the end nut. Therefore, the damage factor of the right-end permanent magnets increases rapidly within 5.4 ms and the maximum value reaches 0.018. In the next time, the permanent magnets are subjected to eddy current damping force, frictional contact force with the moving rod and collision force with the guide ring. The results show that the combined force of the left-end magnetic steel is larger and the damage is greater.

Fig. 14.

Axial strain history curve and damage value at point A.

Fig. 15.

Damage factor distribution in the cross-section of permanent magnets at different times.

Fig. 16.

Distribution and local enlargement of damage factor below 0.2.

In order to show the damage degree of permanent magnets more clearly, according to the definition above, the upper limit of the damage factor is set to 0.2 and the damage distribution map at 135 ms is analyzed. As shown in Fig. 16, it can be seen that the leftmost permanent magnet has a uniform circumferential destruction. At the same time, according to Fig. 16 and the radial damage distribution of the rightmost permanent magnet in Fig. 12, it can be known that the peak damage factor value of the rightmost permanent magnet is 0.189, and no destruction occurs. Therefore, the damage of the left permanent magnet is larger than that of the right permanent magnet at 135 ms. Additionally, as can be seen from the cloud diagram in Fig. 16, the damage cloud colors of the permanent magnets on the left and right sides are basically green, yellow, and red, but the middle part is blue and cyan. According to the numerical comparison on the right, it can be seen that the damage value of the permanent magnets on both sides is greater than the damage value of the permanent magnets in the middle part.

6. The risk analysis and parameter optimization of the permanents

In fact, for different materials and different constitutive model, methods which used in the risk assessment are not the same. The theory and simulation should be compared with the experiment so that a more accurate evaluation method can be obtained. As mentioned in the introduction of this article, most of the previous articles for sintering of NdFeB focused on its magnetic properties. We lack evaluation standards similar to other metal and non-metallic materials. Since the permanent magnets provide the source magnetic field in the buffer, the permanent magnets are equivalent to the power system of the whole device. The damage of the permanent magnet’s mechanical properties will also have a great impact on its magnetic properties. In this paper, for the sake of conservative safety assessment, 0.2 is taken as the critical damage value of permanent magnet and the damage of degree (DOD) is defined according to the critical damage value: $\begin{eqnarray}\displaystyle DOD=\frac{\text{Number of grids that exceed the critical damage value}}{\text{Number of grids for all permanent magnets}}. & & \displaystyle\end{eqnarray}$ (37)

According to this definition method, it can be calculated that the DOD of the permanent magnets in this paper is 0.208%. It can be seen that the DOD is small. However, the permanent magnet on the left side has a large damage, which will also have an impact on the whole. The DOD of the left-most permanent magnet is 5.0%. According to this definition of DOD, the geometric dimension is adjusted to optimize the damage degree.

Fig. 17.

Schematic diagram of parameter position.

As presented in Fig. 17, there are three parameters of permanent magnets, namely its inner diameter r, outer diameter R and height H. At the same time, the magnetic shoe acts as a magnetically permeable material in the magnetic field and as a contact piece between permanent magnets in the structural field. Its height h will also affect the magnitude of the eddy current buffering force and the damage degree of the permanent magnets. Therefore, the four parameters of h, r, R, and H are used as design variables and the objective function is to minimize the damage. At the same time, it is considered that both the buffer can weaken the impact load and the structure of the buffer should be as compact as possible. After preliminary trial calculation, the value range and initial value of the four parameters are determined in Table 2.

Table 2

The initial value and value range of four parameters

Parameters	h/mm	H/mm	r/mm	R/mm
Initial value	10	20	20	66.5
Lower limit	8	18	17.5	65
Upper limit	12	22	22.5	68

Latin hypercube design is an experimental design method based on space filling. It can make the limited samples fill the whole design space as much as possible and weaken the requirement of sample boundary. Latin hypercube design has the characteristics of super nonlinear response fitting ability, which can better reflect the physical essence of the research object. The optimal Latin hypercube design greatly improves the uniformity of the Latin hypercube design by adding a criterion and it is especially suitable for multi-factor, multi-layer experiments and the situations where the system model is completely unknown.

In view of the strong nonlinear characteristics in the working process of the buffer, this paper selects the optimal Latin hypercube experimental design to obtain sample points. A total of 20 training samples are generated and the samples are shown in Table 3.

In this paper, the method of constructing the proxy model is used to simplify the complicated calculation process. We optimize the proxy model to improve computing efficiency. Common methods for constructing proxy models include response surface method, radial basis function method, and neural network method. The neural network [20] has strong non-linear mapping capabilities. A 3-layer BP neural network can fit any complex non-linear function. So this paper uses neural network to build a proxy model. BP neural network modeling process is shown in Fig. 18.

Fig. 18.

BP neural network modeling process.

Genetic algorithm (GA) [21] is a method to search the optimal solution by simulating the biological evolution process in nature. According to the determined fitness function, we use selection, crossover, and mutation as the three main genetic operators to operate on individuals in the population. Populations evolve through continuous exchange of chromosomal information between individuals. Eventually, individuals with good fitness values are retained and individuals with poor fitness values are continuously eliminated during the evolution process. After constant evolution, individuals in the population gradually approach the optimal solution. Moreover, since genetic algorithm does not care about the mathematical relationship between design variables and objective functions, nor does it need information such as the gradient of objective functions to determine the search direction, it is very suitable for parameter optimization.

Table 3

Tables of traning samples

Group	h/mm	H/mm	r/mm	R/mm	Group	h/mm	H/mm	r/mm	R/mm
1	8.00	20.95	18.03	66.74	11	10.11	21.79	22.50	66.90
2	8.21	21.58	19.08	66.58	12	10.32	19.68	20.13	66.42
3	8.42	22.00	20.92	66.27	13	10.53	18.00	18.29	67.69
4	8.63	19.05	19.87	65.32	14	10.74	19.47	19.61	65.47
5	8.84	20.53	17.77	67.37	15	10.95	19.26	22.24	66.11
6	9.05	20.11	17.50	67.84	16	11.16	20.74	20.66	65.94
7	9.26	18.84	21.71	65.79	17	11.37	19.89	18.56	68.00
8	9.47	20.32	21.19	67.53	18	11.58	18.42	18.82	65.63
9	9.68	21.37	20.40	65.16	19	11.79	18.21	19.34	67.06
10	9.89	18.63	21.98	65.00	20	12.00	21.16	21.45	67.21

Table 4

Parameters’ initial value and optimized value

Parameters	h/mm	H/mm	r/mm	R/mm
Initial value	10.00	20.00	20.00	66.50
Optimized value	12.00	22.00	18.37	65.33

In this paper, the neural network model is optimized by genetic algorithm, and the optimized parameters are shown in Table 4.

Fig. 19.

The optimized result cloud map.

According to the optimization results in Fig. 19, it can be seen that the peak value of the damage factor is 0.819, which is 5.5% lower than the initial 0.867. At the same time, according to the statistical results of the grid, the optimized DOD is 0.105%, which is 49.5% lower than the previous. The peak value decreases less, but the overall grid damage degree decreases more obviously. The safety of the device has been greatly improved. It also illustrates the effectiveness of the optimization method.

7. Conclusion

In this paper, the strain-rate dependent dynamic constitutive equation of N38 is established. The incremental form is established and embedded into the dynamic simulation model through the VUMAT subroutine. The damage evolution and distribution of the permanent magnets in the specific device are obtained and the specific conclusions are as follows

(1) Compared with the one-dimensional elastic-brittle damage-modified constitutive model, the damage-modified ZWT constitutive model has a higher degree of fitting with the stress-strain curve of sintered NdFeB. Therefore, the damage-modified ZWT constitutive model can better describe the dynamic constitutive characteristics of sintered NdFeB at different strain rates.

(2) It can be seen from the uniaxial compression simulation of a single permanent magnet model that the simulated values are in good agreement with the fitted values. This shows that the established incremental form is relatively correct.

(3) From the radial damage distribution, the damage of NdFeB decreases gradually from the inner ring to the outer ring. From the lateral damage distribution, the NdFeB on both sides are damaged more and the NdFeB in the middle part are less damaged. Therefore, the contact parts at both ends between the permanent magnets and the moving rod are damaged greatly, which should be paid attention to in design and use.

Footnotes

Acknowledgements

The work was primarily supported by the National Natural Science Foundation of China (grant number 301070603).

References

Kou

B.Q.

Jin

Y.X.

Zhang

and Zhang

H.L.

, Development and application prospects of the electromagnetic damper, Proceedings of the CSEE 35(12) (2015), 3132–3143.

Ebrahimi

Khamesee

M.B.

and Golnaraghi

, A novel eddy current damper: Theory and experiment, Journal of Physics D-Applied Physics 42(7) (2009), 6.

Canova

and Vusini

, Analytical modeling of rotating eddy-current couplers, IEEE Transactions on Magnetics 41(1) (2005), 24–35.

Chen

Z.Q.

Tian

J.Y.

Huang

Z.W.

and Wang

J.H.

, Calculation and test correction method of plane type eddy current damping coefficient, China Journal of Highway and Transport 29(10) (2016), 46–53.

L.Z.

Liu

Y.P.

and Li

D.S.

, Optimization design and test of dual air-gaps and liquid-cooled eddy current retarder, International Journal of Applied Electromagnetics and Mechanics 61(2) (2019), 157–170.

A.H.

Dong

S.Z.

and Li

, Mechanical properties of RE permanent magnetic materials, Metallic Functional Materials (04) (2002), 7–10.

Wang

H.R.

Wan

Chen

D.N.

Lei

G.H.

and Ren

C.Y.

, Dynamic fracture of sintered Nd-Fe-B magnet under uniaxial compression, Journal of Magnetism and Magnetic Materials 456 (2018), 358–367.

Tao

Ahmad

Zhang

P.Y.

Zheng

X.M.

and Zhang

S.Y.

, Structural, magnetic and mechanical properties of Nd6.8Fe63.3Nb4.0Co1.3Cr1.1Zr0.9B22.6 nanocrystalline hard magnet, Aip Advances 8(7) (2018), 8.

Rojek

Nosewicz

and Chmielewski

, Micro-macro relationships from discrete element simulations of sintering, International Journal for Multiscale Computational Engineering 15(4) (2017), 323–342.

10.

Song

Chen

J.K.

Qian

Y.F.

and Cao

Y.H.

, Rate-dependent characteristic of relaxation time of concrete, Acta Mechanica Solida Sinica 32(1) (2019), 69–80.

11.

L.L.Wang

X.L.Dong

and Sun

Z.J.

, Dynamic constitutive behavior of materials at high strain rate taking account of damage evolution, Explosion and Shock Waves (03) (2006), 193–198.

12.

Zhang

X.Q.

Y.X.

Ning

J.G.

Zhao

L.M.

and Yang

G.T.

, A study on the strain-rate dependent dynamic constitutive equation of Al2O3 ceramics, Explosion and Shock Waves (03) (2004), 226–232.

13.

Ahmed

Wang

and Zhang

, High-speed bird impact analysis of aircraft windshield by using a nonlinear viscoelastic model, Applied Mechanics and Materials 290 (2013), 85–90.

14.

Chang

B.N.

Wang

X.Y.

Long

Z.Q.

J.F.

Ruan

S.L.

and Shen

C.Y.

, Constitutive modeling for the accurate characterization of the tension behavior of PEEK under small strain, Polymer Testing 69 (2018), 514–521.

15.

Wang

Y.J.

and Zhang

W.H.

, Finite element simulation of PMMA aircraft windshield against bird strike by using a rate and temperature dependent nonlinear viscoelastic constitutive model, Composite Structures 108 (2014), 21–30.

16.

Lei

G.H.

Wan

and Wang

H.R.

, Study on fracture behavior of sintered Nd-Fe-B magnet under dynamic compression, Chinese Rare Earths 39(06) (2018), 63–68.

17.

Lili

Xinlong

Dong

and Zijian

, Dynamic constitutive behavior of materials at high strain rate taking account of damage evolution, Explos Shock Waves (03) (2006), 193–198.

18.

Liu

W.S.

Z.P.

and Qing

G.H.

, A progressive damage model introducing temperature field for bolted composite joint with preload, Modelling and Simulation in Materials Science and Engineering 27(6) (2019), 19.

19.

Zheng

and Qingying

, Development of a practical uniaxial cyclic constitutive model for concrete, IOP Conference Series: Earth and Environmental Science 295 (2019), 042057 (10 pp.).

20.

Alexandrino

P.D.L.

Gomes

G.F.

and Cunha

S.S.

, A robust optimization for damage detection using multiobjective genetic algorithm, neural network and fuzzy decision making, Inverse Problems in Science and Engineering 28(1) (2020), 21–46.

21.

Y.H.

Chen

Bai

J.T.

and Zuo

W.J.

, Shape optimization of thin-walled cross section for automobile body considering stamping cost, manufacturability and structural stiffness, International Journal of Automotive Technology 21(2) (2020), 503–512.