Dynamic characteristics of Mn-Cu high damping alloy subjected to impact load

Abstract

Mn-Cu high damping alloy is a twin-type damping alloy. Owing to its martensite twin structure at room temperature, it can convert vibration energy into heat energy, thereby reducing vibration. Although essentially a nonlinear elastic material, Mn-Cu damping alloys are treated as linear elastic materials in current engineering practice. However, introducing a constant damping coefficient alone will produce significant errors when modeling vibration reduction characteristics of the material, especially under impact loading. In this study, vibration test was performed on Mn-Cu damping alloy cantilever beam subjected to an impact load and deviation between the test result and the one of existing modeling method was analyzed. A generalized fractional-order Maxwell model was established to describe the nonlinear constitutive relation of the Mn-Cu damping alloy. Then, the model was extended to the three-dimensional state and a secondary development was performed. Finally, various applications for the damping alloy were explored and the effects of the damping alloy on the vibration characteristics of composite cantilever beam structures subjected to impact loads were investigated with the aim of better understanding the dynamic characteristics and improving the effectiveness of vibration reduction applications using Mn-Cu damping alloy in the future.

Introduction

The Mn-Cu high damping alloy is a typical twin-type damping alloy and its main elements are Mn, Cu, Al, Fe, and Ni. Commercially available Mn-Cu based damping alloys include Sonoston, ¹ Incramute, ² ABPOPA, ³ and M2052 (Mn-20Cu-5Ni-2Fe). ⁴ Among them, the M2052 alloy is characterized by high damping characteristics, high strength (comparable to structural steel), high ductility, and fine processing performance, and therefore, offers the best performance. In addition, M2052 can be directly used as a raw material in mechanical parts and has been extensively applied in vibration reduction applications for machinery, particularly in the military industry. Wang et al. ⁵ used the Mn-Cu damping alloy in a vibration reduction component for a high-speed rotary fuse-testing machine and tested its dynamic characteristics. The novel Mn-Cu damping alloy was found to have a high stiffness coefficient and effectively reduced the vibration amplitude under high-frequency service conditions and can be used for high-speed rotary equipment. Xin and Mao ⁶ applied Mn-Cu damping alloy in a remote control weapon station. A simulation analysis was performed on a semi-gear structure and three different vibration reduction schemes were designed. By conducting acceleration attenuation tests on a real semi-gear structure containing the damping alloy, it was concluded that the damping gears effectively absorb vibrations. Wang et al. ⁷ designed a tribological system with a damping alloy component and demonstrated through experiments and simulations that damping alloy components with a grooved surface can significantly improve vibration suppression, noise reduction, and wear resistance of the system. However, in all of the above studies, the Mn-Cu damping alloy was treated as a linear elastic material and the damping characteristics were modeled by introducing a constant damping coefficient, which does not reflect the nonlinear constitutive characteristics of the Mn-Cu damping alloy.

The high damping characteristics of Mn-Cu damping alloys are derived from two sources: First, during high-temperature melting and cooling, above the Néel temperature, the paramagnetism is transformed into anti-ferromagnetism and anti-ferromagnetic martensite twins are generated; Second, the temperature for producing martensite is higher than room temperature, therefore, stress-assisted martensite is generated under the action of the external force, resulting in interface slip between martensite and austenite parent phases.^8–10 These special microstructures and motions absorb vibration energy, resulting in highly nonlinear behavior that is dependent upon temperature, frequency, strain, strain rate, and other factors. Thus, assuming linear behavior limits the application of these models. In general, damping alloys are characterized by an elliptic hysteresis loop, accompanied by more complex nonlinear damping. Many studies on Mn-Cu damping alloys have reported that the damping characteristics are dependent on the vibration amplitude and frequency.^9,11,12 Generally, Mn-Cu damping alloys are used in structural components that do not undergo plastically deformation. Nonetheless, the Mn-Cu damping alloy is a nonlinear elastic material and it behavior differs from steel and other linear elastic materials.

To accurately analyze the dynamic characteristics of the Mn-Cu damping alloy subjected to an impact load, the first step is to build a model to describe its nonlinear constitutive relation. To date, only a few studies have investigated the nonlinear constitutive relation of the Mn-Cu damping alloy, which have mainly focused on shape memory alloys and other types of damping alloys. Satisfactory results can be achieved by treating damping alloys as a special type of viscoelastic materials and describing the constitutive relation using the theory of viscoelasticity.^13–15

Viscoelastic constitutive models mainly include the Maxwell model, Kelvin model, generalized Maxwell model, and generalized Kelvin model, which are all linear viscoelastic models and are not capable of fully expressing the constitutive relation of nonlinear viscoelastic materials. With the introduction of fractional derivatives, the fractional-order Maxwell model and the fractional-order Kelvin model emerged, making it possible to describe the constitutive relation of nonlinear viscoelastic materials. The normal form of the generalized fractional-order Maxwell model combines multiple fractional-order Maxwell models in parallel, thus further improving the ability to describe nonlinear materials. This has performed well in describing the constitutive relation of nonlinear viscoelastic materials. Sun and Chen ¹⁶ described a general linear viscoelastic damper using the generalized fractional-order Maxwell model. The generalized relation between force and deformation of the viscoelastic damper (material) can be depicted within a wide frequency range. In addition, an explicit formula to determine the maximum attainable damping and corresponding optimal model parameters were derived. Costa and Ribeiro ¹⁷ tested the dynamic mechanical properties of a carbon fiber reinforced polymer. By fitting the results, a generalized fractional-order viscoelastic model with high correlation and numerical stability was obtained. Based on the generalized fractional-order Maxwell model, Fang ¹⁸ proposed a nonlinear constitutive model for a thermotropic shape memory polymer, and described its static viscoelastic behavior in detail. However, the fractional-order Maxwell models described above mainly describe the behavior of viscoelastic fluids and cannot be directly applied to viscoelastic solids. Generalized fractional-order Maxwell models for viscoelastic solids are rarely reported.

Current practice for designing structures used in complex loading environments is to first build a finite element model. The finite element software ANSYS contains a variety of material models and elements that are capable of solving a variety of problems related to nonlinear materials. However, a constitutive model for describing the Mn-Cu damping alloy is not currently available. While the stress relation of damping alloys during loading and unloading is assumed to be linear, hysteresis phenomena occur in the loading and unloading processes. Thus, the constitutive model in ANSYS cannot be used to effectively simulate the mechanical properties of Mn-Cu damping alloy and it is necessary to develop a suitable model.

In this study, to investigate the response of the Mn-Cu damping alloy subjected to an impact load, a cantilever beam was first considered and vibration tests were performed on Mn-Cu damping alloy cantilever beam. Deviations in the results of existing approaches using linear constitutive model and constant damping coefficient from the test results were analyzed. Then, a constitutive equation for the damping alloy was established based on the generalized fractional-order Maxwell model and a secondary development of the constitutive model was performed in ANSYS. Finally, eight different applications for the damping alloy were proposed and their responses to an impact load were investigated, thus providing a theoretical reference for subsequent design of structures incorporating the damping alloy.

Vibration tests of Mn-Cu damping alloy cantilever beam subjected to impact load

Cantilever beam structures are basic structures commonly used in the fields of structural and mechanical engineering. The simple form is convenient for analysis and responses of damping alloy cantilever beam structures subjected to an impact load are representative, to some extent. The effect of temperature on the damping alloy was not considered and all tests were performed at room temperature.

Test material and equipment

