Probabilistic Life Calculation Method of NdFeB Based on Brittle Fatigue
Damage Model

This paper proposes a probabilistic life calculation method of NdFeB based on brittle fatigue damage model. Firstly, Zhu-Wang-Tang (ZWT) constitutive model considering strain rate is established, and based on this, a numerical co-simulation model for NdFeB life calculation is constructed. The life distribution diagram of NdFeB under different stress levels is obtained after simulation. Secondly, a new model of brittle fatigue damage based on brittle damage mechanism is proposed. Then the parameters in the model are identified according to the life distribution diagram of NdFeB and the parameter distribution of the damage evolution model when applied to NdFeB is obtained. Finally, the probability density evolution equation of NdFeB life calculation is established and solved using the probability density evolution method. Probability density function (PDF) of NdFeB life under different stress levels is obtained and provides theoretical basis for the reliability of NdFeB in engineering applications.


Introduction
Neodymium-iron-boron (NdFeB) is a third-generation rare earth permanent magnet material that was developed in the 1980s and successfully used in production [1]. Sintered neodymium-iron-boron is made by using powder metallurgy and goes through the process of making the calcined material into fine powder, pressing it into a blank, and then sintering. Because of its excellent magnetic properties, it is widely used in computer, communication, medical, machinery, aerospace, national defense and other fields. At present, most of the research on sintered NdFeB focuses on improving its magnetic property, so that it can obtain higher coercive force and magnetic energy product [2,3], but few studies focus on mechanical and fatigue damage properties. However, as NdFeB is increasingly used in motors, eddy current dampers and other fields, it is often subjected to cyclic loads, and research on its fatigue performance cannot be ignored.
Damage caused by repeated action of a fixed or variable amplitude load is called fatigue damage. Fatigue damage is one of the most common failure forms of engineering structures. The evolution of fatigue damage is an irreversible change in the microstructure of materials due to the action of external forces. The accurate prediction of fatigue life depends on a reasonable criterion of fatigue cumulative damage [4]. With the gradual deepening of fatigue research, many scholars at home and abroad have proposed a variety of life prediction theories and analysis approaches based on a large number of experiments and theoretical studies for different research fields, industries and objects [5][6][7][8]. Typical approaches include stress-based approach, strain-based approach, energy-based approach, and continuum damage mechanics approaches. The above life prediction theories and approaches can be basically divided into two categories. One is that fatigue damage evolution function is constructed directly by fitting the fatigue damage curve based on the test data, which cannot reflect the microscopic damage accumulation process of materials. The other is to construct a fatigue damage model based on the continuous damage theory and the dissipation potential based on the irreversible thermodynamics framework, and then to fit and verify the fatigue damage model based on the experimental data. The damage in this model is related to the plastic strain. However, due to the difficulty in producing slippage and dislocation inside the sinter, NdFeB is a typical elastic-brittle material [9] and hardly has plastic deformation during loading. Therefore, neither of the above two types of models can well describe the fatigue damage process of NdFeB. In this paper, a new brittle fatigue damage model is established to describe the fatigue characteristics of NdFeB based on the mechanism of irreversible thermodynamics and damage micromechanics [10].
The fatigue studies mentioned above belong to the category of deterministic analysis and the dispersion of fatigue life is seldom considered. In real life, some mechanical parameters of the same batch of materials are also different due to manufacturing technology and other reasons. At the same time, it is difficult to ensure the absolute consistency of the loading conditions during the fatigue life test, and it is also difficult to ensure the complete uniformity of the use occasions in specific use cases. A variety of reasons ultimately lead to the dispersion of fatigue life in fatigue test or specific use. Although deterministic analysis can describe the basic law of fatigue damage failure, it cannot reflect the probabilistic characteristics of fatigue life. Therefore, the introduction of probability analysis methods in fatigue life analysis and calculation has become inevitable [11][12][13]. Li et al. [14] established a probability density evolution equation for the nonlinear dynamic response of random structures based on the principle of preservation of probability. The difference method is then used to solve the equation and the probability density function (PDF) of the response is obtained. The proposed probability density evolution method for nonlinear dynamic response analysis of random structures is considered to be an effective method for reliability analysis, error identification, random vibration and seismic safety assessment [15][16][17]. Zhang et al. [18] applied the probability density evolution method to the fatigue life analysis of concrete and obtained the probability density function of fatigue life with different loading levels, which provided a new method for the uncertainty analysis of fatigue life.
As a further exploration of this approach, this paper endeavors to obtain a probabilistic life calculation method based on brittle fatigue damage model. The main research work of this paper is summarized as follows: firstly, on the acquisition of NdFeB life data, this paper adopts the co-simulation method of ABAQUS and nCode. The dynamic constitutive model of the NdFeB is determined and embed in ABAQUS in the form of a VUMAT subroutine. Then in nCode, the calculation process is built according to the fatigue five block diagram. After calculation, NdFeB fatigue life distribution map can be obtained which takes the uncertainty of strain rate, the stress amplitude of load and the compressive strength limit of the material into account. Subsequently, based on the irreversible thermodynamic theoretical framework and damage micromechanics, a new brittle fatigue damage model based on the mechanism of brittle damage is established to describe the fatigue characteristics of NdFeB. This is followed by identifying the distribution of parameters in the model based on the obtained fatigue life distribution map. Finally, the probability density evolution method is applied to fulfill the probabilistic assessment of NdFeB fatigue life, through which a universe probability density evolution equation for NdFeB fatigue life is developed. By solving this equation with TVD finite difference method, the PDF of NdFeB fatigue life due to different fatigue loading levels can be obtained. This paper mainly does threefold contributions. Initially NdFeB is a strain rate-related material, and accurate establishment of its dynamic constitutive model is the key to obtain life data. The ZWT constitutive model adopted in this paper can well describe the dynamic characteristics of NdFeB. After that, the present fatigue damage model cannot take into account the microscopic damage process of brittle sinter. The fatigue damage model based on the brittle damage mechanism used in this paper is expressing the evolution of damage by the growth of microcracks. It is a model derived from Meso-scale Representative Volume Element (RVE) and established to represent the macro-scale brittle fatigue damage model. Finally, the probability density evolution method is applied to the calculation of the life of NdFeB. The probability density evolution equation for calculating the life of NdFeB is established and the probability density distribution of NdFeB with the stress level is obtained.

