A Modified Fractional Maxwell Numerical Model for Constitutive Equation of Mn-Cu Damping Alloy

To better describe its constitutive relation, we need a new constitutive equation for an important nonlinear elastic material, Mn-Cu damping alloy. In this work, we studied the nonlinear and hysteretic characteristics of the stress-strain curve of the M2052 alloy with the uniaxial cyclic tensile test with constant strain rate. The strain rate and amplitude correlations of M2052 resembled those of nonlinear viscoelastic material. Therefore, we created a new constitutive equation for the M2052 damping alloy by modifying the fractional Maxwell model, and we used the genetic algorithm to carry out numerical fitting with MATLAB. By comparing with the experimental data, we confirmed that the new constitutive equation could accurately depict the nonlinear constitutive relation and hysteretic property of the damping alloy. Taken together, this new constitutive equation for Mn-Cu damping alloy based on the fractional Maxwell model can serve as an effective tool for further studies of the constitutive relation of the Mn-Cu damping alloys.


Introduction
Mn-Cu damping alloy is a typical twin-damping alloy. In the 1940s, Zener [1] first developed Mn-20Cu alloy. The Sonoston alloy (Mn-36.2Cu-3.49Al-3.04Fe-1.17Ni wt. %) was developed in the UK for submarine propellers, which effectively reduced the noise of the propeller [2]. International Copper Research Association (INCRA) developed the Incramute alloy (Mn-48.1Cu-1.55Al-0.27Si wt. %), which is mainly used in machines and pedestals [3]. Russia also developed the ABPOPA alloy (Mn-Cu-Al-Fe-Ni-Zn) in 1975 and achieved a good performance in the field of precision machinery manufacturing and ship equipment [4]. However, the machinability and applicable temperature range of the above Mn-Cu alloys were relatively low, until the M2052 (Mn-20Cu-5Ni-2Fe) alloy was developed by Kawahara et al. [5], who added Ni, Fe, Co, Ti and other elements to Mn-20Cu alloy. M2052 alloy has strong stiffness and machinability. Especially, it has high damping at room temperature, which broadens its application range.
At present, M2052 alloy has been widely used in many industries. Xin [6] preliminarily explored the application of M2052 damping alloy on cradle of remote control weapon station. Wang [7] applied M2052 damping alloy to the motor support and fan housing. Baochang Liu et al. [8] added the powder of Mn-Cu damping alloy into the diamond drill, and found that the 40% Mn-Cu damping alloy drill had the same performance as the original one, but the vibration performance was reduce by 4.3%, and the drilling process was more stable. Komatsu et al. [9] described how apply the M2052 alloy as structural material in spacecraft vibration reduction. Yan, Shan et al. [10] mainly analyzed various characteristics of Mn-Cu damping alloy, and designed a turbine containing this alloy, which was proved to be effective in vibration reduction through experiments. On M2052 steel plate after heat treatment, the required 150 mm × 20 mm × 2 mm tensile samples were cut by wire electrical discharge machining, and the detailed size are shown in Figure 1. At last, we simulated the nonlinear constitutive relationship of damping alloy and compared with experimental curves.

Uniaxial Cyclic Tensile Test with Constant Strain Rate
Firstly, in order to study the nonlinear constitutive relationship of Mn-Cu damping alloy, the uniaxial cyclic tensile test at constant strain rate was be used on the damping alloy by MTS universal testing machine. Then the hysteresis stress-strain curve of the material at different strain rates and strain amplitude values can be obtained, and will provide data for the model simulation and verification.

Test Materials and Equipment
The composition of M2052 alloy used in the test is Mn-22.1Cu-5.24Ni-1.93Fe (mass%), and the main performance indexes are shown in the following Table 1. On M2052 steel plate after heat treatment, the required 150 mm × 20 mm × 2 mm tensile samples were cut by wire electrical discharge machining, and the detailed size are shown in Figure 1. The test equipment choose the MTS-Landmark 810 universal testing machine produced by USA MTS Company (Minneapolis, MN, USA), as shown in Figure 2. The bearing capacity range is 0-100 kN, the test sample frequency is 10 Hz, the length of the extensometer is 50 mm, and the installation details are shown in Figure 3. The test equipment choose the MTS-Landmark 810 universal testing machine produced by USA MTS Company (Minneapolis, MN, USA), as shown in Figure 2. The bearing capacity range is 0-100 kN, the test sample frequency is 10 Hz, the length of the extensometer is 50 mm, and the installation details are shown in Figure 3.

Experiment Scheme and Result Analysis
Due to the high strength and damping capacity of M2052 alloy, it is often used in industry as anti-vibration structural parts. It should not produce plastic deformation, so the cyclic strain amplitude of test should not exceed 0.2%. Three types of uniaxial cyclic tensile tests were carried out. Type 1: A constant strain rate of 0.0025%/s was performed on the first tensile sample, with three different cycle strain amplitudes of 0.05%, 0.1% and 0.15%. Type 2: A constant strain rate of 0.005%/s was performed on the second tensile sample, with three different cycle strain amplitudes of 0.05%, 0.1% and 0.15%. Type 3: A constant strain rate of 0.01%/s was performed on the third tensile sample, with three different cycle strain amplitudes of 0.05%, 0.1% and 0.15%. All tests' initial preload was 0.8 kN, and the loading and unloading process was controlled by computer program, at last outputting the stress-strain curve.

Experiment Scheme and Result Analysis
Due to the high strength and damping capacity of M2052 alloy, it is often used in industry as anti-vibration structural parts. It should not produce plastic deformation, so the cyclic strain amplitude of test should not exceed 0.2%. Three types of uniaxial cyclic tensile tests were carried out. Type 1: A constant strain rate of 0.0025%/s was performed on the first tensile sample, with three different cycle strain amplitudes of 0.05%, 0.1% and 0.15%. Type 2: A constant strain rate of 0.005%/s

