# Analytical solution for nonlinear vertical vibration model of mill roll system based on improved complexification averaging method

## Yong Zhu1, Wan-lu Jiang2, Xiang-dong Kong3, Hao-nan Wang4

1, 2, 3, 4Hebei Provincial Key Laboratory of Heavy Machinery Fluid Power Transmission and Control, Yanshan University, Qinhuangdao 066004, Hebei, China

1, 2, 3, 4Key Laboratory of Advanced Forging and Stamping Technology and Science (Yanshan University), Ministry of Education of China, Qinhuangdao 066004, Hebei, China

2Corresponding author

Journal of Vibroengineering, Vol. 18, Issue 8, 2016, p. 5521-5536. https://doi.org/10.21595/jve.2016.17194
Received 21 May 2016; received in revised form 14 August 2016; accepted 17 August 2016; published 31 December 2016

Copyright © 2016 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Views 113
Abstract.

Rolling mill is the core equipment in modern iron and steel industry, and the reliability of its mill roll system (MRS) is the key to ensure the rolling process with high precision, high speed, continuity and stability. However, the MRS possesses some features such as high nonlinearity, time variability and strong coupling. The vertical vibration easily happens in its working process. Nevertheless, the mathematical model of MRS is difficult to be established and hard to be solved. In this paper, the nonlinear dynamics theory and modern signal processing method were introduced to solve this difficult problem. A two degree of freedom (DOF) nonlinear vertical vibration model of MRS was established. And the model was analytically solved by using complexification averaging (CA) method. The solution error of CA method was thoroughly analyzed. Moreover, the CA method was improved by combining the fast empirical mode decomposition (FEMD) method. Research results indicate that the improved CA method presents a significant advantage in improving the solution precision, and can be used to solve strong nonlinear vibration system (NVS) with two DOF.

Keywords: mill roll system, nonlinear vibration, complexification averaging, fast empirical mode decomposition.

#### 1. Introduction

Rolling is a metal forming process in which a metal strip is laminated through a pair of rolls. In modern iron and steel industry, rolling mill is the core equipment in rolling process. And the operational reliability of MRS is the key to guarantee the rolling process with high precision, high speed, continuity and stability . In steady rolling process, the force balance relation will be destroyed when MRS is disturbed by interference force, and the vertical vibration will be produced [2, 3]. The vertical vibration of MRS has a significant effect on rolling precision and product quality. So, the vertical vibration of MRS is worthy of attention. In addition, MRS is a complex dynamics system composed of many parts. It possesses some special characteristics such as high nonlinearity, time-variability and strong coupling. Hence, it is difficult to precisely establish a mathematical model . At present, the nonlinear factors are not taken into consideration in most of modeling, or linearized by the small deviation theory. The linearized model is frequently used to replace the actual nonlinear system [5, 6]. In sometimes, however, it will cause a large error if some small nonlinear factors are ignored in the actual calculation and analysis. With the development of nonlinear science, researchers pay high attention to nonlinear factors, and begin to explore in this field [7, 8]. Therefore, it is worth further research about introducing the nonlinear dynamics theory and method for the modeling of MRS.

The traditional quantitative analysis methods for nonlinear system mainly include Lindstedt-Poincaré method , multi-scale method , average method  and KBM method . The common point of these methods is to treat the nonlinear factors as a small perturbation attached to linear system . These traditional analysis methods require that the nonlinear terms must be very small. They only apply to solve the weak NVS [14, 15]. In 1999, the complex-variable quantization thought was applied to nonlinear coupling vibration by Manevitch L. I. . And a method named complexification averaging (CA) was proposed, which was gradually evolved into a kind of analytical method for NVS. The CA method is not limited to weak NVS. It not only can be utilized to solve the system with strong nonlinear terms, but also can be used to solve the system with multiple degree of freedom (MDOF) . However, the CA method still has some errors when it is used to solve the strong NVS with MDOF. Therefore, the solution accuracy still needs to be further improved.

In recent years, an adaptive signal processing method named empirical mode decomposition (EMD) was proposed by professor Huang etc. . The method changed the helpless predicament for dealing with nonlinear and non-stationary signal. It can decompose the non-stationary signal into a set of steady data series. So, it has aroused great attention in the international academic circles. Since EMD was put forward, it has been applied to many researches such as physics , biomedicine , machinery fault diagnosis , image analysis , pattern recognition  and so on. A series of valuable research results were obtained. However, in the process of researches, it was found that although EMD was very suitable for processing non-stationary signal, the mode mixing phenomenon exists . In order to overcome the defects of EMD, a novel method called fast empirical mode decomposition (FEMD) was put forward by professor Wang etc. . And it was proved that FEMD could significantly improve computation speed and effectively restrain mode mixing.

In this work, we attempt to utilize the nonlinear dynamics theory and modern signal processing method to explore the analytical solution for nonlinear vertical vibration model of MRS. In Section 2, the nonlinear vertical vibration model of MRS is established. In Section 3, the model is analytically solved by CA method. In Section 4, in order to improve the solution accuracy, an improved CA method is proposed by combining FEMD method, so that the improved CA method can be used to solve strong NVS. Finally, some conclusions are provided in Section 5.

#### 2. Nonlinear vertical vibration model of mill roll system

MRS is a distributed mass system with MDOF. So, the analysis for the system are very complicated. At present, in order to facilitate the analysis, MRS is mainly divided into two load models: the lumped load model with single DOF and the distributed load model with MDOF. Moreover, researches show that the performance of the asymmetric distribution model with two DOF more correspond with the actual working condition [26, 27].

In this paper, the four rollers mill is employed as a research object, and the distributed mass model with two DOF is adopted as foundation. Then, the upper backup roll and work roll are seen as a mass system, and the lower backup roll and work roll are served as another mass system. Moreover, the nonlinear elastic force and nonlinear damping force are considered. The structural form of Duffing oscillator and Van der pol oscillator are respectively defined as the nonlinear stiffness term and nonlinear damping term [28, 29]. Furthermore, a nonlinear vertical vibration model of MRS is established, as shown in Fig. 1.

The dynamic equilibrium equation of MRS at vertical vibration can be expressed by:

