Abstract
According to the movement mechanism of strip and rollers during the continuous rolling process, the main drive system of each stand was simplified to a single degree of freedom discrete model, and the strip was simplified to an axially moving Euler beam. Then, a nonlinear continuousdiscrete coupled vibration model between transverse and longitudinal vibrations of strip and torsional vibration of main drive system was established. According to Hamilton’s principle, the nonlinear differential equations were established. Moreover, modified iteration method and Kantorovich averaging method were used to solve the differential equations. Depending on numerical calculation, the amplitudefrequency responses of strip vibration coupled with torsional vibration of main drive system were obtained. Finally, the influences of the axial velocity, the strip tension, the torsional stiffness, and the rotational inertia on the vibration characteristics were discussed. The results would provide a theoretical reference for control and analysis of rolling mill vibration in engineering practice.
1. Introduction
The rolling mill vibration has brought serious troubles to iron and steel industry for many years, which is a significant technical problem. The vibration forms are various, and characteristics and causes are different. At present, the research on the different forms of rolling mill vibration mainly includes vertical vibration and horizontal vibration of rollers, transverse and longitudinal vibration of strip [1, 2], torsional vibration of main drive system [3] and the axial vibration, etc. As early as in 1967, the research of various vibrations of rolling mill was begun. Moller and Hoggart et al. [4] found the existence of torsional vibration of the rolling mill when they were employing the two rolls test machine, and the vibration behavior, which was considered as the selfexcited vibration. Lawrence and Thomas [5] analyzed the chattermark on the strip surface, which caused torsional vibration by defect of transmission gear. Swiatoniowski and Bar [6] analyzed the chattering phenomenon during the rolling process, and the mathematical model of selfexcited vibration was established. Wang et al. [7] analyzed the influence of rolling force on horizontal stiffness of rolling mill, and the functional relationship among the vertical vibration of rolling mill, the horizontal vibration of work roll and the torsional vibration of main drive system was established. Sun et al. [8] investigated the torsional vibration on the dynamic model of six rollers system of 1100 rolling mill, and the influence of torsional vibration on strip shape was acquired. Xu et al. [9, 10] established the hybrid system model of strip coupled with rollers, and the coupled vibration mechanism of roller and strip with tension fluctuation was obtained. Zou et al. [11] analyzed the different influences of cold rolling and hot rolling on axial force, and a mathematical model was established. Du et al. [12] established the nonlinear dynamic model of coupled model between rollers and strip, and the inertial boundary condition was proposed.
However, there are few researches studying on coupled vibration of multiple systems. In the earlier stage, the research of strip vibration coupled with torsional vibration of main drive system has been investigated by our group, such as the vibration analysis of main drive system of each stand by finite element software. Researchers found the energy was transferred from one to another main drive system through the strip when resonance occurred, which resulted to the great fluctuation of rollers, and the vibration amplitude may exceed selfresonance [13]. Therefore, it is significant to analyze the coupled vibration between strip and main drive system of rolling mill furtherly.
In this paper, the main drive system of rolling mill is simplified to a single degree of freedom model of discrete system, and the strip is simplified to an axially moving Euler beam [14, 15]. Then, a continuousdiscrete coupled vibration model between transverse and longitudinal vibrations of strip and torsional vibration of rolling mill main drive system is established. The research results can provide important theoretical reference for control and analysis of strip vibration coupled with torsional vibration during the continuous rolling process.
2. Mechanical and mathematical models
The main drive system of each stand can be regarded as a massspring system in the tandem mill, which mainly composed of inertial components of electromotor, reduction gears and rollers etc., and elastic components of connecting shafts, a dimensional geometric model is shown in Fig. 1. Due to the traction behavior of motor drive, the motor of main drive system could be regarded as a fixed end [16]. And the rotational inertia of roller is much bigger than the other parts, so the main drive system can be simplified to a single degree of freedom discrete model of springmass system. The rollers can be equivalent to the rigid inertial components, $j$ is the rotational inertia and $\theta $ is the rotation angle; the connecting shafts can be equivalent to the elastic element, $k$ is the torsional stiffness, as shown in Fig. 2. When transverse and longitudinal vibrations of strip are considered, the mechanical geometry characteristic is similar to an axially moving beam [17, 18]. Also, the strip is an elastic continuum whose thickness is much less than the length, the strip is equivalent to the isotropic Euler beam based on the theory of axially moving beam. The equivalent mechanical model is shown in Fig. 3, and assumes that the strip with uniform motion is running, ${v}_{0}$is the axial velocity. It hypothesizes that there is no relative motion between rollers and strip, and the upper and lower rollers are symmetric discs of fixed axis rotation along the width of beam. The transverse displacement and the longitudinal displacement of Euler beam are $w({x}_{0},{y}_{0},t)$ and $u({x}_{0},{y}_{0},t)$ respectively; $l$ is the distance of Euler beam between two stands; the left tension and the right tension are ${P}_{1}$ and ${P}_{2}$ respectively, and ${P}_{1}={P}_{2}={P}_{0}$.
Fig. 1A dimensional geometric model of main drive system of one stand: 1 – electromotor, 2 – intermediate shaft, 3 – reduction gears, 4 – profile spindle, 5 – working roll, 6 – backup roll, 7 – screwdown device, 8 – balancing device, 9 – unjamming gear, 10 – mill housing, 11 – platform
Based on Hamilton principle, the mathematical model of strip vibration coupled with torsional vibration of main drive system is established.
The kinetic energy ${T}_{1}$ of axially moving Euler beam can be written as:
where $\rho $ is the density of strip, $A$ is the crosssectional area of Euler beam, ${u}_{,t}$, ${w}_{,t}$ and ${u}_{,{x}_{0}}$, ${w}_{,{x}_{0}}$ are the $u({x}_{0},{y}_{0},t)$ and $w({x}_{0},{y}_{0},t)$ on variables $t$ and ${x}_{0}$ of the first partial derivatives respectively.
Fig. 2A simplified model of main drive system of rolling mill
Fig. 3A mechanical model of strip between main drive systems of two stands
The kinetic energy ${T}_{2}$ of roller can be written as:
where $r$ is the radius of roller.
The potential energy ${U}_{1}$ of Euler beam can be written as:
where ${\epsilon}_{{x}_{0}}$ is the strain of Euler beam, $E$ is the Young’s modulus and $I$ is the moment of inertia.
The potential energy ${U}_{2}$ of tension can be written as:
Then according to Hamilton principle, the following equation can be obtained:
Substitution of Eqs. (14) into Eq. (5), the motion equations can be expressed as:
In general, the axial kinetic energy caused by transverse vibration is relatively small. Let ${u}_{,t}={u}_{,tt}=0$, So, Eq. (6) is simplified as:
The relationship between torque and torsional stiffness is expressed as: ${M}_{i}={k}_{i}\mathrm{\Delta}{\theta}_{i}$, where ${\theta}_{i}$ is relative torsional angle between two shaft sections, ${k}_{i}$ is the torsional stiffness of shaft segment. Therefore, boundary conditions can be written as:
When ${x}_{0}=0$:
When ${x}_{0}=l$:
The time variable of equation is separated from the space variable by substitution of $w={\varphi}_{0}\left({x}_{0}\right)\mathrm{c}\mathrm{o}\mathrm{s}{\omega}_{0}t$ into Eq. (8), the following equation can be obtained as: $u={\phi}_{0}\left({x}_{0}\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}{\omega}_{0}t$, where ${\omega}_{0}$ is vibration frequency.
Based on Kantorovich averaging method on the interval $\left[0,2\pi /{\omega}_{0}\right]$, the time variable can be eliminated. The motion equations can be simplified as:
$\frac{3}{4}EA\left({{\phi}_{0}}_{,{x}_{0}{x}_{0}}{{\varphi}_{0}}_{,{x}_{0}}+{{\varphi}_{0}}_{,{x}_{0}{x}_{0}}{{\phi}_{0}}_{,{x}_{0}}+\frac{3}{2}{{\phi}_{0}}_{,{x}_{0}}^{2}{{\phi}_{0}}_{,{x}_{0}{x}_{0}}\right)=0.$
Boundary conditions can be written as, when ${x}_{0}=0$:
When ${x}_{0}=l$:
The dimensionless quantities are given by:
$S=\frac{A{l}^{2}}{I},J=\frac{j}{\rho A{r}^{2}l},P=\frac{{P}_{0}{l}^{2}}{EI},K=\frac{kl}{EA{r}^{3}}.$
Then, the dimensionless form of Eqs. (1318) can be obtained respectively:
When $x=0$:
When $x=1$:
3. Analytical solution of vibration equations
Due to more complicated solving process of the Eqs. (1924), the modified iteration method is used to solve the equations in this section.
3.1. The firstorder approximate solution
Firstly, all the nonlinear terms of the Eq. (20) are omitted, the simplified motion equation can be obtained as follow:
The series solution of Eq. (25) is:
where:
By substitution of Eq. (26) into Eqs. (2324), one has ${\omega}_{1}=$ 18.71. Then, the coefficients of Eq. (26) have:
Then:
where:
And substituting of Eq. (27) into Eq. (19), one can be got:
The coefficients ${c}_{1}$ and ${c}_{2}$ can be determined by Eqs. (2122), and the results are as follows:
Then, by substitution of ${c}_{1}$ and ${c}_{2}$ into ${\varphi}_{1}\left(x\right)$, the firstorder approximate solution can be obtained.
3.2. The secondorder modifiediterative solution
In this subsection, the secondorder modifiediterative solution can be given.
By substitution of Eqs. (2728) into Eq. (20), one has:
That is:
where:
where:
${A}_{n}^{\left(1\right)}=\sum _{{n}_{1}=0}^{n}\frac{{\mu}_{1}^{2}{\omega}_{1}^{2n}}{\left(4{n}_{1}\right)!\left(4n4{n}_{1}\right)!}+\sum _{{n}_{1}=0}^{n1}\frac{{\mu}_{2}^{2}{\omega}_{1}^{2n2}}{\left(4{n}_{1}+2\right)!\left(4n4{n}_{1}2\right)!},\left(n=\mathrm{1,2}\dots \right),$
${B}_{n}^{\left(1\right)}=\sum _{{n}_{1}=0}^{n}\frac{2{\mu}_{1}{\mu}_{2}{\omega}_{1}^{2n}}{\left(4{n}_{1}\right)!\left(4n4{n}_{1}+2\right)!},\left(n=\mathrm{0,1}\dots \right),$
${C}_{n}^{\left(1\right)}=\frac{{\mu}_{1}{\omega}_{1}^{2n}}{\left(4n1\right)!}\begin{array}{ll},& \left(n=\mathrm{1,2}\dots \right)\end{array},{D}_{n}^{\left(1\right)}=\frac{{\mu}_{2}{\omega}_{1}^{2n}}{(4n+1)!}\begin{array}{ll},& \left(n=\mathrm{0,1}...\right)\end{array}.$
With the above equations, that is obtained as follow:
where:
${B}_{n}^{\left(2\right)}=\sum _{m=0}^{n}\left({A}_{n}^{\left(1\right)}{D}_{nm}^{\left(1\right)}+{B}_{n}^{\left(1\right)}{C}_{nm}^{\left(1\right)}\right),\left(n=\mathrm{1,2}\dots \right),$
${C}_{n}^{\left(2\right)}=\sum _{m=0}^{n}{B}_{n}^{\left(1\right)}{D}_{nm}^{\left(1\right)},\left(n=\mathrm{0,1}...\right).$
Due to the property of series solution, the solution of Eq. (32) can be written as:
where ${\zeta}_{1}$ and ${\zeta}_{2}$ denote the undetermined coefficients:
By substitution of Eq. (34) into Eqs. (2324), one can be get:
where:
where:
${d}_{13}=\left(\sum _{n=1}^{\infty}{A}_{n}+\sum _{n=0}^{\infty}{B}_{n}\right)\left.+{\phi}_{m}^{2}\left(\sum _{n=1}^{\infty}{C}_{n}+\right.\sum _{n=0}^{\infty}{D}_{n}+\sum _{n=0}^{\infty}{E}_{n}\right),$
${d}_{21}=\sum _{n=0}^{\infty}\frac{{\omega}^{2n}}{\left(4n+1\right)!}{\left(\frac{1}{2}\right)}^{4n+1},{d}_{22}=\sum _{n=0}^{\infty}\frac{{\omega}^{2n}}{\left(4n+3\right)!}{\left(\frac{1}{2}\right)}^{4n+3},$
$\begin{array}{l}{d}_{23}=\left(\sum _{n=1}^{\infty}{A}_{n}{\left(\frac{1}{2}\right)}^{4n1}+\sum _{n=0}^{\infty}{B}_{n}{\left(\frac{1}{2}\right)}^{4n+1}\right)\\ +{{\phi}_{m}}^{2}\left[\sum _{n=1}^{\infty}{C}_{n}{\left(\frac{1}{2}\right)}^{4n1}+\sum _{n=0}^{\infty}{D}_{n}{\left(\frac{1}{2}\right)}^{4n+1}\right.\left.+\sum _{n=0}^{\infty}{E}_{n}{\left(\frac{1}{2}\right)}^{4n+3}\right]1,\end{array}$
${d}_{31}=\sum _{n=1}^{\infty}\frac{{\omega}^{2n}}{\left(4n1\right)!},{d}_{32}=\sum _{n=0}^{\infty}\frac{{\omega}^{2n}}{\left(4n+1\right)!},$
The frequency ${\omega}_{2}$ can be obtained by $\mathrm{d}\mathrm{e}\mathrm{t}D=$ 0. The coefficients ${\zeta}_{1}$ and ${\zeta}_{2}$ can be determined by Eq. (35), and then, the secondorder modifiediterative solution can be obtained.
4. Numerical simulation and analysis
To verify the proposed method, some numerical simulations are carried out by GLE36 strip. The main parameters of tandem mill and strip are listed as follows: the distance between stand F2 and stand F3 $l=$ 2.5 m, the thickness of strip $h=$ 0.018 m, the torsional stiffness $k=$ 9×10^{3} N·m/rad, and the strip tension ${P}_{0}=$ 8×10^{3} N. When the rotational inertia of roller $j$ is 200 kg·m^{2}, 400 kg·m^{2}, 600 kg·m^{2} and 800 kg·m^{2}, respectively, the relationship curves between axial velocity and frequency of strip are shown in Fig. 4. It can be seen that the vibration frequency decreases with the increasing axial velocity gradually, and the greater the rotational inertia of roller, the faster the frequency decreases. In particular, the vibration frequency becomes more obvious, when ${v}_{0}<$ 0.01 m/s, greater rotational inertia and larger vibration frequency, and vice versa. That is to say, the axial velocity has strong influence on vibration characteristic of strip while the rotational inertia of roller is larger than the others.
Fig. 4Relationship between axial velocity and frequency
Fig. 5Relationship between tension and frequency
Fig. 5 shows the relationship curves between strip tension and frequency of strip with different values of rotational inertias, which are 200 kg·m^{2}, 400 kg·m^{2 }and 600 kg·m^{2}, respectively. The torsional stiffness $k=$ 9×10^{3} N·m/rad, the axial velocity ${v}_{0}=$ 0.008 m/s. From Fig. 5, the vibration frequency decreases with the increase of the strip tension. When ${P}_{0}\le $ 5×10^{3} N, as long as the rotational inertia becomes larger, the frequency of strip gets stronger. Whereas, when ${P}_{0}>$ 5×10^{3} N, the ${P}_{0}{\omega}_{0}$ curves present fluctuation, then stabilize gradually. Thus, it follows that a smaller strip tension change results in a greater effect of frequency. Additionally, it is shown that the bigger the rotational inertia is, the faster the trend of ${P}_{0}{\omega}_{0}$ curve declines. In other words, a further proof is given that a bigger variation of frequency caused by a greater rotational inertia.
Fig. 6 depicts the relationship curves between torsional stiffness of drive system and frequency of strip with different values of rotational inertias, which are 200 kg·m^{2}, 400 kg·m^{2}, 600 kg·m^{2}, and 800 kg·m^{2}, respectively. The strip tension ${P}_{0}=$ 8×10^{3} N, the axial velocity ${v}_{0}=$ 0.008 m/s. As can be seen from the Fig. 6, with the increase of the torsional stiffness, the $k{\omega}_{0}$ curves show a uniform change, and the vibration frequency decreases gradually. And a larger rotational inertia leads to a larger downtrend of $k{\omega}_{0}$ curve. The result shows that the boundary of strip is close to the clamped condition with the infinite increase of rotational inertia. Compared with a model of moving strip without main drive system, the overall frequency curves tend to be more stable. Moreover, the greater the torsional stiffness is, the smaller the rolling forces exerts, which reduces the vibration frequency of strip to some extent.
From the frequencyresponse curves in Figs. 46, a common feature can be found that a larger rotational inertia of roller gives rise to a greater variation tendency of frequency, which causes a more obvious vibration. In other words, a larger rotational inertia has a greater effect on the relationships from the influencing factors that is the axial velocity, strip tension and torsional stiffness to the frequency. Therefore, it is very important to select the appropriate rotational inertia of roller in the research and practice.
Fig. 7 illustrates the amplitudefrequency responses under the condition of different axial velocities: 2×10^{–3} m/s, 2×10^{–3} m/s, 10×10^{–3} m/s, 14×10^{–3} m/s and 18×10^{–3} m/s. The torsional stiffness $k=$ 4×10^{4} N·m/rad, the strip tension ${P}_{0}=$ 8×10^{3} N, and the rotational inertia $j=$ 800 kg·m^{2}. From Fig. 7, it is obtained that with the raise of the axial velocity, the vibration performance is transformed from sclerotic type into intenerate type gradually. When ${v}_{0}\le $ 10×10^{–3} m/s, with the increase of the amplitude, the vibration frequencies increase, and the hardening degree increases with the decreasing axial velocity. When ${v}_{0}\ge $ 14×10^{–3} m/s, the curve shows a trending of decrease, namely, the vibration performance is intenerated type. That is, being oversized or excessively small of the axial velocity will cause the great influence of the amplitude on frequency.
Fig. 6Relationship between torsional stiffness and frequency
Fig. 7Amplitudefrequency response with different axial velocities
Fig. 8 demonstrates the amplitudefrequency responses under the condition of different strip tensions: 4×10^{3} N, 6×10^{3} N, 8×10^{3} N and 10×10^{3} N. The torsional stiffness $k=$ 4×10^{4} N·m/rad, the rotational inertia $j=$ 800 kg·m^{2}, and the axial velocity ${v}_{0}=$ 0.01 m/s. In Fig. 8, with the increase of the amplitude, the vibration frequency increases, and the vibration performance keeps consistent. It also can be known that the changing of vibration frequency is small and stable under the condition of lower amplitude. However, at the point of amplitude ${\phi}_{0}=$ 0.006, there is an intersection. After that, with the raise of the amplitude, the vibration frequency increases fast. That is to say, a smaller strip tension results in a larger variation range of the amplitudefrequency response. Based on the mechanics of vibration, the smaller the tensile force, the more unstable the system, thereby, it makes strip vibration quite strong.
Fig. 8Amplitudefrequency response with different strip tensions
Fig. 9Amplitudefrequency response with different torsional stiffness
Fig. 9 displays the amplitudefrequency responses under the condition of different torsional stiffness: 2×10^{4} N·m/rad, 4×10^{4} N·m/rad, 6×10^{4} N·m/rad and 8×10^{4} N·m/rad. The rotational inertia $j=$ 800 kg·m^{2}, the strip tension ${P}_{0}=$ 8×10^{3} N, and the axial velocity ${v}_{0}=$ 0.01 m/s. From Fig. 9, as the torsional stiffness increases, the vibration performance changes from positive correlation to negative correlation gradually. When $k\le $ 4×10^{4} N·m/rad, the frequencies increase with the increasing amplitude, and the vibration performance shows positive correlation. Keep on increasing the size of torsional stiffness, the amplitude can be decreased, and the vibration performance shows negative correlation. Due to the above analysis, unsuitable size of the torsional stiffness will lead to inconsistent movement between rollers and strip, therefore, it leads to the torsional vibration of main drive system. Consequently, the torsional stiffness should be selected appropriately in practical design to reduce the rolling instability.
5. Conclusions
In this paper, the strip vibration coupled with torsional vibration of main drive system of rolling mill is investigated. It is given that the influences of the axial velocity, strip tension and torsional stiffness of drive system on vibration frequency of strip, in which all of the influencing factors are related to the rotational inertia of roller. A larger rotational inertia of roller gives rise to a greater changing trend of frequency, and leads to a more obvious vibration of strip. Moreover, the axial velocity, the strip tension and the torsional stiffness of drive system have also great influences on the amplitudefrequency characteristic of strip vibration. The changing trends of the amplitudefrequency response curves are accordantly increasing with the different strip tension, and a larger tension results in a smaller increasing gradient. Whereas, the amplitudefrequency characteristic of strip vibration will transform from sclerotic type into intenerate type with the increasing axial velocity or with the increasing torsional stiffness. Therefore, the analysis results of this paper show that the importance of choosing appropriate axial velocity, strip tension, torsional stiffness and rotational inertia on control and optimization of strip vibration of rolling mill.
References