Experiment Scheme and Result Analysis
Due to the high strength and damping capacity of M2052 alloy, it is often used in industry as anti-vibration structural parts. It should not produce plastic deformation, so the cyclic strain amplitude of test should not exceed 0.2%. Three types of uniaxial cyclic tensile tests were carried out. Type 1: A constant strain rate of 0.0025%/s was performed on the first tensile sample, with three different cycle strain amplitudes of 0.05%, 0.1% and 0.15%. Type 2: A constant strain rate of 0.005%/s was performed on the second tensile sample, with three different cycle strain amplitudes of 0.05%, 0.1% and 0.15%. Type 3: A constant strain rate of 0.01%/s was performed on the third tensile sample, with three different cycle strain amplitudes of 0.05%, 0.1% and 0.15%. All tests' initial preload was 0.8 kN, and the loading and unloading process was controlled by computer program, at last outputting the stress-strain curve. The tests were carried out only under strain control. In each type test, the strain rate was fixed. Firstly, stretching the sample to 0.05% strain amplitude, and then unloading it. After returning to the initial value, stretching it to 0.1% strain amplitude, and then unloading it. Finally, stretching the sample to 0.15% strain amplitude and then unloading it. The schematic diagram of test loading-unloading process is in Figure 4.
The tests were carried out only under strain control. In each type test, the strain rate was fixed. Firstly, stretching the sample to 0.05% strain amplitude, and then unloading it. After returning to the initial value, stretching it to 0.1% strain amplitude, and then unloading it. Finally, stretching the sample to 0.15% strain amplitude and then unloading it. The schematic diagram of test loadingunloading process is in Figure 4.  The raw data contained some noise, which led hardly to show the characteristics of hysteresis curves with different conditions. Therefore, this work used the wavelet filtering method to filter and smooth the experimental data of uniaxial cyclic tensile tests. The result is shown in  In addition, the results of the original test data were put into the supplementary material (Figures S1-S3 and Tables S1 and S2).  The raw data contained some noise, which led hardly to show the characteristics of hysteresis curves with different conditions. Therefore, this work used the wavelet filtering method to filter and smooth the experimental data of uniaxial cyclic tensile tests. The result is shown in  In addition, the results of the original test data were put into the supplementary material (Figures  S1-S3 and Tables S1 and S2). The tests were carried out only under strain control. In each type test, the strain rate was fixed. Firstly, stretching the sample to 0.05% strain amplitude, and then unloading it. After returning to the initial value, stretching it to 0.1% strain amplitude, and then unloading it. Finally, stretching the sample to 0.15% strain amplitude and then unloading it. The schematic diagram of test loadingunloading process is in Figure 4.  The raw data contained some noise, which led hardly to show the characteristics of hysteresis curves with different conditions. Therefore, this work used the wavelet filtering method to filter and smooth the experimental data of uniaxial cyclic tensile tests. The result is shown in Figures 5-7. In addition, the results of the original test data were put into the supplementary material (Figures S1-S3 and Tables S1 and S2). . Hysteresis curves with different strain amplitudes at strain rate of 0.0025%/s (  is stress,  is strain, a  is strain amplitude).

Figure 5.
Hysteresis curves with different strain amplitudes at strain rate of 0.0025%/s (σ is stress, ε is strain, ε a is strain amplitude).  . Hysteresis curves with different strain amplitudes at strain rate of 0.005%/s (σ is stress, ε is strain, ε a is strain amplitude). ε(%) Figure 6. Hysteresis curves with different strain amplitudes at strain rate of 0.005%/s ( σ is stress, ε is strain, a ε is strain amplitude). Figures 5-7 are hysteresis curves with the same strain rate but different strain amplitude. It can be seen from Tables 2 and 3 that at the same strain rate, the hysteretic area increased with strain amplitude, and the slope of the curve decreased with strain amplitude. The hysteretic area represents the damping performance of the damping alloy. This indicates that the damping capacity is positively correlated with the amplitude within the elastic range, which has been proved by the literature [26]. However, it is found from Figures 5-7 that the strain amplitude of hysteresis curve was generally not equal to the set strain value, which is due to the error caused by the response gain of the hydraulic servo valve, but it would not affect trend analysis. The strain amplitude error between the set value and the measured value of each cycle is listed as shown in Table 4. . Hysteresis curves with different strain amplitudes at strain rate of 0.01%/s (σ is stress, ε is strain, ε a is strain amplitude).
Figures 5-7 are hysteresis curves with the same strain rate but different strain amplitude. It can be seen from Tables 2 and 3 that at the same strain rate, the hysteretic area increased with strain amplitude, and the slope of the curve decreased with strain amplitude. The hysteretic area represents the damping performance of the damping alloy. This indicates that the damping capacity is positively correlated with the amplitude within the elastic range, which has been proved by the literature [26]. However, it is found from Figures 5-7 that the strain amplitude of hysteresis curve was generally not equal to the set strain value, which is due to the error caused by the response gain of the hydraulic servo valve, but it would not affect trend analysis. The strain amplitude error between the set value and the measured value of each cycle is listed as shown in Table 4. From Table 4, we can see the measured value of strain amplitude was not equal to the set value. The error value between the measured and the set value was not the same, but all little. The relative error ranged between 0.039% and 0.175%. The measured strain amplitude at the strain rate of 0.0025%/s and 0.005%/s were less than the set value, and it was greater than the set value at the strain rate of 0.01%/s. Because the error caused by response gain of hydraulic servo valve was small, the measured strain amplitude at the strain rate of 0.0025%/s and 0.005%/s were also close. However, for the measured strain amplitude at the strain rate of 0.01%/s, the error caused by response gain of the hydraulic servo valve was large. Although the measured value of strain amplitude is inconsistent with the set value, the measured strain amplitude at each cycle with the same strain rate was still monotonically increasing. It can be concluded from Tables 2 and 3 that with the increase of strain amplitude, the hysteresis area increased correspondingly and the slope decreased correspondingly. Therefore, the inconsistency between measured and set value will not influence the trend analysis of hysteresis curve versus strain amplitude.
Since the strain amplitudes of all tests are different, it is usually impossible to analyse the effect of the strain rate change on the hysteretic curve. However, the measured strain amplitudes at the strain rates of 0.0025%/s and 0.005%/s are close to each other, so their measured strain amplitudes can be approximately equal. It can be seen from the data in Tables 2 and 3 that the strain rate has a certain influence on the hysteretic curve, but the specific law still needs further experimental study.
Because the error between the measured strain amplitude and the set value can be neglected, the original set value of strain amplitude is still used to describe the problem easily. Moreover, the use of measured data does not affect the fitting of the numerical model, but can avoid the loss of useful data points due to operations such as filtering and smoothing. Therefore, this paper only uses the measured data after smoothing to analyze the results of test, and the measured data without filtering is still used in the subsequent research.