A thermally treated M2052 damping alloy plate, mainly composed of Mn–22.1Cu-5.24Ni–1.93Fe, was adopted. Properties of the M2052 damping alloy are listed in Table 1. Spark cutting was used to construct 250 mm × 20 mm × 5 mm cantilever beams from the alloy plate. The length of each cantilever beam (outside of the fixture) was 200 mm and an impact hammer was used for excitation. Due to the large vibration created by the impact load, the accelerometer was fixed to the sample via ∅5 mm mounting holes. The dimensions of the cantilever beams are illustrated in Figure 1. A representative photograph of the cantilever beam is presented in Figure 2.

Table 1.

Properties of M2052 damping alloy.

Name	Value
Tensile strength	540 MPa
Elastic modulus	68.5 GPa
Poisson’s ratio	0.338
Shrinkage rate	50.5%
Yield strength	205 MPa
Density	7250 kg/m3
Elongation at break	31.5%
Brinell hardness	138 HBW

Figure 1.

Dimensions of cantilever beam.

Figure 2.

Photograph of Mn-Cu damping alloy cantilever beam.

The main test equipment were provided by Donghua Testing Technology Co., Ltd., including a DHVTC vibration test and control bench, impact hammer, accelerometer, and data collector. Details of the specific equipment models and parameters are provided in Table 2. The sampling frequency was 2 kHz, the data collection mode was 2-channel data collection, and an integral electronics piezoelectric (IEPE) power supply was adopted.

Table 2.

Test equipment.

Name	Model	Sensitivity	Range
Impact hammer	LC02-3A102E	4.36 Pc/N	1000 N
Accelerometer	1A101	4.945 mV/m·s−2	100 g
Data collector	DH5922N	–	–

One end of the cantilever beam was fixed to the bench using a fixture and a length of 200 mm was reserved. The accelerometer was bolted to the test bench via mounting holes. The equipment set up is shown in Figure 3.

Figure 3.

Equipment set up.

Testing protocol and analysis of results

One end of the cantilever beam was fixed and an impact load of 100 N was applied to the other end using the impact hammer. The free-end vibration acceleration of the cantilever beam was collected using the accelerometer and quadratic integration was performed to obtain vibration displacement-time curves. Since it is impossible to guarantee load stability with an impact hammer, the test was repeated ten times and average values were calculated. The force spectrum and test results are presented in Figures 4 and 5, respectively.

Figure 4.

Force spectrum of impact hammer.

Figure 5.

Free-end vibration displacement of Mn-Cu damping alloy beam subjected to 100 N impact load.

As seen in Figure 5, the Mn-Cu damping alloy has a vibration amplitude of 1.352 mm under an impact load of 100 N and regains stability within approximately 0.4 s. To determine the error produced using the ANSYS constitutive model with a constant damping coefficient for calculating the vibration displacement, a finite element model was established for the Mn-Cu damping alloy cantilever beam. One end of the beam was modeled as a fixed constraint and an impact load of 100 N was applied to the other end of the beam. The duration of load application was set to 0.002 s according to the force spectrum. The finite element model of the beam is shown in Figure 6.

Figure 6.

Finite element model of Mn-Cu damping alloy cantilever beam.

The simulation time step was 0.001 s with a total duration of 1 s. The constant damping coefficient was 0.06. Material properties used in the simulation are presented in Table 1. The mode superposition method was adopted. A comparison of the simulation results and test results is presented in Table 7.

As seen in Figure 7, the vibration displacement amplitude obtained using the linear elastic constitutive model with a constant damping coefficient is far larger than the real value for the Mn-Cu damping alloy. The vibration displacement amplitude obtained from simulations using ANSYS linear constitutive model was 5.864 mm, with an error of 334%. However, at 0.4 s, the vibration displacement was close to 0, and attenuation of the vibration amplitude was greater than the real value because although the Mn-Cu damping alloy is treated as a liner elastic material in ANSYS, its elastic modulus is only 1/3 that of steel. According to the theory of linear elasticity, structures with a smaller elastic modulus will undergo larger deformation when loaded. Thus, under an impact load of 100 N, the simulated vibration displacement amplitude was much higher than the real value. Furthermore, since a relatively large damping coefficient was set in the simulations, attenuation of the vibration amplitude was also larger compared to the real Mn-Cu damping alloy. This will affect the accuracy of design and application of the Mn-Cu damping alloy. Neglecting nonlinear constitutive relation of Mn-Cu damping alloy is the root cause of the discrepancies.

Figure 7.

Comparison of test value of vibration displacement of cantilever beam subjected to 100 N impact load and simulation value obtained by linear elastic constitutive model.

Constitutive model of Mn-Cu damping alloy based on generalized fractional-order Maxwell model

Development of generalized fractional-order Maxwell model

According to the mechanism of high damping characteristics of the Mn-Cu damping alloy, relaxation of the twin structure and friction motion of the stress-assisted martensite are assumed to be between those of elastic motion and viscous motion. The basic form of the Maxwell model was adopted and used to connect one spring-pot element and one spring element in series to create the fractional-order Maxwell model of the twin term and the martensite friction term.

At the microscopic level, twin motion and stress-assisted martensite phase change are both lattice motion deformations driven by an external force; therefore, the strains should be the same and the total stress is the difference between the external force and the elastic force. Two fractional-order Maxwell terms, representing the twin damping term and the martensite friction term, were connected in parallel to express the nonlinear constitutive relation of the M2052 damping alloy. The proposed generalized fractional-order Maxwell model is illustrated in Figure 8.

Figure 8.

Generalized fractional-order Maxwell model.

In Figure 8, ${\sigma }_1$ and ${\sigma }_3$ are the stresses of the spring-pot element; ${\epsilon }_1$ and ${\epsilon }_3$ are the strains of the spring-pot element; ${\sigma }_2$ and ${\sigma }_4$ are the stresses of the spring element; ${\epsilon }_2$ and ${\epsilon }_4$ are the strains of the spring element; $\kappa$ and $\phi$ are the quasi-state characteristics of the spring-pot element; ¹⁹ $\alpha$ and $\beta$ are the fractional-order coefficients, $0\le \alpha ,\beta \le 1$ ; E is the Young’s modulus of the M2052 damping alloy; $\sigma$ is the total stress borne by the damping alloy; $\epsilon$ is the total strain borne by the damping alloy.

Combining the generalized Maxwell model relation and constitutive relation of the fractional-order spring pot ²⁰ yields

$$\epsilon =\frac{1}{\kappa }{D}_t^{-\alpha }{\sigma }_1+\frac{{\sigma }_1}{E}=\frac{1}{\phi }{D}_t^{-\beta }{\sigma }_3+\frac{{\sigma }_3}{E}$$

(1)

where $D_t$ is the fractional-order differential operator. ²¹ Assuming ${\sigma }_1=w\sigma$ and ${\sigma }_3=(1-w)\sigma$ , the constitutive equation of the generalized Maxwell model can be obtained as

$$\epsilon =\frac{w}{\kappa }{D}_t^{-\alpha }\sigma +\frac{w\sigma }{E}=\frac{1-w}{\phi }{D}_t^{-\beta }\sigma +(1-w)\frac{\sigma }{E}$$

(2)

Control equation

From a conservation of energy perspective, the work done by the external force that is converted into strain energy at a constant temperature is the sum of the elastic potential energy, twin energy, and martensite-austenite and martensite-martensite friction energies, expressed as

$$\begin{array}{c}\multicolumn{1}{c}{W_{\mathrm{strain}\mathrm{energy}}=W_{\mathrm{elastic}\mathrm{potential}\mathrm{energy}}}\\ +W_{\mathrm{twin}\mathrm{energy}}+W_{\mathrm{martensite}\mathrm{friction}\mathrm{energ}\mathrm{y}}\hfill \end{array}$$

(3)