Strain Rate Dependent Brittle Constitutive Model of NdFeB
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 under pressure. This will cause a reduction in material carrying capacity. Modified constitutive model is to introduce bearing capacity reduction factor B into the existing constitutive model. B is a strain-related coefficient and can explain the nonlinear relationship of stress and strain of elastic materials. 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 modified constitutive model and a modified ZWT constitutive model are introduced to fit dynamic mechanical properties of NdFeB and compare.

Two Constitutive Models of NdFeB
(1) One-dimensional elastic-brittle modified constitutive model where r is the stress and e is the strain of the material; E is the elastic modulus of the material; B is the bearing capacity reduction factor.
where m is the parameter related to the strain rate; n is the strain index.
(2) Modified ZWT constitutive model The dynamic and static mechanical properties of the sinter are quite different. The 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 [19][20][21].
The modified ZWT constitutive model consists of a nonlinear spring and a linear Maxwell element connected in parallel (Fig. 1), and the equation of the constitutive model is expressed as below.
where E 0 , a and β are elastic constants of the nonlinear spring. E 1 and θ 1 are elastic coefficient and stress relax time of the Maxwell element respectively, and _ e 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 0 ε, ignoring αε 2 and βε 3 . Then the modified ZWT constitutive model can be written as

Comparison of Two Constitutive Models
According to the NdFeB dynamic compression experimental curve obtained by Lei et al. [22], 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 is 1456, 1774, 2764 s −1 , the correlation index R 2 of the curve fitting with the onedimensional elastic-brittle modified constitutive model is 0.832, 0.845 and 0.744, respectively. The correlation index R 2 of the curve fitting with the modified ZWT constitutive model is 0.984, 0.986 and 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.