Establishment of Fractional Maxwell Model
The basic component of the fractional-order viscoelastic model is the spring-pot element, which is a fractional-order model between the spring that represent pure elasticity and the Newton's viscoelastic model that represent pure viscosity. Its constitutive relation can be written as [20]: where, κ is the quasi-property [27], whose unit is Pa · s α , and its expression is Equation (2) [28].
This τ = η E is called relaxation time, η is the viscosity and E is the spring stiffness [29]. Quasi-property is the numerical measurement of a dynamic process, that it is not a simple material property. It relates with the elastic modulus, stress, time, viscosity and other physical quantities, and also associated with mathematical quantities such as fractional index (α). As shown in Figure 8, if α = 0, the spring-pot element can be simplified as a spring element (σ = Eε, κ= E). If α = 1, the spring-pot element simplified to a dashpot element (σ = η . ε, κ = η). It can be seen that the quasi-property makes the model parameters more closely to the actual research object.
One spring-pot mechanism element (σ 1 , ε 1 ) in series with one spring element (σ 2 , ε 2 ) form the fractional Maxwell model of M2052 damping alloy (as shown in Figure 9). κ represents the quasi-property of spring-pot. α is the fractional order and 0 ≤ α ≤ 1. E represent the Young's modulus of M2052 alloy. σ 1 is the stress of spring-pot. ε 1 is the strain of spring-pot. σ 2 is the stress of spring. ε 2 is the strain of spring. σ is the stress of M2052 alloy. ε is the strain of M2052 alloy. It can be obtained by Maxwell model relation: is called relaxation time,  is the viscosity and E is the spring stiffness [29].
Quasi-property is the numerical measurement of a dynamic process, that it is not a simple material property. It relates with the elastic modulus, stress, time, viscosity and other physical quantities, and also associated with mathematical quantities such as fractional index (  ). As shown in Figure 8 Figure 8. Schematic diagram of spring-pot element.
One spring-pot mechanism element 11 ( , )  in series with one spring element 22 ( , )  form the fractional Maxwell model of M2052 damping alloy (as shown in Figure 9).  represents the quasi-property of spring-pot.  is the fractional order and 01   . E represent the Young's modulus of M2052 alloy. 1  is the stress of spring-pot. 1  is the strain of spring-pot. Substitute Equation (1) into Equation (3) and Equation (4) to get: The constitutive equation of fractional Maxwell model is:

The Reduces Vibration Mechanism of Mn-Cu Damping Alloy
In order to establish the governing equation of M2052 damping alloy according to the law of energy conservation, the mechanism of energy consumption of M2052 damping alloy must be understood. The high damping characteristics of Mn-Cu alloy are derived from the fact that in the process of melting cooling at high temperature, paramagnetism turns to anti-ferromagnetism after Neel point, and anti-ferromagnetic martensite twins are generated. The relaxation process of the twin boundary leads to the consumption of external energy. The parent phase of Mn-Cu damping alloy is austenite, which contains martensite twins at room temperature, so it has high damping characteristics at room temperature.
The microstructure of twin type damping alloy under external force as shown in Figure 10, from 20-micron scale under the electron microscope figure can see a horse twin throughout the austenite The constitutive equation of fractional Maxwell model is:

The Reduces Vibration Mechanism of Mn-Cu Damping Alloy
In order to establish the governing equation of M2052 damping alloy according to the law of energy conservation, the mechanism of energy consumption of M2052 damping alloy must be understood. The high damping characteristics of Mn-Cu alloy are derived from the fact that in the process of melting cooling at high temperature, paramagnetism turns to anti-ferromagnetism after Neel point, and anti-ferromagnetic martensite twins are generated. The relaxation process of the twin boundary leads to the consumption of external energy. The parent phase of Mn-Cu damping alloy is Materials 2020, 13, 2020 9 of 36 austenite, which contains martensite twins at room temperature, so it has high damping characteristics at room temperature.
The microstructure of twin type damping alloy under external force as shown in Figure 10, from 20-micron scale under the electron microscope figure can see a horse twin throughout the austenite grain size. The 1-micron scale of damping alloy under the electron microscope can see its internal grain subdivide to many fine twins, the twins stayed parallel and cross state, when under stress.