According to J2 deformation theory and Hooke’s law, the energy per unit volume in each case can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{W_{\mathrm{strain}\mathrm{energy}}={\int }_0^t\sigma \stackrel{·}{\epsilon }\mathrm{dt},W_{\mathrm{elastic}\mathrm{potential}\mathrm{energy}}}\\ =\frac{1}{2}E{\left({\epsilon }_2\right)}^2+\frac{1}{2}E{\left({\epsilon }_4\right)}^2\hfill \end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{W_{\mathrm{twin}\mathrm{energy}}={\int }_0^t{\sigma }_1{\stackrel{·}{\epsilon }}_1\mathrm{dt},W_{\mathrm{martensite}\mathrm{friction}\mathrm{energ}\mathrm{y}}}\\ ={\int }_0^t{{\sigma }_3\stackrel{·}{\epsilon }}_3\mathrm{dt},\hfill \end{array}$$

where $\stackrel{·}{\epsilon }$ is the strain rate; ${\stackrel{·}{\epsilon }}_1$ is the strain rate of the twin term; ${\stackrel{·}{\epsilon }}_3$ is the strain rate of the martensite friction term.

Substituting the above equations into Formula (3) yields the following expression of the control equation:

$$\begin{array}{c}{\int }_0^t{\sigma }_1\stackrel{·}{\epsilon }\mathrm{dt}+{\int }_0^t{\sigma }_3\stackrel{·}{\epsilon }\mathrm{dt}=\frac{1}{2}\frac{{\sigma }_1^2}{E}+\frac{1}{\kappa }{\int }_0^t{\sigma }_1{D}_t^{1-\alpha }{\sigma }_1\mathrm{dt}\hfill \\ +\frac{1}{2}\frac{{\sigma }_3^2}{E}+\frac{1}{\phi }{\int }_0^t{\sigma }_3{D}_t^{1-\beta }{\sigma }_3\mathrm{dt}\hfill \end{array}$$

(4)

Formula (4) can be divided into

$${\int }_0^t{\sigma }_1\stackrel{·}{\epsilon }\mathrm{dt}=\frac{1}{2}\frac{{\sigma }_1^2}{E}+\frac{1}{\kappa }{\int }_0^t{\sigma }_1{D}_t^{1-\alpha }{\sigma }_1\mathrm{dt}$$

(5)

$${\int }_0^t{\sigma }_3\stackrel{·}{\epsilon }\mathrm{dt}=\frac{1}{2}\frac{{\sigma }_3^2}{E}+\frac{1}{\phi }{\int }_0^t{\sigma }_3{D}_t^{1-\beta }{\sigma }_3\mathrm{dt}$$

(6)

When $t\le 0$ , $\sigma (t)=0,\epsilon (t)=0$ ; when $t>0$ , $\sigma (t)$ , and $\epsilon (t)$ exhibit a monotonic increase during loading and a monotonic decrease during unloading. In this case, the initial value and boundary conditions of the control equation can be defined as follows:

$$\{\begin{array}{c}\sigma (0)=0,t>0\hfill \\ \multicolumn{1}{c}{\epsilon (0)=0,t>0}\\ \multicolumn{1}{c}{\sigma (t_{i+1})>\sigma (t_i),\stackrel{·}{\epsilon }>0}\\ \multicolumn{1}{c}{\sigma (t_{i+1})<\sigma (t_i),\stackrel{·}{\epsilon }<0}\end{array}$$

(7)

Numerical solution

By taking the derivatives of both sides of the control equations, equations (5) and (6), with respect to time, we obtain

$$\stackrel{·}{\epsilon }=\frac{{\stackrel{·}{\sigma }}_1}{E}+\frac{1}{k}{D}_t^{1-\alpha }{\sigma }_1$$

(8)

$$\stackrel{·}{\epsilon }=\frac{{\stackrel{·}{\sigma }}_3}{E}+\frac{1}{\phi }{D}_t^{1-\beta }{\sigma }_3$$

(9)

To solve the equations, the fractional derivative is discretized using the finite difference method.

After defining $t_m=m\mathrm{\Delta }t$ and $m=0,1,2,...,K$ , we assume $0\le k\le K$ where $\mathrm{\Delta }t=\frac{T}{K}$ is the step size. According to Caputo’s definition of derivatives, ²² the finite difference of the fractional differential of stress can be obtained as

$$\begin{array}{c}{D}_t^h\sigma (t_k)=\frac{{\left(\mathrm{\Delta }t\right)}^{-h}}{\mathrm{\Gamma }(2-h)}\left[\begin{array}{c}\sigma (t_k)-h_k\sigma (t_0)+\hfill \\ \multicolumn{1}{c}{\sum _{j=2}^k(h_j-h_{j-1})\sigma (t_{k-j+1})}\end{array}\right]\hfill \\ +o((\mathrm{\Delta }t{)}^{2-h})\hfill \end{array}$$

(10)

where h is the fractional-order coefficient, $h_j=j^{1-h}-(j-1{)}^{1-h}$ , $j=2,...,k$ , $o((\mathrm{\Delta }t{)}^{2-h})$ is the higher-order error term.

When $0<\alpha <1$ , $0<\beta <1$ , then $1-\alpha >0$ and $1-\beta >0$ . Substituting $h=1-\alpha$ and $h\prime =1-\beta$ into equation (10) yields

$$\begin{array}{c}{D}_t^{1-\alpha }\sigma (t_k)=\frac{{\left(\mathrm{\Delta }t\right)}^{\alpha -1}}{\mathrm{\Gamma }(1+\alpha )}\left\{\begin{array}{c}\sigma (t_k)-h_k\sigma (t_0)+\hfill \\ \multicolumn{1}{c}{\sum _{j=2}^k(h_j-h_{j-1})\sigma (t_{k-j+1})}\end{array}\right\}\hfill \\ +o((\mathrm{\Delta }t{)}^{1+\alpha })\hfill \end{array}{h}_j=j^{\alpha }-(j-1{)}^{\alpha },j=2,...,k$$

(11)

$$\begin{array}{c}{D}_t^{1-\beta }\sigma (t_k)=\frac{{\left(\mathrm{\Delta }t\right)}^{\beta -1}}{\mathrm{\Gamma }(1+\beta )}\left\{\begin{array}{c}\sigma (t_k)-h_k^{\prime }\sigma (t_0)+\hfill \\ \multicolumn{1}{c}{\sum _{j=2}^k(h_j^{\prime }-h_{j-1}^{\prime })\sigma (t_{k-j+1})}\end{array}\right\}\hfill \\ +o((\mathrm{\Delta }t{)}^{1+\beta })\hfill \end{array}{h}_j^{\prime }=j^{\beta }-(j-1{)}^{\beta },j=2,...,k$$

(12)

By substituting equations (11) and (12) into equations (8) and (9) and omitting the higher-order error terms, we obtain the finite difference expression for strain rate with respect to twin stress and friction-induced martensite stress:

$$\begin{array}{c}\stackrel{·}{\epsilon }=\frac{{\sigma }_1(t_{k+1})-{\sigma }_1(t_k)}{E\mathrm{\Delta }t}\hfill \\ +\frac{1}{\kappa }\frac{{\left(\mathrm{\Delta }t\right)}^{\alpha -1}}{\mathrm{\Gamma }(1+\alpha )}\left\{\begin{array}{c}{\sigma }_1(t_k)-h_k{\sigma }_1(t_0)+\hfill \\ \multicolumn{1}{c}{\sum _{j=2}^k(h_j-h_{j-1}){\sigma }_1(t_{k-j+1})}\end{array}\right\}\hfill \\ (k=1:m,h_j=j^{\alpha }-(j-1{)}^{\alpha },j=2,...,k)\hfill \end{array}$$

(13)

$$\begin{array}{c}\stackrel{·}{\epsilon }=\frac{{\sigma }_3(t_{k+1})-{\sigma }_3(t_k)}{E\mathrm{\Delta }t}\hfill \\ +\frac{1}{\phi }\frac{{\left(\mathrm{\Delta }t\right)}^{\beta -1}}{\mathrm{\Gamma }(1+\beta )}\left\{\begin{array}{c}{\sigma }_3(t_k)-h_k^{\prime }{\sigma }_3(t_0)+\hfill \\ \multicolumn{1}{c}{\sum _{j=2}^k(h_j^{\prime }-h_{j-1}^{\prime }){\sigma }_3(t_{k-j+1})}\end{array}\right\}\hfill \end{array}ç(k=1:m,h{\prime }_j=j^{\beta }-(j-1{)}^{\beta },j=2,...,k)$$

(14)

Through a simplified derivation of Formula (2), equations (13) and (14) can be simplified as

$$\begin{array}{c}\stackrel{·}{\epsilon }=\frac{w\sigma (t_{k+1})-w\sigma (t_k)}{E\mathrm{\Delta }t}\hfill \\ +\frac{1}{\kappa }\frac{{\left(\mathrm{\Delta }t\right)}^{\alpha -1}}{\mathrm{\Gamma }(1+\alpha )}\left\{\begin{array}{c}w\sigma (t_k)-h_{k}w\sigma (t_0)\hfill \\ \multicolumn{1}{c}{+\sum _{j=2}^k(h_j-h_{j-1})w\sigma (t_{k-j+1})}\end{array}\right\}\hfill \\ (k=1:m,h_j=j^{\alpha }-(j-1{)}^{\alpha },j=2,...,k)\hfill \end{array}$$

(15)

The results of the uniaxial cyclic tensile tests at constant strain rate can be fitted to equation (15) to define various unknown parameters.

Simulation analysis and verification of numerical model

The Mn-Cu damping alloy plate was cut into 150 mm × 20 mm × 2 mm samples and the cyclic tensile tests were conducted on an MTS-Landmark 810 universal testing machine at various strain rates and at the maximum strain amplitude, according to a previously described protocol. ²³ Strain was measured with a 50 mm extensometer at a sampling frequency of 10 Hz. The samples were stretched at three fixed strain rates of 0.0025%/s, 0.005%/s, and 0.01%/s up to maximum strain amplitudes of 0.05%, 0.1%, and 0.15%, respectively, and the corresponding stress-strain curves were obtained. The stress-strain curves served as a basis for the simulation analysis and were also used to verify the numerical model. Properties of the M2052 damping alloy are provided in Table 1.

Based on data from the uniaxial cyclic tensile tests at a constant strain rate, the genetic algorithm was adopted to determine the objective function, design variables, and constraint conditions. For instance, the range of the fractional-order coefficient was $0<\alpha <1$ ; ranges of the quasi-state coefficients were $\kappa =C\times {10}^9$ and $68.5<C<\infty$ ; and the range of the weight coefficient w was $0\le w\le 1$ . The test data were fitted, as shown in Figures 9 to 11.

Figure 9.

Stress–strain curves at various strain amplitudes with a strain rate of 0.0025%/s ( ${\epsilon }_{\mathrm{a}}$ is the strain amplitude):(a) ${\epsilon }_{\mathrm{a}}=0.05\%$ , (b) ${\epsilon }_{\mathrm{a}}=0.1\%$ , and (c) ${\epsilon }_{\mathrm{a}}=0.15\%$ .

Figure 10.

Stress–strain curves at various strain amplitudes with a strain rate of 0.005%/s ( ${\epsilon }_{\mathrm{a}}$ is the strain amplitude):(a) ${\epsilon }_{\mathrm{a}}=0.05\%$ , (b) ${\epsilon }_{\mathrm{a}}=0.1\%$ , and (c) ${\epsilon }_{\mathrm{a}}=0.15\%$ .

Figure 11.

Stress–strain curves at various strain amplitudes with a strain rate of 0.01%/s ( ${\epsilon }_{\mathrm{a}}$ is the strain amplitude): (a) ${\epsilon }_{\mathrm{a}}=0.05\%$ , (b) ${\epsilon }_{\mathrm{a}}=0.1\%$ , and (c) ${\epsilon }_{\mathrm{a}}=0.15\%$ .

The data presented in Figures 9 to 11 show obvious hysteretic characteristics of the M2052 damping alloy. The hysteresis curves are spindle-shaped with highly significant nonlinear characteristics. At the same strain rate, the damping characteristics increase with increasing strain amplitude and hysteresis becomes more evident. At the same strain amplitude, changes in the damping characteristics with strain rate are not significant but further experimental studies are needed to define a specific law. Fitted data from the generalized fractional-order Maxwell model are consistent with the test data, and the loading and unloading curves are symmetrical about the centerline, which clearly reflects the nonlinear constitutive relation of the damping alloy. The fitted curves were smoother than the test data, and therefore, more suitable for further analyses. Indices used to assess the goodness of fit of the generalized fractional-order Maxwell model are presented in Table 3.

Table 3.

Indices for assessing goodness of fit of generalized fractional-order Maxwell model.

No.	Strain rate (%/s)	Strain (%)	Genetic algebra	Optimal fitness value	$\alpha$	C	w	Mean square error	Variance	Coefficient of determination
1	0.0025	0.05	1500	0.0249	0.0953	102.3557	0.5126	0.4684	173.7645	0.9941
2		0.1	870	0.0653	0.0578	85.9133	0.5	2.651	2157.9	0.9929
3		0.15	1023	0.5098	0.1099	100.6552	0.5	0.9805	1163.9	0.9989
4	0.005	0.05	1756	0.3817	0.0401	84.6854	0.5036	0.5215	100.1304	0.9939
5		0.1	1100	0.4018	0.1315	97.7623	0.5002	0.5286	219.9084	0.9985
6		0.15	1000	0.5150	0.0769	87.1256	0.5	0.9238	564.4336	0.9991
7	0.01	0.05	1562	0.3498	0.0260	84.2119	0.5380	0.7741	105.2740	0.9959
8		0.1	1200	0.3126	0.0665	83.9569	0.5112	1.2688	324.8074	0.9985
9		0.15	1180	0.4059	0.1937	110.9619	0.5001	1.6425	591.2967	0.9990

The generalized fractional-order Maxwell model fit well with the test data, with a mean square error of 0.4684–2.651 and a coefficient of determination greater than 0.9929 (Table 3). Furthermore, to more conveniently determine parameter values of the model under different loading conditions within a certain elastic range and without test data, the coefficient averaging method for fractional-order Maxwell models was adopted. ²³ Average values of the various coefficients were calculated under three different amplitudes based on the obtained fitted coefficients. Then, expressions for the average coefficient values relative to changes in strain amplitude were fitted and the applicability of the generalized fractional-order Maxwell model was extended. Table 4 presents expressions for the various coefficients of the generalized fractional-order Maxwell model relative to strain.

Table 4.

Expressions of coefficients relative to strain.

No.	Name	Expression
1	$\alpha$	$\alpha =0.07421\times e^{(\frac{\epsilon }{0.0018})}-0.04423$
2	$C$	$C=e^{(4.64128-397.08966*\epsilon +246810.18178*{\epsilon }^2)}$
3	$w$	$w=0.07342\times e^{(-\frac{\epsilon }{3.7545\times {10}^{-4}})}+0.49868$

Secondary development of constitutive model for M2052 damping alloy using user subroutine

The self-installed material library in ANSYS contains no description of the M2052 damping alloy, therefore, as a prerequisite, user-programmable features (UPFs) of the software were used to write a secondary development using the custom user material (USERMAT) subroutine and to describe the physical properties, damping characteristics, and constitutive relation of the M2052 damping alloy. Then, the finite element software was used to perform a simulation analysis to investigate the dynamic characteristics of the Mn-Cu damping alloy subjected to an impact load.

The UPFs in ANSYS refer to the self-installed secondary development platform. The platform is based on the syntax rules of Fortran and designed to offer reprogrammable subroutines and functions to users. In addition, users can create material or element properties according to specific applications, define failure criteria, develop optimization algorithms, etc. More specifically, the primary function of a USERMAT subroutine is to define the constitutive relation of a material and to provide the strain increment and to calculate the stress increment using a user-defined main routine. Thus, a new stress can be obtained and the stress update process can be completed. A flowchart of the USERMAT development process is presented in Figure 12.

Figure 12.

Flowchart of USERMAT development process.

Consistent tangent stiffness matrix

The USERMAT subroutine must accomplish two tasks: (1) Calculate the stress increment from a given strain increment, that is, the stress update process; (2) Define the Jacobi matrix of the material, that is, the consistent tangent stiffness matrix. ²⁴

The consistent tangent stiffness matrix can be defined as

$$D_{\mathrm{eq}}=\frac{\partial \mathrm{\Delta }\sigma }{\partial \mathrm{\Delta }\epsilon }$$

(16)

where $\sigma$ is the Cauchy stress tensor; $\mathrm{\Delta }\sigma$ is the stress increment; $\epsilon$ is the Cauchy strain tensor; $\mathrm{\Delta }\epsilon$ is the strain increment.

To solve the consistent tangent stiffness matrix, the main problem of the generalized fractional-order Maxwell model must first be solved, that is, the inconvenience of expressing the finite difference form as the relation between stress increment and strain increment. The solution is to transform the constitutive model and extend the generalized fractional-order Maxwell model from a one-dimensional state to a three-dimensional state.

According to the constitutive equation of the generalized fractional-order Maxwell model, presented in equation (2), we have

$${\mathrm{wD}}_t^{\alpha }\sigma +\frac{\mathrm{Ew}}{\kappa }\sigma =D_t^{\alpha }{\sigma }_1+\frac{E}{\kappa }{\sigma }_1={\mathrm{ED}}_t^{\alpha }\epsilon$$

(17)

When $0<h<1$ , $a=0,n=1$ , then $D_t^h\sigma (t)$ ( $0<h<1$ ) can be written as

$$D_t^h\sigma (t)=\frac{1}{\mathrm{\Gamma }(1-h)}{\int }_0^t{(t-\xi )}^{-h}\sigma \prime (\xi )d\xi$$

(18)

After defining $k_l\left(t\right)=\frac{t^{l-1}}{\mathrm{\Gamma }\left(l\right)},t>0$ , the h^th-order fractional Caputo derivative can be expressed in convolution form as

$$D_t^h\sigma (t)=\sigma \prime (\mathrm{t})*k_{1-h}(t)$$

(19)

The Laplace transform can be applied to both sides of equation (19), as follows (for the convenience, “–” denotes the Laplace transform):

$$s^{\alpha }\overline{{\sigma }_1}-s^{\alpha -1}{\sigma }_1\left(0\right)+\frac{\overline{{\sigma }_1}}{{\tau }^{\alpha }}=\mathrm{E}(s^{\alpha }\overline{\epsilon }-s^{\alpha -1}\epsilon \left(0\right))$$

(20)

${\sigma }_1\left(0\right)=\mathrm{E}\epsilon \left(0\right)$ , therefore

$$s^{\alpha }\overline{{\sigma }_1}+\frac{\overline{{\sigma }_1}}{{\tau }^{\alpha }}=\mathrm{E}{s}^{\alpha }\overline{\epsilon }$$

(21)

Combined with the transformation presented in equation (2), we obtain

$$w\overline{\sigma }=\frac{\mathrm{E}{s}^{\alpha }}{s^{\alpha }+\frac{1}{{\tau }^{\alpha }}}\overline{\epsilon }$$

(22)

Under stress relaxation conditions, the expression of the relaxation modulus $Y\left(\mathrm{t}\right)$ can be obtained by substituting $\overline{\epsilon }=\frac{{\epsilon }_0}{s}$ into equation (22) and applying the inverse Laplace transform, expressed as

$$Y\left(\mathrm{t}\right)=\frac{E}{w}{L}^{-1}\left[\frac{s^{\alpha -1}}{s^{\alpha }+\frac{1}{{\tau }^{\alpha }}}\right]$$

(23)

To perform an inverse Laplace transform on the relaxation modulus, $E_{\alpha ,\beta }(z)$ is introduced as the generalized Mittag-Leffler (M-L) function,^25,26 defined as

$$E_{\alpha ,\beta }(z)=\sum _{i=0}^{\infty }\frac{z^i}{\mathrm{\Gamma }(\alpha i+\beta )};(\alpha ,\beta >0)$$

(24)

According to Wang, ²⁷ the Laplace transform of the two-parameter M-L function is

$$L[t^{\alpha k+\beta -1}{E}_{\alpha ,\beta }^{(k)}(\pm a{t}^{\alpha })]=\frac{k!s^{\alpha -\beta }}{{(s^{\alpha }\mp a)}^{k+1}}$$

(25)

where

$$E_{\alpha ,\beta }^{(k)}(\mathrm{z})=\frac{d^k}{d{z}^k}{E}_{\alpha ,\beta }(\mathrm{z})=\sum _{j=0}^{\infty }\frac{(j+k)!z^j}{j!\mathrm{\Gamma }(\alpha (j+k)+\beta )}$$

(26)

Then, assuming $k=0,\alpha =\alpha ,\beta =1,a=(\frac{1}{\tau }{)}^{\alpha }$ , we have

$$L^{-1}\left[\frac{s^{\alpha -1}}{s^{\alpha }+\frac{1}{{\tau }^{\alpha }}}\right]=E_{\alpha ,1}[-{\left(\frac{t}{\tau }\right)}^{\alpha }]$$

(27)

Substituting Formula (27) into Formula (23) yields

$$Y\left(\mathrm{t}\right)=\frac{\mathrm{E}}{w}{E}_{\alpha ,1}[-{\left(\frac{t}{\tau }\right)}^{\alpha }]$$

(28)

Furthermore, according to the definition of Stieltjes convolution: ²⁸

It can be assumed that function $f\left(t\right)$ is continuous on $0\le t<\infty$ and when $t<0$ , $f\left(t\right)=0$ ; $g\left(t\right)$ is a function defined on $-\infty <t<+\infty$ ; when $t\to -\infty$ , $g\left(t\right)\to 0$ , then

$$f*\mathrm{dg}=f(t)g(0)+{\int }_0^{t}f(t-\tau )\frac{\partial g(\tau )}{\partial \tau }d\tau$$

(29)

$$L[f*\mathrm{dg}]=s\overline{f}(s)\overline{g}(s)$$

(30)

Thus, the Laplace transform of the Stieltjes convolution can be used to define equation (30) and to solve the inverse Laplace transform of equation (22). The definition of Stieltjes convolution, Formula (29), can be used to ultimately obtain the constitutive equation of the relaxation modulus expression of the damping alloy based on the generalized fractional-order Maxwell model. This also solves another problem associated with the original constitutive equation, that is, it is inconvenient to express its finite difference form as the relation between stress increment and strain increment, thereby limiting its application in the secondary development of the finite element program.

$$\begin{array}{c}\multicolumn{1}{c}{\sigma \left(t\right)=Y(t)*d\epsilon }\\ =Y(t){\epsilon }_0+{\int }_0^{t}Y(t-t\prime )\frac{\partial \epsilon (t\prime )}{\partial t\prime }\mathrm{dt}\prime \hfill \\ =\frac{\mathrm{E}}{w}{E}_{\alpha ,1}[-{\left(\frac{t}{\tau }\right)}^{\alpha }]{\epsilon }_0+{\int }_0^t\frac{\mathrm{E}}{w}{E}_{\alpha ,1}[-{\left(\frac{t-t\prime }{\tau }\right)}^{\alpha }]\frac{\partial \epsilon (t\prime )}{\partial t\prime }\mathrm{dt}\prime \hfill \end{array}$$

(31)

The working range of the damping alloy is within the elastic strain range, and only minor deformation occurs. Thus, the Cauchy stress tensor ${\sigma }_{\mathrm{ij}}$ and the Cauchy strain tensor ${\epsilon }_{\mathrm{ij}}$ are adopted to extend the constitutive equation, equation (31), to the three-dimensional state, as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{\mathrm{ij}}=\frac{\mathrm{E}}{w}{E}_{\alpha ,1}[-{\left(\frac{t}{\tau }\right)}^{\alpha }]\left[A\right]{\epsilon }_{\mathrm{ij}}(0)}\\ +{\int }_0^t\frac{\mathrm{E}}{w}{E}_{\alpha ,1}[-{\left(\frac{t-t\prime }{\tau }\right)}^{\alpha }]\left[A\right]\frac{\partial {\epsilon }_{\mathrm{ij}}}{\partial t\prime }\mathrm{dt}\prime \hfill \end{array}$$