(1)
${m}_{1}\stackrel{¨}{y}+{c}_{1}\stackrel{˙}{y}+{c}_{3}{y}^{2}\stackrel{˙}{y}+{k}_{1}y+{k}_{3}{y}^{3}+{c}_{12}\left(\stackrel{˙}{y}-\stackrel{˙}{z}\right)+{k}_{12}\left(y-z\right)=0,$
(2)
${m}_{2}\stackrel{¨}{z}+{c}_{2}\stackrel{˙}{z}+{c}_{4}{z}^{2}\stackrel{˙}{z}+{k}_{2}z+{k}_{4}{z}^{3}-{c}_{12}\left(\stackrel{˙}{y}-\stackrel{˙}{z}\right)-{k}_{12}\left(y-z\right)=0,$

where ${m}_{1}$ and ${m}_{2}$ are the equivalent quality of moving parts of the upper roll system (URS) and lower roll system (LRS); ${c}_{1}$ and ${c}_{2}$ are the linear damping coefficient (DC) of moving parts of the URS and LRS; ${c}_{3}$ and ${c}_{4}$ are the nonlinear DC; ${k}_{1}$ and ${k}_{2}$ are the linear stiffness coefficient (SC) between the frame beam and the moving parts of URS or LRS; ${k}_{3}$ and ${k}_{4}$ are the nonlinear SC; ${c}_{12}$ is the DC between work roll and rolled piece; ${k}_{12}$ is the SC between work roll and rolled piece; $y$ and $z$ are the vibration displacement of the URS and LRS.

Fig. 1. Structure of four rollers mill and its nonlinear mechanics model  #### 3.1. Basic principle of CA method

In 1999, a pair of conjugate complex variables and some slow-varying parameters were introduced to solve vibration equation by Manevitch L. I. Therefrom, the CA method was gradually evolved into a kind of analytical method for NVS . The basic steps of CA method are as follows :

(1) Complexification of the nonlinear vibration equation;

(2) Partition the vibration equation into fast and slow components;

(3) Eliminate the fast-varying components by removing secular terms;

(4) Extract the slow-varying components as the approximate solution of system equation.

#### 3.2. Solve the nonlinear model by using CA method

(1) Complexification of the nonlinear vibration equation.

Suppose the vertical vibration responses of MRS are:

(3)

Four complex variables are introduced as:

(4)

where ${\omega }_{1}$ and ${\omega }_{2}$ are the natural angular frequencies of system response.

Then, a complex change is performed to the vibration responses $y\left(t\right)$ and $z\left(t\right)$, there are:

(5)

$\stackrel{¨}{y}={\stackrel{˙}{\psi }}_{1}-i{\omega }_{1}\frac{{\psi }_{1}+{\psi }_{1}^{\mathrm{*}}}{2}+{\stackrel{˙}{\psi }}_{2}-i{\omega }_{2}\frac{{\psi }_{2}+{\psi }_{2}^{\mathrm{*}}}{2},$
(6)

$\stackrel{¨}{z}={\stackrel{˙}{\psi }}_{3}-i{\omega }_{1}\frac{{\psi }_{3}+{\psi }_{3}^{\mathrm{*}}}{2}+{\stackrel{˙}{\psi }}_{4}-i{\omega }_{2}\frac{{\psi }_{4}+{\psi }_{4}^{\mathrm{*}}}{2},$

where $*$ denotes conjugation.

The Eqs. (5) and (6) are substituted into Eqs. (1) and (2), there are:

(7)
${m}_{1}\left({\stackrel{˙}{\psi }}_{1}-i{\omega }_{1}\frac{{\psi }_{1}+{\psi }_{1}^{\mathrm{*}}}{2}+{\stackrel{˙}{\psi }}_{2}-i{\omega }_{2}\frac{{\psi }_{2}+{\psi }_{2}^{\mathrm{*}}}{2}\right)+{c}_{1}\left(\frac{{\psi }_{1}+{\psi }_{1}^{\mathrm{*}}}{2}+\frac{{\psi }_{2}+{\psi }_{2}^{\mathrm{*}}}{2}\right)$

(8)
${m}_{2}\left({\stackrel{˙}{\psi }}_{3}-i{\omega }_{1}\frac{{\psi }_{3}+{\psi }_{3}^{\mathrm{*}}}{2}+{\stackrel{˙}{\psi }}_{4}-i{\omega }_{2}\frac{{\psi }_{4}+{\psi }_{4}^{\mathrm{*}}}{2}\right)+{c}_{2}\left(\frac{{\psi }_{3}+{\psi }_{3}^{\mathrm{*}}}{2}+\frac{{\psi }_{4}+{\psi }_{4}^{\mathrm{*}}}{2}\right)$

(2) Partition the vibration equation into fast and slow components.

The complex variable $\psi \left(t\right)$ is decomposed into two parts: $\phi \left(t\right)$ represents the slow-varying component and ${e}^{i\omega t}$ represents the fast-varying component:

(9)

The Eq. (9) is substituted into Eqs. (7) and (8), after being further trimmed, there are:

(10)
${m}_{1}\left[\left({\stackrel{˙}{\phi }}_{1}+\frac{i{\omega }_{1}{\phi }_{1}}{2}\right){e}^{i{\omega }_{1}t}-\frac{i{\omega }_{1}{\phi }_{1}^{\mathrm{*}}}{2}{e}^{-i{\omega }_{1}t}+\left({\stackrel{˙}{\phi }}_{2}+\frac{i{\omega }_{2}{\phi }_{2}}{2}\right){e}^{i{\omega }_{2}t}-\frac{i{\omega }_{2}{\phi }_{2}^{\mathrm{*}}}{2}{e}^{-i{\omega }_{2}t}\right]$

(11)
${m}_{2}\left[\left({\stackrel{˙}{\phi }}_{3}+\frac{i{\omega }_{1}{\phi }_{3}}{2}\right){e}^{i{\omega }_{1}t}+\frac{i{\omega }_{1}{\phi }_{3}^{\mathrm{*}}}{2}{e}^{-i{\omega }_{1}t}+\left({\stackrel{˙}{\phi }}_{4}+\frac{i{\omega }_{2}{\phi }_{4}}{2}\right){e}^{i{\omega }_{2}t}+\frac{i{\omega }_{2}{\phi }_{4}^{\mathrm{*}}}{2}{e}^{-i{\omega }_{2}t}\right]$

(3) Eliminate the fast-varying components by removing secular terms.