The Reduces Vibration Mechanism of Mn-Cu Damping Alloy
In order to establish the governing equation of M2052 damping alloy according to the law of energy conservation, the mechanism of energy consumption of M2052 damping alloy must be understood. The high damping characteristics of Mn-Cu alloy are derived from the fact that in the process of melting cooling at high temperature, paramagnetism turns to anti-ferromagnetism after Neel point, and anti-ferromagnetic martensite twins are generated. The relaxation process of the twin boundary leads to the consumption of external energy. The parent phase of Mn-Cu damping alloy is austenite, which contains martensite twins at room temperature, so it has high damping characteristics at room temperature.
The microstructure of twin type damping alloy under external force as shown in Figure 10, from 20-micron scale under the electron microscope figure can see a horse twin throughout the austenite grain size. The 1-micron scale of damping alloy under the electron microscope can see its internal grain subdivide to many fine twins, the twins stayed parallel and cross state, when under stress.  Figure 11 is a schematic diagram of the movement of twins under stress [30]. Figure 11a is the lattice structure of the alloy without stress, and the points in the figure represent the metal atoms. When the shear stress was applied, the lattice structure began to deform, as shown in Figure 11b. When the stress reached a certain degree, twins began to generate, as shown in Figure 11c. As the stress increased, the twin band also began to increase and became a large twin crystal, as shown in Figure 11d. The deformed part and the undeformed part of the lattice formed a symmetrical relationship on the twin boundary, and the lattice direction of the twin band all changed, as shown in Figure 11f. Moreover, just because of this lattice redirection way, the energy could be absorbed. As the stress continued to increase, the large twin crystal began to differentiate into several small twin crystal, as shown in Figure 11e.  Figure 11 is a schematic diagram of the movement of twins under stress [30]. Figure 11a is the lattice structure of the alloy without stress, and the points in the figure represent the metal atoms. When the shear stress was applied, the lattice structure began to deform, as shown in Figure 11b. When the stress reached a certain degree, twins began to generate, as shown in Figure 11c. As the stress increased, the twin band also began to increase and became a large twin crystal, as shown in Figure 11d. The deformed part and the undeformed part of the lattice formed a symmetrical relationship on the twin boundary, and the lattice direction of the twin band all changed, as shown in Figure 11f. Moreover, just because of this lattice redirection way, the energy could be absorbed. As the stress continued to increase, the large twin crystal began to differentiate into several small twin crystal, as shown in Figure 11e. M2052 alloy contains martensite twins and austenite at room temperature, and the austenite is the parent phase. When the alloy is subjected to stress at room temperature, the stress-induced martensitic transformation will occur. This new martensite is called the stress-induced martensite. Since the martensitic transformation, there will be friction between the martensitic phase and austenitic phase, also between the original martensitic phase and new martensitic phase, which is M2052 alloy contains martensite twins and austenite at room temperature, and the austenite is the parent phase. When the alloy is subjected to stress at room temperature, the stress-induced martensitic transformation will occur. This new martensite is called the stress-induced martensite. Since the martensitic transformation, there will be friction between the martensitic phase and austenitic phase, also between the original martensitic phase and new martensitic phase, which is also a reason for the high damping of M2052 alloy. The martensitic phase transformation includes two types: stress induced martensitic transformation and strain induced martensitic transformation. The relationship between them is shown in Figure 12.   , for the stress exceeding the yield strength of austenite, the plastic deformation will occur. The martensitic transformation generated by the deformation is called strain-induced martensite (BF). In addition, as the temperature increases, some martensite will begin to transform into austenite (BD). When the temperature is greater than d M , the martensite will no longer be produced, and all phase will become austenite. Since M2052 alloy works at room temperature and plastic deformation is not allowed, stress-induced martensitic transformation will occur under the stress.

Governing Equation
Therefore, from the perspective of conservation of energy, under the condition of constant temperature and thermal insulation, the work done by external force all into strain energy and it should be equal to the elastic potential energy and internal energy. Moreover, the internal energy is due to the twin's relaxation movement and the interface slip between stress-induced martensite and austenite phase [31]. At present, it is generally considered that the friction energy between martensite and austenite, martensite and martensite is small and can be ignored for the time being, so it is assumed that the strain energy is the sum of elastic potential energy and twin energy. The expression is shown in Equation (7).
where W is the strain energy per unit volume, W is the elastic potential energy per unit volume, and W is the twin energy per unit volume. (σ is stress. σ B is the critical stress of stress-induced martensite. σ C is the nominal stress of stress-induced martensite at T 2 . σ F is the stress of strain-induced martensite at T 2 . σ D is the stress required for the transformation of austenite. T is temperature. Ms is the start temperature of the martensitic transformation. Md is the end temperature of the martensitic transformation. T 1 is the room temperature.).
As shown in Figure 12. When M S < T < M S1 and 0 < σ < σ B , the martensitic phase transformation is mainly induced by stress (AB), and the stress and temperature show a linear relationship in the elastic range. When M S1 < T < M d and σ B < σ, for the stress exceeding the yield strength of austenite, the plastic deformation will occur. The martensitic transformation generated by the deformation is called strain-induced martensite (BF). In addition, as the temperature increases, some martensite will begin to transform into austenite (BD). When the temperature is greater than M d , the martensite will no longer be produced, and all phase will become austenite. Since M2052 alloy works at room temperature and plastic deformation is not allowed, stress-induced martensitic transformation will occur under the stress.

Governing Equation
Therefore, from the perspective of conservation of energy, under the condition of constant temperature and thermal insulation, the work done by external force all into strain energy and it should be equal to the elastic potential energy and internal energy. Moreover, the internal energy is due to the twin's relaxation movement and the interface slip between stress-induced martensite and austenite phase [31]. At present, it is generally considered that the friction energy between martensite and austenite, martensite and martensite is small and can be ignored for the time being, so it is assumed that the strain energy is the sum of elastic potential energy and twin energy. The expression is shown in Equation (7).
where W S is the strain energy per unit volume, W E is the elastic potential energy per unit volume, and W T is the twin energy per unit volume. The schematic diagram of energy conversion was shown in Figure 13. σ loading is the stress of loading stage. σ unloading is the stress of unloading stage. σ e is elastic stress. σ t is the stress of twin. Each energy density is the area, which is surround by each curve and x-axis. The area between σ loading and σ unloading is the loss of energy per unit volume. The strain energy per unit volume is equal to the σ e surrounding area and σ t surrounding area combined. Each energy density is the area, which is surround by each curve and x-axis. The area between loading  and unloading  is the loss of energy per unit volume. The strain energy per unit volume is equal to the e  surrounding area and t  surrounding area combined. According to the J2 deformation theory, Hooke's law and Equations (1)-(4), we can obtain blow equations [32] under the unit volume: Then taking Equation (8) into Equation (7) we can receive Equation (9).
The formula can convert into Equation (10): The above equation is the governing equation of fractional Maxwell model of M2052 damping alloy.
In the governing equation, when 0 t  , then ( ) 0, ( ) 0 tt monotonically increasing in the loading segment and monotonically decreasing in the unloading segment. Therefore, the initial value and boundary conditions of the governing equation are shown in Equation (11): According to the J2 deformation theory, Hooke's law and Equations (1)-(4), we can obtain blow equations [32] under the unit volume:
The formula can convert into Equation (10): The above equation is the governing equation of fractional Maxwell model of M2052 damping alloy. In the governing equation, when t ≤ 0, then σ(t) = 0, ε(t) = 0 ; when t > 0, σ(t) and ε(t) is monotonically increasing in the loading segment and monotonically decreasing in the unloading segment. Therefore, the initial value and boundary conditions of the governing equation are shown in Equation (11):