(32)

where $\left[A\right]$ is the stiffness matrix, expressed as

$$\left[A\right]=\left(\begin{array}{cccccc}\frac{1-\nu }{(1+\nu )(1-2\nu )}\hfill & \frac{\nu }{(1+\nu )(1-2\nu )}\hfill & \frac{\nu }{(1+\nu )(1-2\nu )}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{\nu }{(1+\nu )(1-2\nu )}}& \frac{1-\nu }{(1+\nu )(1-2\nu )}\hfill & \frac{\nu }{(1+\nu )(1-2\nu )}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{\nu }{(1+\nu )(1-2\nu )}}& \frac{\nu }{(1+\nu )(1-2\nu )}\hfill & \frac{1-\nu }{(1+\nu )(1-2\nu )}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \frac{1}{2(1+\nu )}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \frac{1}{2(1+\nu )}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{1}{2(1+\nu )}\hfill \end{array}\right)$$

(33)

The integral term in Formula (32) is related to the deformation history, and can be expressed in incremental form according to the Cauchy stress tensor ${\sigma }_{\mathrm{ij}}$ and the Cauchy strain tensor ${\epsilon }_{\mathrm{ij}}$ , as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{\mathrm{ij}}^t=\frac{\mathrm{E}}{w}{E}_{\alpha ,1}[-{\left(\frac{t}{\tau }\right)}^{\alpha }]\left[A\right]{\epsilon }_{\mathrm{ij}}(0)+{\int }_0^t\frac{\mathrm{E}}{w}{E}_{\alpha ,1}}\\ \multicolumn{1}{c}{[-{\left(\frac{t-t\prime }{\tau }\right)}^{\alpha }]\left[A\right]\frac{\partial {\epsilon }_{\mathrm{ij}}^t}{\partial t\prime }\mathrm{dt}\prime }\end{array}$$