The fast-varying components are eliminated by removing the secular terms ${e}^{i{\omega }_{1}t}$ and ${e}^{i{\omega }_{2}t}$. The coefficients of ${e}^{i{\omega }_{1}t}$ and ${e}^{i{\omega }_{2}t}$ derived from Eqs. (10) and (11) are set to zero, there are:

(12)
${m}_{1}\left({\stackrel{˙}{\phi }}_{1}+\frac{i{\omega }_{1}{\phi }_{1}}{2}\right)+{c}_{1}\frac{{\phi }_{1}}{2}+{c}_{3}\left(\frac{{\left|{\phi }_{1}\right|}^{2}{\phi }_{1}}{8{\omega }_{1}^{2}}+\frac{{\left|{\phi }_{2}\right|}^{2}{\phi }_{1}}{4{\omega }_{2}^{2}}\right)$
(13)
${m}_{1}\left({\stackrel{˙}{\phi }}_{2}+\frac{i{\omega }_{2}{\phi }_{2}}{2}\right)+{c}_{1}\frac{{\phi }_{2}}{2}+{c}_{3}\left(\frac{{\left|{\phi }_{2}\right|}^{2}{\phi }_{2}}{8{\omega }_{2}^{2}}+\frac{{\left|{\phi }_{1}\right|}^{2}{\phi }_{2}}{4{\omega }_{1}^{2}}\right)$
(14)
${m}_{2}\left({\stackrel{˙}{\phi }}_{3}+\frac{i{\omega }_{1}{\phi }_{3}}{2}\right)+{c}_{2}\frac{{\phi }_{3}}{2}+{c}_{4}\left(\frac{{\left|{\phi }_{3}\right|}^{2}{\phi }_{3}}{8{\omega }_{1}^{2}}+\frac{{\left|{\phi }_{4}\right|}^{2}{\phi }_{3}}{4{\omega }_{2}^{2}}\right)$
(15)
${m}_{2}\left({\stackrel{˙}{\phi }}_{4}+\frac{i{\omega }_{2}{\phi }_{4}}{2}\right)+{c}_{2}\frac{{\phi }_{4}}{2}+{c}_{4}\left(\frac{{\left|{\phi }_{4}\right|}^{2}{\phi }_{4}}{8{\omega }_{2}^{2}}+\frac{{\left|{\phi }_{3}\right|}^{2}{\phi }_{4}}{4{\omega }_{1}^{2}}\right)$

(4) Extract the slow-varying components as the approximate solution of system equation.

The slow-varying components $\phi \left(t\right)$ of the nonlinear equations are expressed in the form of amplitude $a$ and phase angle $\theta$:

(16)

Then, the Eq. (16) is respectively substituted into Eqs. (12)-(15), there are:

(17)
${m}_{1}\left({\stackrel{˙}{a}}_{1}{e}^{i{\theta }_{1}}+{a}_{1}{\stackrel{˙}{\theta }}_{1}{e}^{i{\theta }_{1}}i+\frac{{\omega }_{1}{a}_{1}{e}^{i{\theta }_{1}}}{2}i\right)+{c}_{1}\frac{{a}_{1}{e}^{i{\theta }_{1}}}{2}+{c}_{3}\left(\frac{{a}_{1}^{3}{e}^{i{\theta }_{1}}}{8{\omega }_{1}^{2}}+\frac{{a}_{2}^{2}{a}_{1}{e}^{i{\theta }_{1}}}{4{\omega }_{2}^{2}}\right)$

(18)
${m}_{1}\left({\stackrel{˙}{a}}_{2}{e}^{i{\theta }_{2}}+{a}_{2}{\stackrel{˙}{\theta }}_{2}{e}^{i{\theta }_{2}}i+\frac{{\omega }_{2}{a}_{2}{e}^{i{\theta }_{2}}}{2}i\right)+{c}_{1}\frac{{a}_{2}{e}^{i{\theta }_{2}}}{2}+{c}_{3}\left(\frac{{a}_{2}^{3}{e}^{i{\theta }_{2}}}{8{\omega }_{2}^{2}}+\frac{{a}_{1}^{2}{a}_{2}{e}^{i{\theta }_{2}}}{4{\omega }_{1}^{2}}\right)$

(19)
${m}_{2}\left({\stackrel{˙}{a}}_{3}{e}^{i{\theta }_{3}}+{a}_{3}{\stackrel{˙}{\theta }}_{3}{e}^{i{\theta }_{3}}i+\frac{{\omega }_{1}{a}_{3}{e}^{i{\theta }_{3}}}{2}i\right)+{c}_{2}\frac{{a}_{3}{e}^{i{\theta }_{3}}}{2}+{c}_{4}\left(\frac{{a}_{3}^{3}{e}^{i{\theta }_{3}}}{8{\omega }_{1}^{2}}+\frac{{a}_{4}^{2}{a}_{3}{e}^{i{\theta }_{3}}}{4{\omega }_{2}^{2}}\right)$

(20)
${m}_{2}\left({\stackrel{˙}{a}}_{4}{e}^{i{\theta }_{4}}+{a}_{4}{\stackrel{˙}{\theta }}_{4}{e}^{i{\theta }_{4}}i+\frac{{\omega }_{2}{a}_{4}{e}^{i{\theta }_{4}}}{2}i\right)+{c}_{2}\frac{{a}_{4}{e}^{i{\theta }_{4}}}{2}+{c}_{4}\left(\frac{{a}_{4}^{3}{e}^{i{\theta }_{4}}}{8{\omega }_{2}^{2}}+\frac{{a}_{3}^{2}{a}_{4}{e}^{i{\theta }_{4}}}{4{\omega }_{1}^{2}}\right)$

The real parts and imaginary parts of Eqs. (17)-(20) are respectively set to zero. The mathematical transformation ${e}^{i\theta }=\mathrm{c}\mathrm{o}\mathrm{s}\theta +i\mathrm{s}\mathrm{i}\mathrm{n}\theta$ is also substituted. Then, the slow-varying model of the original system can be obtained as:

(21)
${m}_{1}{\stackrel{˙}{a}}_{1}+{c}_{1}\frac{{a}_{1}}{2}+{c}_{3}\left(\frac{{a}_{1}^{3}}{8{\omega }_{1}^{2}}+\frac{{a}_{2}^{2}{a}_{1}}{4{\omega }_{2}^{2}}\right)+{c}_{12}\frac{{a}_{1}-{a}_{3}\mathrm{cos}\left({\theta }_{3}-{\theta }_{1}\right)}{2}-{k}_{12}\frac{{a}_{3}\mathrm{sin}\left({\theta }_{3}-{\theta }_{1}\right)}{2{\omega }_{1}}=0,$
(22)
${m}_{1}\left(2{\omega }_{1}{\stackrel{˙}{\theta }}_{1}+{\omega }_{1}^{2}\right)-{k}_{1}-{k}_{3}\left(\frac{3{a}_{1}^{2}}{4{\omega }_{1}^{2}}+\frac{3{a}_{2}^{2}}{2{\omega }_{2}^{2}}\right)-{c}_{12}\frac{{a}_{3}{\omega }_{1}\mathrm{sin}\left({\theta }_{3}-{\theta }_{1}\right)}{{a}_{1}}$
(23)
${m}_{1}{\stackrel{˙}{a}}_{2}+{c}_{1}\frac{{a}_{2}}{2}+{c}_{3}\left(\frac{{a}_{2}^{3}}{8{\omega }_{2}^{2}}+\frac{{a}_{1}^{2}{a}_{2}}{4{\omega }_{1}^{2}}\right)+{c}_{12}\frac{{a}_{2}-{a}_{4}\mathrm{cos}\left({\theta }_{4}-{\theta }_{2}\right)}{2}-{k}_{12}\frac{{a}_{4}\mathrm{sin}\left({\theta }_{4}-{\theta }_{2}\right)}{2{\omega }_{2}}=0,$
(24)
${m}_{1}\left(2{\omega }_{2}{\stackrel{˙}{\theta }}_{2}+{\omega }_{2}^{2}\right)-{k}_{1}-{k}_{3}\left(\frac{3{a}_{2}^{2}}{4{\omega }_{2}^{2}}+\frac{3{a}_{1}^{2}}{2{\omega }_{1}^{2}}\right)-{c}_{12}\frac{{a}_{4}{\omega }_{2}\mathrm{sin}\left({\theta }_{4}-{\theta }_{2}\right)}{{a}_{2}}$
(25)
${m}_{2}{\stackrel{˙}{a}}_{3}+{c}_{2}\frac{{a}_{3}}{2}+{c}_{4}\left(\frac{{a}_{3}^{3}}{8{\omega }_{1}^{2}}+\frac{{a}_{4}^{2}{a}_{3}}{4{\omega }_{2}^{2}}\right)+{c}_{12}\frac{{a}_{3}-{a}_{1}\mathrm{cos}\left({\theta }_{1}-{\theta }_{3}\right)}{2}-{k}_{12}\frac{{a}_{1}\mathrm{sin}\left({\theta }_{1}-{\theta }_{3}\right)}{2{\omega }_{1}}=0,$
(26)
${m}_{2}\left(2{\omega }_{1}{\stackrel{˙}{\theta }}_{3}+{\omega }_{1}^{2}\right)-{k}_{2}-{k}_{4}\left(\frac{3{a}_{3}^{2}}{4{\omega }_{1}^{2}}+\frac{3{a}_{4}^{2}}{2{\omega }_{2}^{2}}\right)-{c}_{12}\frac{{a}_{1}{\omega }_{1}\mathrm{sin}\left({\theta }_{1}-{\theta }_{3}\right)}{{a}_{3}}$
(27)
${m}_{2}{\stackrel{˙}{a}}_{4}+{c}_{2}\frac{{a}_{4}}{2}+{c}_{4}\left(\frac{{a}_{4}^{3}}{8{\omega }_{2}^{2}}+\frac{{a}_{3}^{2}{a}_{4}}{4{\omega }_{1}^{2}}\right)+{c}_{12}\frac{{a}_{4}-{a}_{2}\mathrm{cos}\left({\theta }_{2}-{\theta }_{4}\right)}{2}-{k}_{12}\frac{{a}_{2}\mathrm{sin}\left({\theta }_{2}-{\theta }_{4}\right)}{2{\omega }_{2}}=0,$
(28)
${m}_{2}\left(2{\omega }_{2}{\stackrel{˙}{\theta }}_{4}+{\omega }_{2}^{2}\right)-{k}_{2}-{k}_{4}\left(\frac{3{a}_{4}^{2}}{4{\omega }_{2}^{2}}+\frac{3{a}_{3}^{2}}{2{\omega }_{1}^{2}}\right)-{c}_{12}\frac{{a}_{2}{\omega }_{2}\mathrm{sin}\left({\theta }_{2}-{\theta }_{4}\right)}{{a}_{4}}$

According to Eqs. (21)-(28), the slow-varying amplitude $a$ and slow-varying phase angle $\theta$ can be acquired. Furthermore, the responses of the original system can be obtained by:

(29)
$y\left(t\right)={y}_{1}\left(t\right)+{y}_{2}\left(t\right)=\frac{{a}_{1}}{{\omega }_{1}}\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t+{\theta }_{1}\right)+\frac{{a}_{2}}{{\omega }_{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{2}t+{\theta }_{2}\right),$
(30)
$z\left(t\right)={z}_{1}\left(t\right)+{z}_{2}\left(t\right)=\frac{{a}_{3}}{{\omega }_{1}}\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t+{\theta }_{3}\right)+\frac{{a}_{4}}{{\omega }_{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{2}t+{\theta }_{4}\right).$

#### 3.3. Numerical experiment of CA method

The numerical experiment is performed by taking the actual structure parameters of a rolling mill. The parameters are:150.9×103 kg, 120.7×103 kg, 10.1×1010 N/m, ${c}_{1}=$2.6×106 N·s/m, 6.9×1010 N/m, 1.5×106 N·s/m, 5.6×1010 N/m, 5×104 N·s/m.

The initial vibration displacement is taken as 1 mm, 0. The initial vibration velocity is taken as $\stackrel{˙}{y}\left(0\right)=0$, $\stackrel{˙}{z}\left(0\right)=0$. The nonlinear stiffness coefficients ${k}_{3}$, ${k}_{4}$ and the nonlinear damping coefficients ${c}_{3}$, ${c}_{4}$ are taken different values. Then, the theoretical solution and CA analytical solution (CAAS) of Eqs. (1) and (2) are calculated by using MATLAB, the results are displayed in Fig. 2.

It can be seen from Fig. 2, for the linear vibration system shown in Fig. 2(a) and the weak NVS displayed in Fig. 2(b), the CAAS well coincides with the theoretical solution. And, for the strong NVS represented in Fig. 2(c), the CAAS also shows good agreement with the theoretical solution. However, as for the strong NVS revealed in Fig. 2(d), because of the nonlinear elastic force $k{y}^{3}$, $k{z}^{3}$, and the nonlinear damping force $c{y}^{2}\stackrel{˙}{y}$, $c{z}^{2}\stackrel{˙}{z}$ containing high-order terms, so the nonlinear factors will have a dominant effect on vibration response. Hence, for the strong NVS which troubled by dominant nonlinear factors, there is a large error exist.

#### 3.4. Error analysis of CA method

In order to evaluate the calculation accuracy of the CA method, it is necessary to search proper evaluation indices. Based on the previous study , this paper takes the average peak error, the average maximum relative error and the error of square sum as evaluation indices.

The average peak error presents the correlation degree of waveform maximum range between the CAAS and the theoretical solution, which is expressed as follows:

(31)
${E}_{ap}=\frac{1}{2}\left\{\left|\frac{\mathrm{m}\mathrm{a}\mathrm{x}\left[X\left(t\right)\right]-\mathrm{m}\mathrm{a}\mathrm{x}\left[S\left(t\right)\right]}{\mathrm{m}\mathrm{a}\mathrm{x}\left[S\left(t\right)\right]}\right|+\left|\frac{\mathrm{m}\mathrm{i}\mathrm{n}\left[X\left(t\right)\right]-\mathrm{m}\mathrm{i}\mathrm{n}\left[S\left(t\right)\right]}{\mathrm{m}\mathrm{i}\mathrm{n}\left[S\left(t\right)\right]}\right|\right\}.$

The average maximum relative error presents the local maximum error of calculation results, which is represented as follows:

(32)
${E}_{ar}=\frac{1}{2}\left\{\left|\frac{\mathrm{m}\mathrm{a}\mathrm{x}\left[X\left(t\right)-S\left(t\right)\right]}{\mathrm{m}\mathrm{a}\mathrm{x}\left[S\left(t\right)\right]}\right|+\left|\frac{\mathrm{m}\mathrm{i}\mathrm{n}\left[X\left(t\right)-S\left(t\right)\right]}{\mathrm{m}\mathrm{i}\mathrm{n}\left[S\left(t\right)\right]}\right|\right\}.$

Fig. 2. Comparison between theoretical solution and CA analytical solution a)0, 0, 0, 0 b)0.01${k}_{1}$, 0.01${k}_{2}$, 0.01${c}_{1}$, 0.01${c}_{2}$ c)0.1${k}_{1}$, ${k}_{4}=$ 0.1${k}_{2}$, ${c}_{3}=$ 0.1${c}_{1}$, ${c}_{4}=$ 0.1${c}_{2}$ d)${k}_{3}={k}_{1}$, ${k}_{4}={k}_{2}$, ${c}_{3}={c}_{1}$, ${c}_{4}={c}_{2}$

The error of square sum can describe the energy error between the CAAS and the theoretical solution, which is depicted as follows:

(33)
${E}_{qs}=\frac{\sum _{i=1}^{N}\left[X\left(i\right){\right]}^{2}-\sum _{i=1}^{N}\left[S\left(i\right){\right]}^{2}}{\sum _{i=1}^{N}\left[S\left(i\right){\right]}^{2}}.$

In Eqs. (31)-(33), $X\left(t\right)$ denotes the sampling value of CAAS, $S\left(t\right)$ is the sampling value of theoretical solution, and $N$ is the number of the sampling points.

The accuracy of the CAAS displayed in Fig. 2 is evaluated by using the above indices, and the results are depicted in table 1. Seen from table 1, as for the linear vibration system displayed in Fig. 2(a) and the weak NVS shown in Fig. 2(b), the CAAS have a high accuracy. Moreover, for the strong NVS described in Fig. 2(c), the results still have high precision. But for the strong NVS represented in Fig. 2(d), which troubled by dominant nonlinear factors, a large deviation appears. Especially, the values of ${E}_{ar}$ and ${E}_{qs}$ are large, and they need to be amended.

Table 1. Errors of CAAS

 Error value Vibration displacement $y$ Vibration displacement $z$ ${E}_{ap}$ (%) ${E}_{ar}$ (%) ${E}_{qs}$ (%) ${E}_{ap}$ (%) ${E}_{ar}$ (%) ${E}_{qs}$ (%) 0, 0, 0, 0 1.11 8.39 9.58 0.91 8.36 9.60 0.01${k}_{1}$, 0.01${k}_{2}$, 0.01${c}_{1}$, 0.01${c}_{2}$ 0.80 8.26 9.66 0.86 8.25 9.68 0.1${k}_{1}$, ${k}_{4}=$ 0.1${k}_{2}$, ${c}_{3}=$ 0.1${c}_{1}$, ${c}_{4}=$ 0.1${c}_{2}$ 0.22 8.44 10.30 1.14 8.62 10.33 ${k}_{3}={k}_{1}$, ${k}_{4}={k}_{2}$, ${c}_{3}={c}_{1}$, ${c}_{4}={c}_{2}$ 0.79 35.34 12.22 4.97 37.02 12.43

#### 4.1. Basic principle of FEMD method

EMD can decompose non-stationary signals into a set of steady data series named intrinsic mode function (IMF). And IMF components must satisfy the following conditions :

(1) Among the entire data set, the number of extrema and the number of zero-crossings must be equal or the difference is no more than one.

(2) At any moment, the mean value of the upper envelope formed by the local maxima and the lower envelope formed by the local minima is zero.

The purpose of the first condition is to remove the superimposed wave. The second condition changes the traditional global limit into the localized limit so that the waveform can be more symmetrical. Virtually, EMD is a screening process. Its detailed steps are as follows:

(1) All of the extreme points of the original signal $s\left(t\right)$ are confirmed. Then, the upper and lower envelopes are calculated according to the extreme points. These two envelopes should envelop all data points.

(2) The mean of the two envelopes is denoted as $\stackrel{-}{\mu }$. Ideally, the first IMF component of $s\left(t\right)$ is achieved as follows:

(34)
${I}_{1}\left(t\right)=s\left(t\right)-\stackrel{-}{\mu }.$