Hou F. X., Zhang J., Cao J. G. Review of chatter studies in cold rolling. Journal of Iron and Steel Research, Vol. 19, Issue 10, 2007, p. 611, (in Chinese).

Wright J. Mill drive system to minimize torque amplification. Iron and Steel Engineer, Vol. 53, Issue 7, 1976, p. 5660.

Fan X. B., Zang Y., Wang F., et al. Hot strip mill nonlinear torsional vibration with multistand coupling. Journal of Vibroengineering, Vol. 17, Issue 4, 2015, p. 16231633.

Moller R. H., Hoggart J. S. Periodic surface finish and torque effects during cold strip rolling. Journal of the Australian Institute of Metals, Vol. 12, Issue 2, 1967, p. 155164.

Lawrence A. B., Thomas S. R. Winding reel involvement in temper mill chatter. Iron and Steel Engineer, Vol. 71, Issue 11, 1994, p. 2729.

Swiatoniowski A., Bar A. Parametrical excitement vibration in tandem millsmathematical model and its analysis. Journal of Materials Processing Technology, Vol. 134, Issue 2, 2003, p. 214224.

Wang R. P., Peng Y., Zhang Y. Mechanism research of rolling mill coupled vibration. Journal of Mechanical Engineering, Vol. 49, Issue 12, 2013, p. 6671, (in Chinese).