Numerical Solution
Taking the derivatives of both sides of the governing Equation (10) with respect to time, we can get the below equations.
In order to solve the above equation, the fractional derivative is discretized by the finite difference method.
Define t m = m∆t, m = 0, 1, 2, . . . , K, 0 ≤ k ≤ K; ∆t = T K stands the time step. According to the definition of fractional calculus of Caputo [33], when 0 < h < 1, then a = 0, n = 1, D h t σ(t) (0 < h < 1) can be written as: where, according to the finite difference method, where In addition, because of the idea of difference, t k+1 = (k + 1)∆t, by the same token t k , t i+1 , t i , the above equation can be changed into: Taking Equations (15) and (17) into Equation (14), and defining j = k − i we can obtain: The coefficient in the higher-order term o((∆t) r ) is the higher-order error term.
So we can get the expression of D h t σ(t k ), as shown in Equation (19).
Take Equations (15) and (21) into Equation (13), and the higher-order error terms can be omitted:

Model Parameter Analysis
Taking the fractional Maxwell numerical model and boundary conditions into MATLAB software, and then combining with the test data determined the time step ∆t = 0.1 s. We selected the more regular curve which was under strain rate 0.0025%/s and strain amplitude 0.1% as the basic experiment data and σ(0)= 4.5901 MPa, ε(0)= 0.00903%. Then we analyzed the influence of α and κ values on the model. α is the fractional order coefficient, and its range is 0 < α < 1, and the range of quasi-properties κ can be determined by E and relaxation time according to Equation (2). Although no data on the relaxation time of M2052 alloy have been found (which can be measured by future relaxation tests), we can know when α = 0, the minimum value of κ is 68.5 GPa.
When C = 100 and α = 0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, the stress-time fitting curve of the loading stage with strain rate of 0.0025%/s and strain amplitude of 0.1% is shown in Figure 14. As can be seen from Figure 14, fractional Maxwell can simulate the nonlinear curve of the convex function in the loading section. Moreover, with the increase of the  , the curvature is higher and higher, but the maximum stress is lower. In order to observe the influence of the  , choose  that at the largest curvature as an invariant. When =0.9  and =100,500, 2500,125000 C , the stress-time fitting curve of the loading stage with the strain rate of 0.0025%/s and the strain amplitude of 0.1% is shown in Figure 15. It can be seen from Figure 15 that the increase of  will reduce the curvature of the curve, but can significantly increase the increment of the stress over time. Therefore, there is an optimal solution between  and  .

Genetic Algorithm Setting
Since the governing Equation (21) is a difference iterative form with time as the iterative variable, the least square method or linear regression method cannot be used, so the fitting problem can be transformed into a multi-objective optimization problem. The objective function, design variables and constraint conditions are prepared by using the GA toolbox of MATLAB 2016R. The optimization process is shown in Figure 16. As can be seen from Figure 14, fractional Maxwell can simulate the nonlinear curve of the convex function in the loading section. Moreover, with the increase of the α, the curvature is higher and higher, but the maximum stress is lower. In order to observe the influence of the κ, choose α that at the largest curvature as an invariant.
When α = 0.9 and C = 100, 500, 2500, 12500, the stress-time fitting curve of the loading stage with the strain rate of 0.0025%/s and the strain amplitude of 0.1% is shown in Figure 15. As can be seen from Figure 14, fractional Maxwell can simulate the nonlinear curve of the convex function in the loading section. Moreover, with the increase of the  , the curvature is higher and higher, but the maximum stress is lower. In order to observe the influence of the  , choose  that at the largest curvature as an invariant. When =0.9  and =100,500, 2500,125000 C , the stress-time fitting curve of the loading stage with the strain rate of 0.0025%/s and the strain amplitude of 0.1% is shown in Figure 15. It can be seen from Figure 15 that the increase of  will reduce the curvature of the curve, but can significantly increase the increment of the stress over time. Therefore, there is an optimal solution between  and  .

Genetic Algorithm Setting
Since the governing Equation (21) is a difference iterative form with time as the iterative variable, the least square method or linear regression method cannot be used, so the fitting problem can be transformed into a multi-objective optimization problem. The objective function, design variables and constraint conditions are prepared by using the GA toolbox of MATLAB 2016R. The optimization process is shown in Figure 16. It can be seen from Figure 15 that the increase of κ will reduce the curvature of the curve, but can significantly increase the increment of the stress over time. Therefore, there is an optimal solution between α and κ.

Genetic Algorithm Setting
Since the governing Equation (21) is a difference iterative form with time as the iterative variable, the least square method or linear regression method cannot be used, so the fitting problem can be transformed into a multi-objective optimization problem. The objective function, design variables and constraint conditions are prepared by using the GA toolbox of MATLAB 2016R. The optimization process is shown in Figure 16. n is the test data number. Population size: 200, elite count: 10, crossover fraction: 0.85, the end condition is that the two optimal fitness errors of individuals are less than 1 × 10 −15 .

Test Data Fitting
Based on the loading section of uniaxial cyclic tensile test data and genetic algorithm, there are nine sets of data corresponding to three strain rates and three strain amplitudes can fit out the different values of  and C. Then can take them into unloading program to get the loading and unloading cycle tensile curve.
The fitting results of fractional Maxwell model are shown in Figures 17-25. Some parameters are set as follows: y i stands stress value of the fitting curve corresponding to the test strain point, y i is the stress value of the test data, and n is the test data number.
Population size: 200, elite count: 10, crossover fraction: 0.85, the end condition is that the two optimal fitness errors of individuals are less than 1 × 10 −15 .