(3) If ${I}_{1}\left(t\right)$ does not meet the IMF conditions, it will be viewed as new $s\left(t\right)$. And the above steps (1)-(2) are repeated until ${I}_{1}\left(t\right)$ corresponds with the IMF conditions.

(4) ${I}_{1}\left(t\right)$ is separated from $s\left(t\right)$, and the remaining portion is a differential signal ${r}_{1}\left(t\right)$ obtained by removing the high-frequency component from $s\left(t\right)$. The expression is as follows:

(35)
${r}_{1}\left(t\right)=s\left(t\right)-{I}_{1}\left(t\right).$

(5) The ${r}_{1}\left(t\right)$ will be treated as a new input signal $s\left(t\right)$. The above steps (1)-(4) are repeated, and the second IMF component ${I}_{2}\left(t\right)$ can be obtained.

(6) When the above steps (1)-(5) are repeated, then the functions from ${I}_{3}\left(t\right)$ to ${I}_{n}\left(t\right)$ can be gained. This process cannot be stopped until ${r}_{n}\left(t\right)$ fulfils a given termination condition when ${r}_{n}\left(t\right)$ is a monotonic function. Then, the decomposition formula of $s\left(t\right)$ can be described as follows:

(36)
$s\left(t\right)=\sum _{j=1}^{n}{I}_{j}\left(t\right)+{r}_{n}\left(t\right),$

among them, ${I}_{1}\left(t\right)$,…, ${I}_{n}\left(t\right)$ are IMF components; ${r}_{n}\left(t\right)$ is the residual component.

The above steps indicate that original signal can be decomposed by EMD to finite IMF components and a residual component. Each IMF component has different frequency-band features, and the residual component is an average trend of original signal. Meanwhile, the IMFs depend on the signal itself. Thus, EMD is a self-adaptive method.

The basic idea of FEMD method is : the number of decomposition mode is beforehand defined; moreover, the data series is divided by partition windows ${W}_{win}$; then, each sub-series thus generated would be sequentially processed by EMD from left to right, as demonstrated in Fig. 3. To ensure the continuity of the decomposed IMFs, the computation window ${W}_{comp}$ is thus defined by extending ${W}_{win}$ on both sides with an overlap percentage $OP$ (01), as described in Eq. (37). FEMD can preset the number of decomposition mode so that it can improve the calculation speed and effectively suppress the mode mixing phenomenon:

(37)
${W}_{comp}=\left(1+2OP\right){W}_{win}.$

Fig. 3. The diagram of FEMD method #### 4.2. Relation between FEMD and CA

Firstly, the theoretical solutions exhibited in Fig. 2(a) are calculated by fast Fourier transform (FFT), the results are shown in Fig. 4. The results indicate that the natural frequencies of system are 126 Hz, 192 Hz, and the corresponding natural angular frequencies are 791.6813 rad/s, 1206.3716 rad/s.

Secondly, the theoretical solutions revealed in Fig. 2(a) are decomposed by FEMD method, and some IMF components are acquired. Then, the phase information ${\varphi }_{j}\left(t\right)$ of each IMF is obtained by using Hilbert transform method. Furthermore, the instantaneous frequency ${f}_{j}\left(t\right)$ of each IMF is achieved according to Eq. (38), and the results are displayed in Fig. 5.

(38)
${f}_{j}\left(t\right)=\frac{1}{2\pi }\frac{d{\varphi }_{j}\left(t\right)}{dt}.$

After calculation, the mean instantaneous frequencies of IMF1, IMF2 decomposed from y are 191.5904 Hz and 125.5786 Hz. The mean instantaneous frequencies of IMF1, IMF2 decomposed from $z$ are 191.9599 Hz and 125.6273 Hz. Compared with the natural frequencies obtained by FFT, the mean instantaneous frequency of IMF1 corresponds to the natural frequency ${f}_{2}$, and the mean instantaneous frequency of IMF2 is corresponding to the natural frequency ${f}_{1}$. So, IMF1 and IMF2 are the inherent mode components of system.

Fig. 4. Amplitude spectrum of vibration displacement signal a) Vibration displacement $y$ b) Vibration displacement $z$

Fig. 5. Time-frequency spectrum of vibration displacement signal a) Vibration displacement $y$ b) Vibration displacement $z$

In order to further explore the connection between FEMD and CA, the IMF components acquired by FEMD are compared with the analytical solutions obtained by CA. The results are represented in Fig. 6 and Fig. 7.

Fig. 6. Comparison between IMF component and CAAS derived from $y$  It can be seen from Fig. 6 and Fig. 7, the CAAS are consistent with the IMF component. So, the amplitude and phase of them should have one-to-one relationship. Hence, when deviation appears in CA method influenced by nonlinear factors, the amplitude and phase of CAAS can be modified by FEMD method.

Fig. 7. Comparison between IMF component and CAAS derived from $z$  The theoretical connection between FEMD method and CA method is further derived. If the IMF components obtained by FEMD are processed by Hilbert transform, there are:

(39)
$H\left[{I}_{j}\left(t\right)\right]=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\frac{{I}_{j}\left(\tau \right)}{t-\tau }d\tau .$

Then, the analytic signal of IMF component can be constructed by ${I}_{j}\left(t\right)$ and $H\left[{I}_{j}\left(t\right)\right]$. Its expression is as follows:

(40)
${Z}_{i}\left(t\right)={I}_{j}\left(t\right)+iH\left[{I}_{j}\left(t\right)\right]={A}_{j}\left(t\right){e}^{i{\varphi }_{j}\left(t\right)},$

where amplitude ${A}_{j}\left(t\right)=\sqrt{{I}_{j}^{2}\left(t\right)+{H}^{2}\left[{I}_{j}\left(t\right)\right]}$, phase ${\varphi }_{j}\left(t\right)=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}H\left[{I}_{j}\left(t\right)\right]/{I}_{j}\left(t\right)$.

According to the above analysis, the relationship between the IMF component and the CAAS can be acquired from Eq. (40), Eq. (29) and Eq. (30). Its expression can be depicted by:

(41)
$A\left(t\right){e}^{i\varphi \left(t\right)}=\frac{a}{\omega }\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega t+\theta \right).$

Moreover, according to Eq. (41) and mathematical transformation , the modified expressions of the slow-varying amplitude $a$ and slow-varying phase angle $\theta$ derived from CAAS can be further obtained. Their expressions can be given by:

(42)
$a=\omega A\left(t\right),$
(43)
$\theta \approx \varphi \left(t\right)-\omega t+\frac{\pi }{2}.$

#### 4.3. Correction for CA analytical solution

Based on the above analysis and theoretical derivation, the detailed steps of the correction for CAAS can be described as follows:

(1) The vibration response of the system is tested to obtain a set of vibration displacement signals.

(2) The vibration displacement signals are individually processed by FFT to acquire the natural frequency $f$ and natural angular frequency $\omega$.

(3) The vibration displacement signals are also individually decomposed by FEMD to achieve a number of IMF components.

(4) The instantaneous frequency of each IMF component is computed and further compared with the natural frequency acquired by FFT. Then, the IMF components which can represent the natural frequency of system are selected.

(5) The selected IMF components are processed by Hilbert transform to obtain the slow-varying amplitude $A\left(t\right)$ and slow-varying phase $\varphi \left(t\right)$.

Fig. 8. Comparison between theoretical solution and CA corrected solution a)0, 0, 0, 0 b)0.01${k}_{1}$, 0.01${k}_{2}$, 0.01${c}_{1}$, 0.01${c}_{2}$ c)0.1${k}_{1}$, ${k}_{4}=$ 0.1${k}_{2}$, ${c}_{3}=$ 0.1${c}_{1}$, ${c}_{4}=$ 0.1${c}_{2}$ d)${k}_{3}={k}_{1}$, ${k}_{4}={k}_{2}$, ${c}_{3}={c}_{1}$, ${c}_{4}={c}_{2}$

(6) The parameters of the slow-varying amplitude $a$ and slow-varying phase angle $\theta$ are corrected by utilizing the modified expressions Eqs. (42) and (43). Afterwards, the corrected $a$, $\theta$ and the natural angular frequency $\omega$ are substituted into the Eqs. (29) and (30). Then, an accurate CAAS can be acquired.

By using the above method, the CAAS described in Fig. 2 are modified, and the corrected results are shown in Fig. 8.

Similarly, the accuracy of the CA corrected solutions depicted in Fig. 8 is evaluated by using the error indices, and the results are displayed in Table 2.

Table 2. Errors of CA corrected solution

 Error value Vibration displacement $y$ Vibration displacement $z$ ${E}_{ap}$ (%) ${E}_{ar}$ (%) ${E}_{qs}$ (%) ${E}_{ap}$ (%) ${E}_{ar}$ (%) ${E}_{qs}$ (%) 0, 0, 0, 0 1.81×10-16 3.07×10-14 3.61×10-16 3.34×10-16 3.20×10-14 2.91×10-16 0.01${k}_{1}$, 0.01${k}_{2}$, 0.01${c}_{1}$, 0.01${c}_{2}$ 0 3.12×10-14 1.98×10-15 6.14×10-17 3.25×10-14 2.18×10-15 0.1${k}_{1}$, ${k}_{4}=$ 0.1${k}_{2}$, ${c}_{3}=$ 0.1${c}_{1}$, ${c}_{4}=$ 0.1${c}_{2}$ 1.11×10-16 3.14×10-14 1.81×10-16 2.50×10-16 3.23×10-14 2.91×10-15 ${k}_{3}={k}_{1}$, ${k}_{4}={k}_{2}$, ${c}_{3}={c}_{1}$, ${c}_{4}={c}_{2}$ 1.89×10-16 3.27×10-14 1.68×10-15 1.64×10-16 3.43×10-14 1.51×10-16

Comparing Table 2 and Table 1, whether for weak NVS or strong NVS, the CA corrected solutions have high precision. The errors between the theoretical solutions and CA corrected solutions are small. Especially, as for the strong NVS (${k}_{3}={k}_{1}$, ${k}_{4}={k}_{2}$, ${c}_{3}={c}_{1}$, ${c}_{4}={c}_{2}$) which is troubled by dominant nonlinear factors, the errors are significantly reduced. The results indicate that the CAAS can be effectively corrected by utilizing FEMD method, and the accuracy of the solution can be significantly improved.

#### 5. Conclusions

On the basis of theoretical research, the nonlinear vertical vibration model of MRS is established. And the analytical solution of the model is explored. Through theoretical research and numerical experiment analysis, the following conclusions are obtained:

(1) The CA method has a high accuracy in solving linear vibration system and weak NVS. As for strong NVS, because of the nonlinear elastic force and the nonlinear damping force containing high-order terms, the nonlinear factors have a dominant effect on system vibration response. So, the CAAS will have a large error while system is troubled by dominant nonlinear factors, which needs to be modified.

(2) The FEMD can decompose non-stationary signals into some IMF components. The amplitude and phase between CAAS and the selected IMF component have specific relationship. Hence, when deviation appears in CA method influenced by nonlinear factors, the amplitude and phase of CAAS can be modified by FEMD method. Results indicate that CAAS can be effectively corrected by utilizing FEMD method, and the accuracy can be significantly improved. The improved CA method can be used to solve strong NVS with two DOF.

#### Acknowledgements

This work is supported by National Natural Science Foundation of China (No. 51475405), National Key Basic Research Program (973 Program) of China (No. 2014CB046405), Innovation Foundation for Graduate Students of Hebei Province (00302-6370002) and Natural Science Foundation of Hebei Province (E2017203039).