Sun J. L., Peng Y., Liu H. M. Dynamic characteristics of cold rolling mill and strip based on flatness and thickness control in rolling process. Journal of Central South University, Vol. 21, Issue 2, 2014, p. 567576.

Wang Y., Xu F., Li Y. Dynamic modeling and coupling vibration analysis of hybrid systems consisting of strip, rolls and flexible supports. Journal of Vibration Engineering, Vol. 26, Issue 4, 2013, p. 599607, (in Chinese).

Xu F., Li Y., Li L. F. Parametric vibration analysis for sheet metal and its supporting system. Mechanical Research and Application, Vol. 26, Issue 3, 2013, p. 2730, (in Chinese).

Zou J. X., Xu L. J. Tandem Mill Vibration Control. Metallurgical Industry Press, Beijing, 1998, (in Chinese).

Gao C. Y., Du G. J., Feng Y., et al. Nonlinear vibration analysis of moving strip with inertial boundary condition. Mathematical Problems in Engineering. Vol. 8, Issue 2015, 2015, p. 19.

Gao C. Y., Du G. J., Li J. X. Numerical analysis on torsional vibration of main driving system of rolling mill with strip steel. Journal of Yanshan University, Vol. 40, Issue 1, 2016, p. 5865, (in Chinese).

Chen L. Q., Yang X. D. Steadystate response of axially moving viscoelastic beam with pulsating speed: comparison of two nonliner models. International Journal of Solids and Structures, Vol. 42, 2005, p. 1015.

Du G. J., Ma J. Q. Nonlinear vibration and buckling of circular sandwich plate under complex load. Applied Mathematics and Mechanics, Vol. 28, Issue 8, 2007, p. 10811091.

Shen Y. Z., Liu H. M., Xiong J., and et al. Analysis of chaotic behavior of the main drive system with clearance of a heavey plate mill. Engineering Mechanics, Vol. 27, Issue 7, 2010, p. 232236, (in Chinese).

Han S. M., Benaroya H., Wei T. Dynamics of transversely vibrating beam using four engineering theories. Journal of Sound and Vibration, Vol. 225, Issue 5, 1999, p. 935988.

Huang J. L., Chen S. H. Study on nonlinear vibration of an axially moving beam with coupled transverse and longitudinal motions. Journal of Vibration and Shock, Vol. 30, Issue 8, 2011, p. 2427, 50, (in Chinese).
Cited by
About this article
This research was supported in part by Natural Science Foundation of Hebei Province (E2017203115), in part by Key Project of Science and Technology of Hebei Higher School (No. ZD2015077), and in part by Doctor Foundation of Yanshan University (No. B992).