Writing and Verification of ABAQUS VUMAT Subroutine
Considering the finite deformation, the second Kirchhoff stress and Green strain are adopted to expand the ZWT model to the three-dimensional form [23].
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.
Therefore the above formula (7) can be rewritten as where μ is Poisson's ratio.
The process of establishing the incremental form of the second term in Eq. (9) is as follows.
In this way, the incremental form of the ZWT constitutive model is obtained.
Kirchhoff stress is also called pseudo-stress and Cauchy stress is real stress. The relationship between Kirchhoff stress and Cauchy stress is where σ ij is the Cauchy stress tensor; F is the deformation gradient tensor; the density ratio ρ/ρ 0 is equal to 1 when the material is incompressible.
The expression of the bearing capacity reduction factor is Since the factor is unrecoverable, it is necessary to determine whether the strain is increasing.
The incremental expression for the bearing capacity reduction factor is The expression of the stress tensor containing the bearing capacity reduction factor is S The incremental form of the stress tensor containing the bearing capacity reduction factor is The Kirchhoff stress is converted to Cauchy stress according to Eq. (16): Although ABAQUS has a powerful library of material constitutive models, it cannot 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. 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 [24,25]. For the modified ZWT constitutive model mentioned in this paper, the calculation flow of ABAQUS calling subroutine is shown in Fig. 3, where the specific process of updating the value of each variable is given by the above formula.
The 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.

Calculation Method of the Fatigue Life in nCode
It can be seen from this classic five block diagrams of the fatigue analysis in Fig. 6 that there are three necessary pieces of information to obtain before a structural fatigue analysis. They are the cyclic load that the  The fitting value at 1456s -1 The simulation value at 1456s -1 The fitting value at 1774s -1 The simulation value at 1774s -1 The fitting value at 2764s -1 The simulation value at 2764s -1   7 shows the co-simulation of ABAQUS and nCode. The first is to input the results of the finite element calculation into nCode. The results include not only the stress and strain information, but also the load curve. The load curve enters nCode step by step as the initial condition for fatigue calculation. There is no need to define an additional load spectrum. According to the five block diagrams of fatigue analysis, the Co-simulation process built in nCode is shown in the figure. Firstly the result file is imported from ABAQUS into nCode. Secondly the fatigue properties of the material are defined. Finally, the results are presented as cloud maps and numerical values after calculations.

Co-Simulation Calculation and Results
The dispersion of fatigue life is an inevitable problem in the actual fatigue test even if the same batch of test pieces is used and under the same loading conditions. The reasons for the dispersion of fatigue life are quite complex. In this paper, the effects of loading rate, stress amplitude and random distribution of compressive strength limit on life are mainly considered. According to the constitutive model in the first section, NdFeB is a strain rate-dependent material. The dispersion of the loading rate and the dispersion of the stress amplitude will cause the dispersion of the finite element results. The compressive strength limit is an inherent property of the material and an important parameter for the fatigue performance of the material. It is difficult to unify product attributes due to the difficulty in completely unifying specific production processes. Comprehensive consideration of the distribution of loading rate, stress amplitude, and compressive strength limit will reflect the causes of fatigue life dispersion to a certain extent and make the life calculation result more realistic.
Random sampling was performed using the optimal Latin hypercube sampling method for strain rate, stress amplitude, and compressive strength limit. The parameter ranges are shown in Tab. 1.   Fig. 10, microcrack growth is used to represent the evolution of damage in RVE and the classic crack growth law (Paris formula) is assumed to be suitable for crack growth in mesoscopic RVE, then a brittle fatigue damage model that can represent the macro scale can be established.
The microscopic RVE damage D is defined as the average value of all microelement damages. At the same time, it is considered that the specimen breaks when the area of all microcracks is equal to the RVE surface area, that is The damage of RVE is simplified as Assuming that the crack surfaces of all microelements are equal, then Using the macro Paris formula: where da=dN means the crack growth rate, a is the crack length, N is the cycles, c and g are material parameters, F is the geometric shape factor, r and Dr are nominal stress amplitude and stress amplitude at crack, K max and K min are the maximum and minimum values of the intensity factor at the crack, respectively. According to the strain equivalence principle, the nominal stress in the formula is replaced with the effective stress which is " r ¼ r=ð1 À DÞ. At the same time, the Paris formula can be regarded as the integral of the rate equation of a cycle, that is The crack width of the micro-element is defined as e and the derivative is obtained L L Micro element Microcrack d Figure 10: The Meso-scale representative volume element Combine formulation (26)-(28): According to formula (28), a ¼ L 2 D ne can be obtained. Substituting A ¼ cp g=2 F g ðneÞ g=2À1 L gÀ2 into Eq. (29) and it can be rewritten as When integrating the above formula, substitute D g=2 approximately as ðN =N f Þ g=2 into it and N f can be formulated as follows:

Parameter Distribution Identification of Brittle Fatigue Damage Model
It can be seen from the comparison between the model fitting value and the simulation mean that fatigue damage evolution model based on brittle damage mechanism has good fitting accuracy for the fatigue life of permanent magnets (Fig. 11). At the same time, the change in the value of g and A causes the dispersion of the fatigue life curve. The parameter distribution of g and A will be used to describe the dispersion of external conditions and the resulting dispersion of fatigue life.
According to the fatigue damage evolution model of brittle materials described in Section 3.1, it can be seen in Fig. 12 that under the same stress level, the fatigue life will change according to the two parameters g and A. In order to identify these two parameters, the identification criteria are defined as follows

Probability Density Evolution Equation for NdFeB Fatigue Life and Its Solution Procedure
Based on the above generalized probability density evolution, the probability density evolution equation for fatigue life is established. The stress amplitude Dr is taken as the independent variable and the fatigue life of NdFeB N f as the dependent variable. The fatigue life distribution will be given based on the stress amplitude Dr. According to the introduction above, the stochastic fatigue life N f can be written as: In which Θ ¼ Θðg; AÞ.
Substituting Formula (42)  (1) Select points in the probability space and determine their assigned probability.
A series of representative discrete points are selected from the distribution space Â of random vector Θ and denoted as Θ q ¼ ðh q;1 ; h q;2 ; …; h q;s Þ; q ¼ 1; 2; …; n s where n s is the total number of the selected points. The assigned probability of each point in the point set is calculated through the following equation where V q represents the volume of the selected point in the distribution space. For a given point Θ ¼ h q , solve deterministic Eq. (31) and we can get the generalized velocity, _ Z j ðh q ; tÞ, at the same time.
(3) Solve the generalized probability density evolution equation.
With the representative point set and the assigned probability from (1), the Eq. (40) can be rewritten as the following discrete form @p ZΘ ðz; h q ; tÞ @t þ X m l¼1 h l ðh q ; tÞ @p ZΘ ðz; h q ; tÞ @z l ¼ 0 Note the initial conditions as Introducing the generalized velocity obtained from (2) into Eq. (46) and the TVD finite difference method is used to solve the equation.
For each discrete representative point, the corresponding @p ZΘ ðz; h q ; tÞ can be achieved by the above solution step and P Z ðz; tÞ can be achieved by integrating q from 1 to n s .

Calculation of NdFeB Fatigue Life Distribution
The analyzed results are shown in Figs. 15-17. Fig. 15 is a 3d cloud image of stress loading level, life and corresponding probability density. Fig. 16 is a top view of Fig. 15. Fig. 16 shows the value of probability density in color. Fig. 17 shows the probability density distribution of different life spans at stress levels of 0.75, 0.78, 0.80 and 0.85.  As can be seen from Figs. 15 and 16, the overall life decreases with the continuous increase of the stress loading level. At the same time, it can be seen that the initial rate of decline is faster, and the decrease range tends to be flat after the stress level is further increased. It is further found from Fig. 17 that when the stress level is high, the probability density distribution of life is more compact, while the life is more dispersed when the stress level is low. The lifetime distribution of NdFeB obtained by solving the probability density evolution equation can provide a theoretical basis for the reliability analysis of the lifetime of permanent magnets in engineering.
From the comparison between the simulation value and the results predicted by the probability density evolution method in Fig. 18, the mean value and standard variance value are all close. The prediction method Figure 19: Full-text computing flow chart in this paper can reflect the simulation data to a certain extent. And from the simulation data, the life under other stress amplitudes can also be effectively predicted, which provides a theoretical basis for the reliability of NdFeB in engineering applications.

Conclusion
This article aimed to explore the probabilistic life calculation method of NdFeB based on brittle fatigue damage model. The computing flow chart of the whole paper is shown in Fig. 19. The following conclusions can be drawn: 1. ZWT constitutive model considering strain rate is established, and based on this, a numerical cosimulation model for NdFeB life calculation is constructed. The life distribution diagram of NdFeB under different stress levels is obtained after simulation.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.