1. Hameed W. I., Mohamad K. A. Strip thickness of cold rolling mill with roll eccentricity compensation using fuzzy neural network. Engineering, Vol. 6, 2014, p. 27-33. [Publisher]
2. Tran D. C., Tardif N., Limam A. Experimental and numerical modeling of ﬂatness defects in strip cold rolling. International Journal of Solids and Structures, Vol. 69-70, 2015, p. 343-349. [Publisher]
3. Zhang Y. J., Shen G. X., Li Y. G., et al. Spatial vibration and its numerical analytical method of four-high rolling mills. International Journal of Iron and Steel Research, Vol. 21, Issue 9, 2014, p. 837-843. [Publisher]
4. Zárate L. E., Bittencout F. R. Representation and control of the cold rolling process through artiﬁcial neural networks via sensitivity factors. Journal of Materials Processing Technology, Vol. 197, Issues 1-3, 2008, p. 344-362. [Publisher]
5. Yang X., Tong C. N., Yue G. F., et al. Coupling dynamic model of chatter for cold rolling. International Journal of Iron and Steel Research, Vol. 17, Issue 12, 2010, p. 30-34. [Publisher]
6. Schausberger F., Steinboeck A., Kugi A. Mathematical modeling of the contour evolution of heavy plates in hot rolling. Applied Mathematical Modelling, Vol. 39, Issue 15, 2015, p. 4534-4547. [Publisher]
7. Zhu Y., Jiang W. L., Kong X. D., et al. Study on nonlinear dynamics characteristics of electro-hydraulic servo system. Nonlinear Dynamics, Vol. 80, Issues 1-2, 2015, p. 723-737. [Publisher]
8. Novak N., Vladimir S., Vladimir D. Optimal control of hydraulically driven parallel robot platform based on firefly algorithm. Nonlinear Dynamics, Vol. 82, Issue 3, 2015, p. 1457-1473. [Search CrossRef]
9. Navarro J. F. On the implementation of the Poincaré-Lindstedt technique. Applied Mathematics and Computation, Vol. 195, Issue 1, 2008, p. 183-189. [Publisher]
10. Lakrad F., Belhaq M. Periodic solutions of strongly non-linear oscillators by the multiple scales method. Journal of Sound and Vibration, Vol. 258, Issue 4, 2002, p. 677-700. [Publisher]
11. Zhang C. L., Harne R. L., Li B., et al. Reconstructing the transient, dissipative dynamics of a bistable Dufﬁng oscillator with an enhanced averaging method and Jacobian elliptic functions. International Journal of Non-Linear Mechanics, Vol. 79, 2016, p. 26-37. [Publisher]
12. Cai J. P., Wu X. F., Li Y. P. Comparison of multiple scales and KBM methods for strongly nonlinear oscillators with slowly varying parameters. Mechanics Research Communications, Vol. 31, Issue 5, 2004, p. 519-524. [Publisher]
13. Farzaneh Y., Tootoonchi A. A. Global error minimization method for solving strongly nonlinear oscillator differential equations. Computers and Mathematics with Applications, Vol. 59, 2010, p. 2887-2895. [Publisher]
14. Hosen M. A., Chowdhury M. S. H. A new reliable analytical solution for strongly nonlinear oscillator with cubic and harmonic restoring force. Results in Physics, Vol. 5, 2015, p. 111-114. [Publisher]
15. Razzak M. A., Molla M. H. U. A new analytical technique for strongly nonlinear damped forced systems. Ain Shams Engineering Journal, Vol. 6, Issue 4, 2015, p. 1225-1232. [Publisher]
16. Manevitch L. I. Complex representation of dynamics of coupled nonlinear oscillators. Mathematical Models of Nonlinear Excitations, Transfer Dynamics and Control in Condensed Systems, 1999, p. 269-300. [Publisher]
17. Kerschen G., Vakakis A. F., Lee Y. S., et al. Toward a fundamental understanding of the Hilbert-Huang transform in nonlinear structural dynamics. Journal of Vibration and Control, Vol. 14, Issues 1-2, 2008, p. 77-105. [Publisher]
18. Huang N. E., Shen Z., Long S. R., et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis. Proceedings of Royal Society of London A, Vol. 454, 1998, p. 903-995. [Publisher]
19. Shaw P. K., Saha D., Ghosh S., et al. Investigation of coherent modes in the chaotic time series using empirical mode decomposition and discrete wavelet transform analysis. Chaos, Solitons and Fractals, Vol. 78, 2015, p. 285-296. [Publisher]
20. Gupta P., Sharma K. K., Joshi S. D. Fetal heart rate extraction from abdominal electrocardiograms through multivariate empirical mode decomposition. Computers in Biology and Medicine, Vol. 68, 2016, p. 121-136. [Publisher]
21. Lei Y. G., Lin J., He Z. J., et al. A review on empirical mode decomposition in fault diagnosis of rotating machinery. Mechanical Systems and Signal Processing, Vol. 35, Issues 1-2, 2013, p. 108-126. [Publisher]
22. Xie X. H. Illumination preprocessing for face images based on empirical mode decomposition. Signal Processing, Vol. 103, 2014, p. 250-257. [Publisher]
23. Fan G. F., Peng L. L., Hong W. C., et al. Electric load forecasting by the SVR model with differential empirical mode decomposition and auto regression. Neurocomputing, Vol. 173, 2016, p. 958-970. [Publisher]
24. Wu Z., Huang N. E. Ensemble empirical mode decomposition: a noise-assisted data analysis method. Advances in Adaptive Data Analysis, Vol. 1, Issue 1, 2009, p. 1-41. [Publisher]
25. Wang Y. H., Yeh C. H., Young H. V., et al. On the computational complexity of the empirical mode decomposition algorithm. Physica A, Vol. 400, Issue 15, 2014, p. 159-167. [Publisher]
26. Hou D. H., Chen H., Liu B., et al. Analyze on nonlinear parametrically exited vibrations of roller system. Journal of Vibration and Shock, Vol. 28, Issue 11, 2009, p. 1-5. [Search CrossRef]
27. Meng L. Q., Du Y., Ma S. B., et al. Non-linearity of vertical vibration of medium plate mill. Journal of Jilin University, Vol. 39, Issue 3, 2009, p. 712-715. [Search CrossRef]
28. Chen Y. H., Shi T. L., Yang S. Z. Study on parametrically excited nonlinear vibrations on 4-H cold rolling mills. Journal of Mechanical Engineering, Vol. 39, Issue 4, 2003, p. 56-60. [Publisher]
29. Liu H. R., Shi P. M., Chen H., et al. Study on nonlinear parametrically exited coupling vibrations of roller system on 4-H rolling mills. China Mechanical Engineering, Vol. 22, Issue 12, 2011, p. 1397-1401. [Search CrossRef]
30. Zhu Y., Jiang W. L., Kong X. D., et al. An accurate integral method for vibration signal based on feature information extraction. Shock and Vibration, Vol. 2015, 2015, p. 1-13. [Publisher]

#### Cited By

 Processes Yong Zhu, Shengnan Tang, Chuan Wang, Wanlu Jiang, Xiaoming Yuan, Yafei Lei 2019 AIP Advances Yong Zhu, Pengfei Qian, Shengnan Tang, Wanlu Jiang, Wei Li, Jianhua Zhao 2019