(34)

$$\begin{array}{c}{\sigma }_{\mathrm{ij}}^{\mathrm{t}+\mathrm{\Delta }t}=\frac{\mathrm{E}}{w}{E}_{\alpha ,1}[-{\left(\frac{t+\mathrm{\Delta }t}{\tau }\right)}^{\alpha }]\left[A\right]{\epsilon }_{\mathrm{ij}}(0)\hfill \\ +{\int }_0^{t+\mathrm{\Delta }t}\frac{\mathrm{E}}{w}{E}_{\alpha ,1}[-{\left(\frac{t+\mathrm{\Delta }t-t\prime }{\tau }\right)}^{\alpha }]\left[A\right]\frac{\partial {\epsilon }_{\mathrm{ij}}^{t+\mathrm{\Delta }t}}{\partial t\prime }\mathrm{dt}\prime \hfill \end{array}$$

(35)

When $\mathrm{\Delta }t$ is very small (that is, within the range of $t\to t+\mathrm{\Delta }t$ ), it can be approximately held that the strain rates at moment t and moment $\mathrm{\Delta }t+t$ are equal, that is, $\frac{\partial {\epsilon }_{\mathrm{ij}}^t}{\partial t\prime }=\frac{\partial {\epsilon }_{\mathrm{ij}}^{t+\mathrm{\Delta }t}}{\partial t\prime }=\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}{\mathrm{\Delta }t}$ . In this case, equations (34) and (35) can be simplified, as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{\mathrm{ij}}^t=\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{[-{(\frac{t}{\tau })}^{\alpha }]}^i}{\mathrm{\Gamma }(\alpha i+1)}\left[A\right]{\epsilon }_{\mathrm{ij}}(0)}\\ +\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\left[A\right]\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+2)}{(-1)}^i\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}{\mathrm{\Delta }t}{t}^{\alpha i+1}\hfill \\ \multicolumn{1}{c}{=\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+1)}{(-1)}^i\left[A\right]}\\ \{{\epsilon }_{\mathrm{ij}}(0)t^{\alpha i}+\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}{(\alpha i+1)\mathrm{\Delta }t}{t}^{\alpha i+1}\}\hfill \end{array}$$

(36)

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{\mathrm{ij}}^{t+\mathrm{\Delta }t}=\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{[-{(\frac{t+\mathrm{\Delta }t}{\tau })}^{\alpha }]}^i}{\mathrm{\Gamma }(\alpha i+1)}\left[A\right]{\epsilon }_{\mathrm{ij}}(0)}\\ +\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\left[A\right]\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+2)}{(-1)}^i\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}{\mathrm{\Delta }t}{(t+\mathrm{\Delta }t)}^{\alpha i+1}\hfill \\ =\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+1)}{(-1)}^i\left[A\right]{(1+\frac{\mathrm{\Delta }t}{t})}^{\alpha i}\hfill \\ \{{\epsilon }_{\mathrm{ij}}(0)t^{\alpha i}+\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}{(\alpha i+1)\mathrm{\Delta }t}{t}^{\alpha i+1}\}\hfill \\ +\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+2)}{(-1)}^i\left[A\right]{(t+\mathrm{\Delta }t)}^{\alpha i}\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}\hfill \end{array}$$

(37)

When $\mathrm{\Delta }t$ →0, then $\frac{\mathrm{\Delta }t}{t}$ →0 and the ${(1+\frac{\mathrm{\Delta }t}{t})}^{\alpha i}$ in the above formula can be approximated to 1. In this case, Formula (37) becomes

$$\begin{array}{c}{\sigma }_{\mathrm{ij}}^{t+\mathrm{\Delta }t}=\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+1)}{(-1)}^i\left[A\right]{(1+\frac{\mathrm{\Delta }t}{t})}^{\alpha i}\hfill \\ \{{\epsilon }_{\mathrm{ij}}(0)t^{\alpha i}+\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}{(\alpha i+1)\mathrm{\Delta }t}{t}^{\alpha i+1}\}\hfill \\ +\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+2)}{(-1)}^i\left[A\right]{(t+\mathrm{\Delta }t)}^{\alpha i}\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}\hfill \\ ={\sigma }_{\mathrm{ij}}^t+\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+2)}{(-1)}^i\left[A\right]{(t+\mathrm{\Delta }t)}^{\alpha i}\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}\hfill \end{array}$$