Test Data Fitting
Based on the loading section of uniaxial cyclic tensile test data and genetic algorithm, there are nine sets of data corresponding to three strain rates and three strain amplitudes can fit out the different values of α and C. Then can take them into unloading program to get the loading and unloading cycle tensile curve.
The fitting results of fractional Maxwell model are shown in Figures 17-25.   ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).  ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).  ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.). ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).  ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).  ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).  From Figures 17-25, it can be seen that the fractional Maxwell model has a good fit for the loading section, but the unloading section has a large deviation, and the overlap of loading and unloading curves cannot reflect hysteresis loop, this is because the maximum stress value of the loading section fitting curve is smaller. Table 5 shows the fitting coefficients of fractional Maxwell model corresponding to each group of data. It can be seen from Table 5 that each fitting iteration has a large number to ensure the credibility of the optimization and avoid falling into local optimization. The optimal fitness value (Fval) ranges from 0.79 to 3.98, which represents the approximation between the fitting curve and the test data, and the smaller Fval, the better the result gets. The values of  and C are not the same in the group 1-9. This phenomenon demonstrates that the damping capacity of M2052 is related to the strain rate and strain amplitude. At the same strain rate, with the increase of the strain amplitude, the combination  with C also makes the slope of the fitting curve larger. At the same strain amplitude, with the change of the strain rate, the combination  with C also change but no special law. The value range of  is generally small, while C value is bigger, in the range of 7811.4-29563. Combined with Section 4.1, since the max-stress of loading stage is too big, the fractional coefficient should be decreased and the  should be increased to achieve the minimum mean square error. This will lead to the weak nonlinear of fitting curve. It cannot conform to the characteristics of nonlinear damping alloy. Therefore, we need to modify model in Section 5. However, the fitting curve of the basic model is concentrated near the symmetry line of the test hysteresis loop, which indicates that the numerical range of the model fitting is close to the actual ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).

(%)
From Figures 17-25, it can be seen that the fractional Maxwell model has a good fit for the loading section, but the unloading section has a large deviation, and the overlap of loading and unloading curves cannot reflect hysteresis loop, this is because the maximum stress value of the loading section fitting curve is smaller. Table 5 shows the fitting coefficients of fractional Maxwell model corresponding to each group of data. It can be seen from Table 5 that each fitting iteration has a large number to ensure the credibility of the optimization and avoid falling into local optimization. The optimal fitness value (Fval) ranges from 0.79 to 3.98, which represents the approximation between the fitting curve and the test data, and the smaller Fval, the better the result gets. The values of α and C are not the same in the group 1-9. This phenomenon demonstrates that the damping capacity of M2052 is related to the strain rate and strain amplitude. At the same strain rate, with the increase of the strain amplitude, the combination α with C also makes the slope of the fitting curve larger. At the same strain amplitude, with the change of the strain rate, the combination α with C also change but no special law. The value range of α is 5.1945 × 10 −4 0.0125 generally small, while C value is bigger, in the range of 7811.4-29563. Combined with Section 4.1, since the max-stress of loading stage is too big, the fractional coefficient should be decreased and the κ should be increased to achieve the minimum mean square error. This will lead to the weak nonlinear of fitting curve. It cannot conform to the characteristics of nonlinear damping alloy. Therefore, we need to modify model in Section 5. However, the fitting curve of the basic model is concentrated near the symmetry line of the test hysteresis loop, which indicates that the numerical Materials 2020, 13, 2020 20 of 36 range of the model fitting is close to the actual value. It can be used as approximate curve or equivalent curve in engineering application or when the accuracy demand is not high.

Establishment of Correction Term
The fractional Maxwell model cannot simulate the hysteresis curve of damping alloys completely. Its maximum stress value of fitting curve is small. This is because of ignoring the friction between martensite and austenite or itself, that leads to an error term. So according to the difference value between the uniaxial cyclic tensile test under constant strain rate and fitting data of fractional Maxwell model, add a correction term. Figure 26 shows the error diagram of uniaxial cyclic tensile test data with the fitting curve of fractional Maxwell model in the loading stage. value. It can be used as approximate curve or equivalent curve in engineering application or when the accuracy demand is not high.

Establishment of Correction Term
The fractional Maxwell model cannot simulate the hysteresis curve of damping alloys completely. Its maximum stress value of fitting curve is small. This is because of ignoring the friction between martensite and austenite or itself, that leads to an error term. So according to the difference value between the uniaxial cyclic tensile test under constant strain rate and fitting data of fractional Maxwell model, add a correction term. Figure 26 shows the error diagram of uniaxial cyclic tensile test data with the fitting curve of fractional Maxwell model in the loading stage. As shown in Figure 26, it can be seen that the curve basically shows a sinusoidal wave peak. So the first-order sinusoidal function is used to fit, and the correction term can be obtained: Then, changing the original fractional Maxwell constitutive Equation (6) into (25), let the stress of original fractional Maxwell model is m  : The constitutive equation of modified fractional order Maxwell model can be obtained by adding the modified term (24):

Simulation Analysis of Fitting Data
Firstly, the correction coefficients of a, b and c were fitted in MATLAB according to Equation (24), and then the fitting curve of modified fractional order Maxwell model was obtained by substituting into Equation (26). The fitting diagram of the correction term is shown in Figure 27. The fitting coefficient and evaluation index of the correction term are shown in Table 6. As shown in Figure 26, it can be seen that the curve basically shows a sinusoidal wave peak. So the first-order sinusoidal function is used to fit, and the correction term can be obtained: Then, changing the original fractional Maxwell constitutive Equation (6) into (25), let the stress of original fractional Maxwell model is σ m : The constitutive equation of modified fractional order Maxwell model can be obtained by adding the modified term (24):

Simulation Analysis of Fitting Data
Firstly, the correction coefficients of a, b and c were fitted in MATLAB according to Equation (24), and then the fitting curve of modified fractional order Maxwell model was obtained by substituting into Equation (26). The fitting diagram of the correction term is shown in Figure 27. The fitting coefficient and evaluation index of the correction term are shown in Table 6. As shown in Figure 26, it can be seen that the curve basically shows a sinusoidal wave peak. So the first-order sinusoidal function is used to fit, and the correction term can be obtained: Then, changing the original fractional Maxwell constitutive Equation (6) into (25), let the stress of original fractional Maxwell model is m  : The constitutive equation of modified fractional order Maxwell model can be obtained by adding the modified term (24):

Simulation Analysis of Fitting Data
Firstly, the correction coefficients of a, b and c were fitted in MATLAB according to Equation (24), and then the fitting curve of modified fractional order Maxwell model was obtained by substituting into Equation (26). The fitting diagram of the correction term is shown in Figure 27. The fitting coefficient and evaluation index of the correction term are shown in Table 6.  It can be seen from Figure 27 and Table 6 that the correction term fits the deviation value well, which basically reflects the change of the deviation value. Root of mean square (RMSE) values are all less than 0.74, and some coefficient of determination (R 2 ) values are smaller, which is caused by too many points of the fitting samples.   It can be seen from Figure 27 and Table 6 that the correction term fits the deviation value well, which basically reflects the change of the deviation value. Root of mean square (RMSE) values are all less than 0.74, and some coefficient of determination (R 2 ) values are smaller, which is caused by too many points of the fitting samples. Figures 28-36  ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.). ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.). ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).  ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).  ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.). ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.). ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.). ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.). From Figures 28-36 we can see modified fractional Maxwell model fitting is better with the test data, and can clearly show the nonlinear constitutive relation of damping alloy. The hysteresis area is obvious, and the fitting curve is smoother than test data. This is convenient to be used in analysis. The loading and unloading curve of modified fractional Maxwell model fitting are symmetry about the center line. However, the stress-strain fitting curve at the strain rate of 0.01%/s and the strain ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.). From Figures 28-36 we can see modified fractional Maxwell model fitting is better with the test data, and can clearly show the nonlinear constitutive relation of damping alloy. The hysteresis area is obvious, and the fitting curve is smoother than test data. This is convenient to be used in analysis. The loading and unloading curve of modified fractional Maxwell model fitting are symmetry about the center line. However, the stress-strain fitting curve at the strain rate of 0.01%/s and the strain amplitude of 0.15% has a small hysteresis area, and the loading section is well fitted, but the unloading section has a large deviation, which is due to the error of experimental data collection. Table 7 is the evaluation index of modified fractional Maxwell model corresponding to each group of data. It can be seen from Table 7 that the modified fractional Maxwell model has a good fitting effect on the experimental data, and the determination coefficients can reach above 0.99, but the mean square error and variance of group 9 are large. Since its fitting deviation is mainly in the unloading section, try to set the unloading period of experiment data as object of correction fitting program. First, we can use the coefficient and the model of the fractional Maxwell to calculate the fitting curve of unloading stage. Then we can use it and the unloading stage data of the test to fit the modified term. Finally, the loading stage was calculated according to the modified fractional Maxwell model. The unloading deviation curve and the unloading correction term fitting curve were shown in Figure 37a The fitting curve in Figure 38 compared with the test data concluded that MSE = 1.363, SSE = 490.673, R 2 = 0.9990, also it can be seen from the picture that effect is superior to the original loading correction fractional fitting. This reason should be caused by the error generated in test collection, suggesting non-linear of unloading stage is better. This has confirmed that modified fractional Maxwell for this strain rate and strain amplitude is feasible. ε(%/s) is strain rate, ε a (%) is strain amplitude, the fitting curve is red and the experiment data is black.).

(%)
The fitting curve in Figure 38 compared with the test data concluded that MSE = 1.363, SSE = 490.673, R 2 = 0.9990, also it can be seen from the picture that effect is superior to the original loading correction fractional fitting. This reason should be caused by the error generated in test collection, suggesting non-linear of unloading stage is better. This has confirmed that modified fractional Maxwell for this strain rate and strain amplitude is feasible.
In order to show the applicability of modified fractional Maxwell model, a comparison of fitting performance between the modified fractional Maxwell model and other models was be done. Choosing the stress-strain curve of experiment at 0.005%/s and 0.1% as the fitting standard.
According to Boltzmann superposition principle [34], the stress of material can be express by Equation (27).
where ε 0 is the Initial stain value, G(t) is the relaxation modulus, . ε is the strain rate, " * " is the symbol of Stieltjes convolution [34].
For the classical Maxwell model (Maxwell two-parameter model), as shown in Figure 39.
where E 1 is the elastic modulus, τ 1 = In order to show the applicability of modified fractional Maxwell model, a comparison of fitting performance between the modified fractional Maxwell model and other models was be done. Choosing the stress-strain curve of experiment at 0.005%/s and 0.1% as the fitting standard.
According to Boltzmann superposition principle [34], the stress of material can be express by Equation (27).
where 0  is the Initial stain value, () Gt is the relaxation modulus,  is the strain rate, " * " is the symbol of Stieltjes convolution [34].  ε is constant, take Equation (28) into (27) can get: For the Maxwell three-parameter model, as shown in Figure 40.
where E 1 and E 2 are the elastic modulus, τ 1 = η 1 E 1 is the relaxation time.
where E1 is the elastic modulus, Since the  is constant, take Equation (28) For the Maxwell three-parameter model, as shown in Figure 40.  ε is constant, take Equation (30) into (27) can get: For the Maxwell four-parameter model, as shown in Figure 41.
where E 1 and E 2 are the elastic modulus, τ 1 = where E1 and E2 are the elastic modulus, Since the  is constant, take Equation (30) For the Maxwell four-parameter model, as shown in Figure 41.
where E1 and E2 are the elastic modulus,  ε is constant, take Equation (32) into (27) can get: Next, least-squares approximations are used, and two-, three-and four-parameter optimal models were determined. The parameters and evaluation index of Maxwell models are given in Table 8. Figure 42 presents experiment data and the optimal Maxwell models.  It is easy to observe, that the classical Maxwell model, Maxwell three-parameter model or Maxwell four-parameter model are inappropriate for description of this nonlinearity of elastic loading process. A better fit to experimental data can be obtained, if the modified fractional Maxwell model is used. In addition to, because the E1 of Maxwell three-parameter model is too small that is almost equal to zero, the fitting curve of Maxwell two-parameter model is close to Maxwell threeparameter model shown in the Figure 39 and Table 8. All this has confirmed that the modified fractional Maxwell can well explain the nonlinear constitutive relation of damping alloy.
Since the parameters of the modified Maxwell model need to be determined according to the different conditions of each load, the application of this model is limited. In order to improve the usability of the model, the relationship between model parameters and loading conditions is  It is easy to observe, that the classical Maxwell model, Maxwell three-parameter model or Maxwell four-parameter model are inappropriate for description of this nonlinearity of elastic loading process. A better fit to experimental data can be obtained, if the modified fractional Maxwell model is used. In addition to, because the E1 of Maxwell three-parameter model is too small that is almost equal to zero, the fitting curve of Maxwell two-parameter model is close to Maxwell three-parameter model shown in the Figure 39 and Table 8. All this has confirmed that the modified fractional Maxwell can well explain the nonlinear constitutive relation of damping alloy.
Since the parameters of the modified Maxwell model need to be determined according to the different conditions of each load, the application of this model is limited. In order to improve the usability of the model, the relationship between model parameters and loading conditions is analyzed. By using the method, which calculated the average value of the parameters at different strain rates but same strain amplitude, the influence of strain rate was eliminated, and the variation rule between strain amplitude and each parameter was obtained. Then, a strain-related formula is proposed for each parameter, so that the parameters of the modified Maxwell model can be determined and the model can maintain at a high precision under the condition that the particular experiment cannot be carried out.
However, the first step is to assess the impact of replacing the original parameter with the mean value of parameter. According to the values of parameters under different strain rates and strain amplitudes, the parameters under the same strain amplitudes but different strain rates were averaged, as shown in Table 9. This way can minimize the impact of strain rates on model parameters, and get the functions of model parameters, which are only related to strain. For evaluating the impact of the average value of each parameter on the accuracy of model fitting, the average value (Table 10) was used to replace the original model parameters, and substituted into MATLAB. The fitting effect is expressed in figures of error, as shown in Table 11.  As can be seen from Table 11, the errors of the mean parameters fitting relative to the test and the original parameters fitting are all within the acceptable range. The R 2 values between mean parameters fitting curve and test are all above 0.9954, while the R 2 values between mean parameters fitting curve and the original parameters fitting curve are all above 0.9996. Therefore, the influence of the strain rate can be ignored and the mean parameters can be used to replace the original parameters. Then, the function of each parameter can be proposed according to the relationship between the mean parameter and the strain, as shown in Figures 43-47.  As can be seen from Table 11, the errors of the mean parameters fitting relative to the test and the original parameters fitting are all within the acceptable range. The R 2 values between mean parameters fitting curve and test are all above 0.9954, while the R 2 values between mean parameters fitting curve and the original parameters fitting curve are all above 0.9996. Therefore, the influence of the strain rate can be ignored and the mean parameters can be used to replace the original parameters. Then, the function of each parameter can be proposed according to the relationship between the mean parameter and the strain, as shown in Figures 43-47.      The fitting function of the relation between each parameter and strain is shown in Table 12.   Table 12 all have a good fitting result for the change of each parameter with strain. What is interesting is that all the functions of parameters can be expressed as exponential functions, except for α. Therefore, the functions in Table 12 can be used to calculate the values of model parameters under other loading conditions within the elastic range, instead of the constant strain rate tensile test, which extends the practicability of the modified Maxwell model.