(38)

Then

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Delta }{\sigma }_{\mathrm{ij}}={\sigma }_{\mathrm{ij}}^{t+\mathrm{\Delta }t}-{\sigma }_{\mathrm{ij}}^t}\\ =\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{\left[A\right]{\epsilon }_{\mathrm{ij}}(0)}{\mathrm{\Gamma }(\alpha i+1)}\{{[-{(\frac{t+\mathrm{\Delta }t}{\tau })}^{\alpha }]}^i-{[-{(\frac{t}{\tau })}^{\alpha }]}^i\}\hfill \\ +\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\left[A\right]\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+2)}{(-1)}^i\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}{\mathrm{\Delta }t}[{(t+\mathrm{\Delta }t)}^{\alpha i+1}-t^{\alpha i+1}]\hfill \\ =\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+2)}{(-1)}^i\left[A\right]{(t+\mathrm{\Delta }t)}^{\alpha i}\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}\hfill \end{array}$$

(39)

Thus, the incremental recursion relation of stress can be obtained as

$${\left({\sigma }_{\mathrm{ij}}\right)}_{n+1}={\left({\sigma }_{\mathrm{ij}}\right)}_n+{\left(\mathrm{\Delta }{\sigma }_{\mathrm{ij}}\right)}_n$$

(40)

The consistent tangent stiffness matrix is

$$D_{\mathrm{eq}}^{\mathrm{ij}}=\frac{\partial \mathrm{\Delta }{\sigma }_{\mathrm{ij}}}{\partial \mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}=\sum _{i=0}^{\infty }\frac{\mathrm{E}}{w}\frac{{\tau }^{-\alpha i}}{\mathrm{\Gamma }(\alpha i+2)}{(-1)}^i\left[A\right]{(t+\mathrm{\Delta }t)}^{\alpha i}$$

(41)

Formula (39) is the three-dimensional incremental form of the constitutive equation required to develop the user-defined subroutine; Therefore, equation (40) can be used to calculate the corresponding Cauchy force. The user-defined subroutine development task is used to perform the stress calculations presented in equations (39), (40), and (41) on elements of the finite element model of the material based on the parameters imported from the main routine. Then, the calculated stresses are returned to the main routine.

Verification of user subroutine

To guarantee the correctness of the subroutine, a finite element model of the uniaxial cyclic tensile test was established and the USERMAT subroutine was used to simulate the test conditions in section 2.2. The simulated 150 mm × 20 mm × 2 mm test piece was comprised of 1520 elements and 9053 nodes, as shown in Figure 13.

Figure 13.

Finite element model of test piece used in uniaxial cyclic tensile tests.

A reference strain rate of 0.005%/s and a reference strain amplitude of 0.1% were adopted. One end of the tensile sample was defined as the fixed constraint and a constant-speed displacement was applied to the other end to control the load. The boundary conditions are illustrated in Figure 14.

Figure 14.

Model constraints.

A comparison of the simulation results and test results is shown in Figure 15.

Figure 15.

Comparison of stresses and strains obtained from the uniaxial cyclic tensile test and USERMAT simulation.

As shown in Figure 15, the simulation results are consistent with the real test results, and errors were within 5%, suggesting the subroutine accurately describes the constitutive properties of the M2052 damping alloy.

Dynamic characteristics of various proposed applications of damping alloy cantilever beams subjected to impact load

The damping alloy can be applied in industrial vibration reduction in two ways: first, replacing the overall structure, which is convenient in terms of design, but consumes more material and changes the overall stiffness and damping characteristics of the original component; attaching a damping layer to the surface of the original structure, which consumes less material and can achieve certain vibration reduction effects without affecting the overall stiffness of the original structure, however, the specific application and dosage of the damping alloy must be considered. After the secondary development of the model, the constitutive model of the Mn-Cu damping alloy was used to conveniently compare different applications and alloy consumption rates.

Eight application schemes were proposed based on finite element models for the cantilever beams established earlier in the paper (Figure 16). Details of each application are listed in Table 5.

Figure 16.

Finite element models of cantilever beams: (a) single-layer beam, (b) double-layer beam, and (c) sandwich beam.

Table 5.

Applications of damping alloy.

No.	Form of cantilever beam	Thickness of damping alloy layer	Thickness of steel layer
1	Single-layer beam	0 mm	5 mm
2	Double-layer beam	1 mm	4 mm
3	Double-layer beam	2 mm	3 mm
4	Double-layer beam	3 mm	4 mm
5	Double-layer beam	4 mm	1 mm
6	Single-layer beam	5 mm	0 mm
7	Sandwich beam	1 mm at both the top and the bottom	3 mm
8	Sandwich beam	2 mm at both the top and the bottom	1 mm

A fixed constraint was imposed on one end of the cantilever beam and an impact load of 100 N was applied to the other end for a duration of 0.002 s. The element type was a solid186 element and the simulation duration was 1 s. The nonlinear Newton-Simpson algorithm was adopted with Rayleigh damping for steel. The damping coefficient of steel was adopted. Table 6 presents the specific material parameters used in the simulations.

Table 6.

Material parameters.

Name	Young’s modulus	Poisson’s ratio	Yield strength	Density	Mass/damping coefficient alpha	Stiffness/damping coefficient beta
M2052	68.5 GPa	0.338	205 MPa	7.25 g/cm3	–	–
Steel	201 GPa	0.3	250 MPa	7.8 g/cm3	0.3113	8.27 × 10−6

The vibration displacement was monitored at the free end of the cantilever beam. The nonlinear transient simulation results for the single-layer beam are shown in Figures 17 and 18. Figure 17 compares the test and simulation results for a single-layer Mn-Cu damping alloy beam.

Figure 17.

Comparison of simulation and test results for single-layer Mn-Cu damping alloy beam.

Figure 18.

Comparison of vibration displacements of single-layer beams subjected to impact load.

Results obtained with the USERMAT subroutine of the constitutive model of the Mn-Cu damping alloy agree well with the test values. The simulated maximum vibration displacement was 1.252 mm and the corresponding test value was 1.352 mm, with a relative error of 7%. Consistency between the vibration attenuation time also demonstrated the correctness of the simulation results. Figure 18 compares the vibration displacements obtained for a single-layer Mn-Cu damping alloy beam and an all-steel beam.

As seen in Figure 18, the all-steel beam had a vibration displacement amplitude of 2.016 mm, whereas the M2052 damping alloy cantilever beam had a maximum amplitude of 1.352 mm under the same impact load, that only 67.06% of all-steel beam. The vibration displacement of the damping alloy stabilized at a significantly faster pace compared to structural steel and almost became stable at around 0.4 s, however, the all-steel beam could not restore the initial state after 1 s. Figure 19 compares the vibration displacements of different combinations of double-layer Mn-Cu damping alloy beams subjected to an impact load.

Figure 19.

Comparison of vibration displacements of double-layer beams subjected to impact load.