Conclusions
In view of the existing study about nonlinear constitutive relation of Mn-Cu damping alloy being less, and that it should not be treated as linear elastic material, this study chose M2052 damping alloy as the research object, through uniaxial cyclic tensile test under constant strain rate to analysis its nonlinear constitutive relation and hysteretic characteristics. By considering the damping alloy as a special viscoelastic material and based on the fractional Maxwell model which describes nonlinear viscoelasticity, a modified fractional Maxwell model suitable for M2052 damping alloy was proposed. Through the numerical simulation, we can get the following conclusions: 1. Through uniaxial cyclic tensile test with constant strain rate, it is concluded that Mn-Cu damping alloy can be considered as a viscoelastic material. In the elastic strain range, as the strain amplitude increases, the slope of the stress-strain curve decreases and the hysteresis loop area and the damping capacity increases with the same strain rate. Under the same strain amplitude, the slope and hysteresis area of the stress-strain curve change with the change of the strain rate, this still needs further study.
2. The Fractional Maxwell model only contains two unknown coefficients and does not have to get analytic solutions. It can combine with uniaxial cyclic tensile test under constant strain rate and relaxation test to determine its parameters. Since the max-stress of loading stage is biggish, the fractional coefficient should be decreased and the quasi-state coefficient should be increased to achieve the minimum mean square error. This leads to weak nonlinear fitting curve, which does not match the nonlinear characteristics of damping alloy. However, the fitting curve of fractional Maxwell model is basically in the middle of the loading and unloading stress-strain curve, which can be used in some engineering application or situations where the accuracy is not high.
3. The error of fractional Maxwell model is due to ignoring the friction of new martensite with austenite and original martensite. Combining with the deviation curve, a five-variable modified fractional Maxwell phenomenological model is proposed. This model can well simulate the nonlinear characteristics and the hysteresis curve of damping alloys. Determine coefficients all above 0.99, and the deviation is small. Data fluctuation and measurement error exist in the test curve. However, the fitting curves do not have these problems and can fully capture the stress-strain variation trend of the damping alloy. Through compare with other Maxwell models, it also shows a better fitting performance. In addition, a strain-related formula is proposed for each parameter, so that the parameters of the modified Maxwell model can be determined without the particular experiment. Therefore, the modified fractional Maxwell model can be used as the constructive equation of Mn-Cu damping alloys and reference for the further analysis.
Supplementary Materials: The following are available online at http://www.mdpi.com/1996-1944/13/9/2020/s1. Figure S1: Hysteresis measured curves with different strain amplitudes at strain rate of 0.0025%/s (σ is stress, ε is strain, ε a is strain amplitude); Figure S2: Hysteresis measured curves with different strain amplitudes at strain rate of 0.005%/s (σ is stress, ε is strain, ε a is strain amplitude); Figure S3: Hysteresis measured curves with different strain amplitudes at strain rate of 0.01%/s (σ is stress, ε is strain, ε a is strain amplitude); Table S1: Hysteresis area of measured under different strain rates and strain amplitudes (unit: 10 kJ/m 3 ); Table S2: The slope of measured curve under different strain rates and strain amplitudes.