Figure 19 shows that the vibration displacements of the four damping alloy composite beams are different. The maximum amplitudes were all within the range of 1.51–1.65 mm, however, the rate of decrease in the vibration amplitude varied and could be ranked in descending order, as follows: 4 mm composite beam >2 mm composite beam >3 mm composite beam >1 mm composite beam. Clearly, the vibration reduction did not increase linearly with increasing damping alloy content. For example, the 3 mm composite beam had a higher damping alloy content, however, the vibration reduction effect was closer to that of a 1mm composite beam. This was likely because the amplitude is affected by matching of the matrix steel plate and damping layer content under the same excitation. A higher damping alloy content led to lower stiffness and larger deformation, whereas a lower damping alloy content resulted in higher stiffness and less deformation. The amount of deformation was positively correlated with the damping characteristics of the damping alloy. Comparing results for the double-layer beams, it can be seen that the 4 mm damping alloy double-layer beam performed best; however, from a cost efficiency and material processing complexity perspective, the double-layer beam with 2 mm layer of damping alloy had the highest cost performance and its vibration reduction effect was close to that of the double-layer beam with 4 mm layer of damping alloy. Figure 20 compares the results of the sandwich beams.

Figure 20.

Comparison of vibration displacement of sandwich beams subjected to impact load.

A comparison of the vibration displacements of the two sandwich beams subjected to an impact load (Figure 20), shows that the M2052 (2 mm) + steel (1 mm) + M2052 (2 mm) sandwich beam exhibits satisfactory damping effects with a maximum vibration displacement of 1.226 mm, which is better than the result for the (1.283 mm) of M2052 (1 mm) + steel (3 mm) + M2052 (1 mm) sandwich beam. The sandwich beams had a higher rate of vibration attenuation and regained stability within about 0.25 s.

To better analyze the vibration displacements of different beam combinations subjected to an impact load, vibration displacement amplitude curves were plotted, as illustrated in Figure 21. Figure 7 shows the vibration response of the cantilever beams in various groups.

Figure 21.

Vibration displacement amplitudes of cantilever beams subjected to impact load.

According to Table 7, for the different beam combinations and consumption rates of the damping alloys, the percentage decrease in the maximum amplitude decline percentage and curvature will change as well. The rate of decrease in the maximum amplitude for the eight application schemes can be ranked in descending order, as follows: 7 > 8 > 6 > 2 > 5 > 3 > 4 > 1. Curvature represents the degree of bending of a curve, and in relation to the angle of vibration amplitude attenuation, represents the rate of attenuation: 7 > 6 > 5 > 3 > 8 > 4 > 2 > 1. In summary, seen from the application forms of the damping alloy, sandwich beams have a larger effect on the vibration attenuation rate than double-layer and single-layer beams. The rates of decrease in the maximum vibration amplitude of the double-layer beams were consistently lower than those of pure damping alloy beams. However, the attenuation rates of M2052 (4 mm) + steel (1 mm) and M2052 (2 mm) + steel (3 mm) were higher than those of pure damping alloy beams, and their vibration reduction effects were superior.

Table 7.

Vibration responses of cantilever beams in various groups.

No.	Type	Maximum amplitude (mm)	Decline percentage of maximum amplitude (%)	Maximum curvature
1	Steel (5 mm)	2.016	0	0.04346
2	M2052 (1 mm) + steel (4 mm)	1.454	27.88	0.0646
3	M2052 (2 mm) + steel (3 mm)	1.626	19.35	0.08292
4	M2052 (3 mm) + steel (2 mm)	1.654	17.96	0.06734
5	M2052 (4 mm) + steel (1 mm)	1.497	25.74	0.0882
6	M2052 (1 mm) + steel (3 mm) + M2052 (1 mm)	1.229	39.04	0.15099
7	M2052 (2 mm) + steel (1 mm) + M2052 (2 mm)	1.283	36.36	0.16165
8	M2052 (5 mm)	1.266	37.20	0.07582

In combination with the above results for double-layer beams and sandwich beams, four superior schemes were identified: M2052 (4 mm) + steel (1 mm), M2052 (2 mm) + steel (3 mm), M2052 (2 mm) + steel (1 mm) +M2052 (2 mm), and M2052 (1 mm) + steel (3 mm) +M2052 (1 mm). However, in terms of both cost efficiency and vibration reduction effect, the preferred schemes are M2052 (1 mm) + steel (3 mm) + M2052 (1 mm) sandwich beam and M2052 (2 mm) + steel (3 mm) double-layer beam. Although the other two schemes, M2052 (2 mm) + steel (1 mm) +M2052 (2 mm) and M2052 (4 mm) + steel (1 mm), achieved satisfactory vibration reduction effects, the higher damping alloy content is more expensive and causes a substantial decline in the overall stiffness of the cantilever beam.

In order to further verify the simulation results, M2052 (2 mm) + Steel (3 mm) and M2052 (1 mm) + Steel (3 mm) + M2052 (1 mm) composite beams were processed according to the simulation dimensions. The materials were fixed by spot welding, and the impact load vibration test was conducted according to the simulation loading conditions. The physical drawings of the two composite beams are shown in Figure 22.

Figure 22.

Photograph of M2052 (2 mm) + steel (3 mm) and M2052 (1 mm) + steel (3 mm) + M2052 (1 mm) composite beams.

Figure 23 shows the comparison between the test and simulation results of M2052 (2 mm) + Steel (3 mm) double-layer beam.

Figure 23.

Comparison of simulation and test results for M2052 (2 mm) + steel (3 mm) double-layer beam.

As can be seen from Figure 23, the simulation curve is basically consistent with the test curve. The maximum vibration displacement of M2052 (2 mm) + Steel (3 mm) double-beam subjected to 100 N impact load is 1.476 mm, while the simulation is 1.626 mm, with an error value of 10.16%. However, the experimental amplitude attenuation time of M2052 (2 mm) + Steel (3 mm) double-beam is 0.4098 s, while the simulation is 0.3733 s, with an error value of 8.91%. Figure 24 shows the comparison between the test and simulation results of M2052 (1 mm) + Steel (3 mm) + M2052 (1 mm) sandwich beam.

Figure 24.

Comparison of simulation and test results for M2052 (1 mm) + steel (3 mm) + M2052 (1 mm) sandwich beam.

As shown in Figure 24, the maximum vibration displacement of M2052 (1 mm) + steel (3 mm) + M2052 (1 mm) sandwich beam subjected to impact load is 1.253 mm, while the simulation is 1.283 mm, with an error value of 2.39%. Although the vibration of test and simulation all have declined to 0 in 0.32 s, the attenuation rate of test is faster than the simulation at early stage and slower than the simulation at later stage. Due to the influence of the air resistance in test environment, the resistance leads to larger amplitude decay at the early stage. However, as the velocity decreases in the later period, the effect of air resistance also decreases, which results in a slower amplitude attenuation speed. But the simulation doesn’t set the air resistance in the environment. Overall, the simulation results are in agreement with the test results, which verifies the feasibility of simulating the vibration of damping alloy under impact load by the secondary development program of constitutive relation of Mn-Cu damping alloy.

Conclusions

In this study, vibration tests were conducted on Mn-Cu damping alloy cantilever beams subjected to an impact load. Deviations of existing finite element simulation methods for damping alloys from reality were first analyzed, the generalized fractional-order Maxwell model for Mn-Cu damping alloy was established. Next, data from uniaxial cyclic tensile tests at a constant strain rate were used to perform parameter fitting and to verify the proposed model. The generalized fractional-order Maxwell model can satisfactorily reflect the nonlinear constitutive relation of the damping alloy. The constitutive model to the three-dimensional state, a secondary development was performed, and dynamic responses of eight different forms of damping alloy cantilever beams were carried out. Finally, it can be concluded that the preferred vibration reduction scheme for cantilever beam structures should be an attached damping layer on both upper and lower sides of the beam (sandwich beam) or a damping layer on one side. Among the eight schemes, a sandwich beam with a 1mm layer of damping alloy and or double-layer beam with a 2 mm layer of damping alloy beam achieved the best effects in terms of both cost efficiency and damping effects. The results offer a theoretical reference for the future design of vibration reduction-oriented applications using damping alloys.