Abstract

In order to understand the nonlinear dynamic behavior of a planar mechanism with clearance, the nonlinear dynamic model of the 2-DOF nine-bar mechanism with a revolute clearance is proposed; the dynamic response, phase diagrams, Poincaré portraits, and largest Lyapunov exponents (LLEs) of mechanism are investigated. The nonlinear dynamic model of 2-DOF nine-bar mechanism containing a revolute clearance is established by using the Lagrange equation. Dynamic response of the slider’s kinematics characteristic, contact force, driving torque, shaft center trajectory, and the penetration depth for 2-DOF nine-bar mechanism are all analyzed. Chaos phenomenon existed in the mechanism has been identified by using the phase diagrams, the Poincaré portraits, and LLEs. The effects of the different clearance sizes, different friction coefficients, and different driving speeds on dynamic behavior are studied. Bifurcation diagrams with changing clearance value, friction coefficient, and driving speed are drawn. The research could provide important technical support and theoretical basis for the further study of the nonlinear dynamics of planar mechanism.

1. Introduction

Clearance joint is unavoidable, thanks to many uncertainties such as manufacturing tolerances, assemblage, material deformation, and wear. The clearance of kinematic pair in mechanical system is universal, and the joint clearance will obviously degrade working performance of the mechanical system. [13]. It is known that the nonlinear characteristic is one of most important dynamic behaviors of the nonlinear dynamic system, and it can well describe the movement rule of the system [46]. Therefore, it is necessary to investigate the nonlinear dynamic behavior of mechanism with clearance, and it also provides important theoretical basis and technical support for the further research of nonlinear dynamics of planar mechanism.

In the recent years, the dynamic response of mechanism containing clearance was studied by a lot of scholars. Reis et al. [7] built the dynamic equation of slider-crank mechanism containing revolute joint clearance. The contact force model based upon the Hertz model considers the dissipative influence related to the contact between the piston and pin. Frictional force is based upon the Coulomb friction, but it is applied to multibody dynamics. Wang et al. [8] investigated dynamic characteristics of a slider-crank with revolute clearances and modified contact force model has been built, and friction force of the clearance joint has been also described by utilizing an improved Coulomb friction model. Li et al. [9] conducted a numerically comparative study on dynamic response of a crank-slider mechanism with clearances thinking about harmonic drive and link flexibility and also investigated the effect of harmonic drive on the dynamic behavior of the mechanism with different locations of single clearance joint and with two clearance joints, respectively. Muvengei et al. [10] explored the influences of differently located frictionless joints with clearance on overall dynamic behavior of the crank-slider mechanism. Marques et al. [11] proposed a novel method to build spatial revolute joint containing the axial clearance and the radial clearance. The motion equation which governs dynamic response of multibody system has been established. The Newton–Euler approach was adapted to solve the dynamic equation. Zhang et al. [12] developed the crank-slider mechanism with a revolute clearance joint at the piston and pin and also presented the simulation model by using the ADAMS. The effects of difference clearance value, crank speed, and contact stiffness on dynamic response are all studied. Varedi et al. [13, 14] put forward an optimization method based on PSO to reduce the adverse effects of the clearance joint. The algorithm solved the nonlinear optimization problem for the crank-slider mechanism containing clearances. Flores et al. [15] proposed and discussed a new method based on the nonsmooth dynamics approach for the dynamic modeling and analysis of the crank-slider system containing translational clearance. Erkaya et al. [16] carried out experimental as well as numerical researches to discuss influences of clearance on conventional articulated and partly compliant mechanisms. Abdallah et al. [17] studied the dynamic behavior for the crank-slider mechanism with the flexible links and the multijoint clearances by using ADAMS software. Yang et al. [18] proposed a novel computing method for dynamic analysis of a four-bar mechanism containing multiple joints with clearance by using the finite element method. Zhang et al. [19] analyzed the dynamic response of the 3-RRR mechanism with multiple revolute clearances. The equations of motion have been built by using the Newton–Euler method with the modified Coulomb friction force model and Lankarani–Nikravesh contact force model. Xu and Li [20] researched the influence of clearances on dynamic response of the 2-DOF pick-and-place mechanism containing revolute clearances. The influence of the revolute clearances on contact forces, driving torques, and the slider accelerations are all studied. Varedi-Koulaei et al. [21] investigated influences of the clearance on dynamics response of a 3-RRR mechanism. At the same time, a comprehensive optimization method for kinematics and dynamics is implemented, which decreases undesirable influences and improves mechanism’s accuracy and performance.

Most previous studies researched nonlinear characteristic of the mechanical system with joint clearances and mainly concentrated upon the dynamics study of the simple mechanism. Flores [22] proposed and discussed a method for dynamic analysis and modeling of rigid multibody systems including clearances. The dynamic response and the Poincaré portraits are studied. Rahmanian and Ghazavi [23] researched nonlinear dynamic behavior of the crank-slider mechanism including clearance joints. Through the Poincaré map relating to nonautonomous systems, nonlinear dynamic behavior of crank-slider has been exhibited in discrete space, which is the standard of the chaos detection. Farahan et al. [24] investigated the nonlinear dynamic behavior of the four-bar mechanism containing clearance between the connecting link and rocker. Tang et al. [25] explored the nonlinear dynamic behavior of the four-bar mechanism including the revolute clearance. It is certified that the chaos and the strange attractors exist in dynamic behavior by adopting the Poincaré portrait. Hou et al. [26] took the two-rotation decoupled parallel mechanism RU-RPR; considering clearance as a studied object, the equation of motion of RU-RPR mechanism has been built, and then the chaos and impact phenomena of the mechanism have been investigated. The effects of friction coefficients and driving speeds on chaos and impact phenomena have been all analyzed. Flores et al. [27] used the general methodology to evaluate the effect of the clearance value and friction coefficient on the dynamic response of the rigid crank-slider mechanism containing clearance joints. By utilizing the Poincaré maps, the periodic and the chaotic responses of the mechanism are all studied. Ma and Qian [28] presented a general procedure for dynamic simulation as well as modeling of a slider-crank mechanism including clearances. Effectiveness of the proposed methodology was validated by comparing with ADAMS simulation results of the crank-slider mechanism containing clearances. The Poincaré maps have been used to study the nonlinear characteristics, and influences of different clearance sizes and different driving speeds on dynamic response have been both investigated.

The main aim of this paper is to study the nonlinear dynamic behavior of the multilink complex mechanism with multi-DOF. Taking 2-DOF nine-bar mechanism as the study object, dynamic response of the kinematics characteristic of the slider, contact force, driving torque, shaft center trajectory, and the penetration depth for 2-DOF nine-bar mechanism are all studied, and chaos phenomenon existed in the mechanism is systematically identified by the phase diagrams, the Poincaré portraits, and LLEs, respectively. The effects of clearance value, friction coefficient, and driving speed on nonlinear characteristics of clearance joint are systematically and deeply studied by Poincaré map, Lyapunov exponent, and bifurcation diagram. The study is meaningful for nonlinear dynamics analysis and control of similar multilink mechanisms.

The mechanism consists of the frame, the crank 1, link 2, link 3, the crank 4, rocking-bar 6, triangular panel 7, link 8, and slider 9. The link 5 is a part of the frame. The structure diagram of the planar nine-bar mechanism containing a revolute clearance is shown in Figure 1. Crank 1 and link 2 are connected by a clearance revolute joint which is denoted by . Local enlarged drawing of is shown in Figure 1. The motion pairs of the planar mechanism are made up of revolute pairs and translational pair, and among them, the number of revolute pairs is much higher than translational pair. Therefore, it is more representative to study the dynamic behaviors of the mechanism with revolute clearance. It is remarkable that revolute clearance joint introduces two extra DOFs which contain horizontal and vertical displacements of center of shaft corresponding to the bearing center. The DOF of the mechanism is changed from 2 to 4.

3. Dynamic Model of the Planar Nine Bars Mechanism with a Revolute Clearance

3.1. The Establishment of the Contact Force Model

As shown in Figure 2, eccentricity vector connecting bearing’s center and shaft’s center is expressed as follows:

The penetration depth of the shaft and the bearing could be expressed as follows:where represents value of eccentricity vector and ; and represent relative displacement of shaft inside bearing in direction and direction; represents clearance value; and , , are radius of bearing and the shaft, respectively.

Then, and represent the position vectors of collision point in the coordinate system:where and represent contact point of bearing and the shaft; is unit vector; and .

The relative velocities in normal and tangential direction can be written aswhere and are speed of contact points of bearing and the shaft; and can be obtained by rotating vector in the counter clockwise direction by 90°.

The joint contact forces are shown in Figure 3. Normal force between the bearing and the shaft can be calculated by using the Lankarani–Nikravesh(L-N) model, which is an easy calculation and has fast convergence [27, 29, 30]:where is stiffness parameter; is hysteresis damping coefficient; and is relative penetration velocity. The stiffness parameter can be given bywhere , are Poisson’s ratios of bearing and shaft; and and represent elastic modulus of bearing and shaft.

Hysteresis damping coefficient could be written aswhere means recovery coefficient; is the initial impact velocity; and if , and , the penetration velocity at moment is .

The modified Coulomb’s friction law presented by Ambrósio, and the dynamic friction force is expressed aswhere is friction coefficient; represents the dynamic correction coefficient; andwhere and are the given bounds for tangential velocity.

Contact force of the revolute pair with revolute clearance could be expressed aswhere and represent the components of contact force of the shaft to the bearing in the and direction, respectively.

3.2. Establishment of Dynamic Model with a Revolute Clearance
3.2.1. Establishment of Kinematics Equation Containing a Revolute Clearance

The coordinate system of the mechanism with a revolute clearance is shown in Figure 4. Local enlargement of the clearance joint A is also shown in Figure 4. For the 2-DOF nine-bar mechanism with clearance, the closed vector equation can be written aswhere and represent relative displacement of shaft inside bearing in direction and direction.

By using the perturbation coordinate approach [25, 26, 31], all the angles of bars and the displacement of slider 9 can be written aswhere and represent the angle of each member without clearance and the disturbance angle, respectively; and and represent the displacement of slider 9 without clearance and the disturbance displacement, respectively.

According to Equations (11) and (12), the following equation could be obtained:where , , , , , , , , , , .

The velocity (Equation (14)) and the acceleration model (Equation (15)) of 2-DOF nine-bar mechanism with clearance can be written as

3.2.2. Establishment of the Dynamic Equation with a Revolute Clearance

Based on the results of kinematics in the section 3.2.1, the dynamic model with clearance is established in this section. According to the Lagrange equation, the dynamic model of nine-bar mechanism with clearance is given bywhere represent kinetic energy, potential energy, and general force, respectively; represents generalized coordinate; and .

Kinetic energy of the mechanism could be expressed aswhere represents mass of the components; represents the velocity of mass center of component ; and represents inertia moment of the component .

Potential energy of the mechanism could be expressed aswhere represents centroid coordinates in the direction.

The generalized force could be expressed aswhere and are external forces and the external torques acting on mass’s center of component , respectively; and and represent moving displacement and rotation angle of mass of component , respectively.

Therefore, the general force could be written aswhere represent driving torque of crank 1 and crank 4, respectively.

According to Equations (16)–(18) and (20), the second-order nonlinear differential equations (Equation (21)) with variable coefficients are both derived:

The expressions of the two driving torques could be obtained by

4. Dynamic Response of Mechanism with a Revolute Clearance

Dynamic equations (Equation (21)) are derived in the previous section, which contain two second-order nonlinear equations. Equations of the system are complicated, and it is difficult to solve equations analytically. Thus, in order to guarantee efficiency and accuracy of the calculation, the fourth-order Runge–Kutta method with variable step sizes programed by MATLAB is utilized, which converts the two second-order differential equations of motion into four first-order equations.

System parameters of the mechanism with clearance are seen in Tables 1 and 2.

It is assumed that driving speeds of two cranks are set as and , the revolute clearance value is 0.5 (mm), and friction coefficient is 0.01, respectively.

Dynamic response of the mechanism is shown in Figure 5. Displacement of slider, velocity of slider, acceleration of slider, the contact force, and the driving torques are, respectively, shown in Figures 5(a)5(f). The red solid line represents the virtual simulation by ADAMS, the blue solid line represents the numerical calculation by MATLAB, and the black dotted line represents ideal state.

From the Figure 5(a), the displacement curve of the slider with a revolute clearance is basically coincided with the displacement curve of slider in ideal state. The maximum value of displacement in numerical calculation appears at 0.3648 s, and the peak value is 0.2718 m. The maximum value of displacement in virtual simulation appears at 0.3638 s, and the peak value is 0.2682 m. The peak difference between ADAMS and MATLAB is 0.0036 m. From the Figure 5(b), it is shown that the clearance has a certain effect on the velocity of the slider. The maximum value of velocity in numerical calculation appears at 0.2731 s, and the peak value is . The maximum value of velocity in virtual simulation appears at 0.2771 s, and the peak value is . The peak difference between ADAMS and MATLAB is . As shown in Figure 5(c), the acceleration fluctuates violently and occurs a great peak, while the mechanism contains a revolute clearance. The peak value of acceleration in numerical calculation appears at 0.5007 s, and its value is . The peak value of acceleration in virtual simulation appears at 0.5016 s, and its value is . The peak difference between ADAMS and MATLAB is . From the Figure 5(d), the peak value of contact force in numerical calculation appears at 0.5007 s, and its value is . The peak value of contact force in virtual simulation appears at 0.5016 s, and its value is . The peak difference between ADAMS and MATLAB is . As shown in Figures 5(e) and 5(f), when mechanism contains a revolute clearance, driving torque is influenced greatly, and there occurs a violent oscillation. According to the numerical calculation results by MATLAB, the peak value of driving torque of crank 1 and crank 4 are both at 0.5007 s, peak values are , and , respectively. Based on the virtual simulation result by ADAMS, the peak value of driving torques of crank 1 and crank 4 are both at 0.5016 s, peak values are , and , respectively. The peak difference of driving torque of crank 1 between ADAMS and MATLAB is . The peak difference of driving torque of crank 4 between ADAMS and MATLAB is . It is shown that the virtual simulation result by ADAMS and numerical calculation result by MATLAB are slightly deviant. Because of the difference of modeling and solving methods and the influence of integral error [32], ADAMS uses classical algorithm; compared with MATLAB programming, the accuracy is higher, and the results of virtual simulation are different from the numerical calculation simulation. Though the comparison curves are slightly different, they have similar regularity and are close to the vibration amplitude and vibration frequency. Therefore, it could be proved that theoretical model is correct.

Based on the numerical calculation result by MATLAB, kinematics characteristic errors of the slider are shown in Figure 6. Displacement error, velocity error, and acceleration error of slider are, respectively, shown in Figures 6(a)6(c). The maximum displacement error of slider occurs in 0.3979 s, and its value is . Maximum velocity error of the slider occurs in 0.4986 s, and its valve is . Maximum acceleration error occurs in 0.5007 s, and its value is .

The penetration depth and shaft center trajectory are shown in Figures 7(a) and 7(b), respectively. It can be seen from the Figures 7(a) and 7(b), the revolute pair changes among the three states: separation, contact, and impact. When penetration depth is negative, the revolute joint is in separation state. While penetration depth is equal to 0, the revolute pair is in contact state, and while penetration depth is positive, revolute joint is in impact state.

5. Chaos Identification of Mechanism Containing a Revolute Clearance

According to the numerical solutions of MATLAB, the phase diagrams and Poincaré portraits of the revolute joint with clearance are plotted, and preliminary qualitative judgement of the chaotic motion of the mechanism is presented. And then, according to the numerical results, the LLE of the mechanism is also solved, and chaotic motion of the mechanism is quantitatively analyzed.

Phase diagrams of 2-DOF nine-bar mechanism are shown in Figure 8. The periodic motion repeats previous motion every other cycle, so its phase diagram is a closed curve. Chaotic motion does not have the periodicity, so its phase diagram is an unclosed curve. According to the Figure 8, the phase diagrams are all in a mess, so mechanism is in chaotic state at the moment.

By observing distribution characteristic of the nonlinear system reflected on the Poincaré cross section, we could judge the motion form of the system. When the Poincaré cross section has only one isolated Poincaré mapping point or very few discrete points, the motion of system is periodic. When a continuous or discrete closed curve appears on the Poincaré cross section, system does quasiperiodic motion. When the Poincaré cross section has some patchy distribution of point set, the system is in a chaos.

Poincaré portraits of 2-DOF nine-bar mechanism are shown in Figure 9. From the Figures 9(a)9(c), there are lots of mapping points on Poincaré cross section, showing irregular distribution, and there is no law between the points. According to characteristics of the Poincaré portrait, it is known that the 2-DOF nine-bar mechanism with clearance is in a chaotic state at this time.

The Lyapunov exponent can effectively describe the sensitivity of initial conditions of the system dynamics, so it can be considered as a criterion for determining strength of chaotic systems. When the LLE is negative, system corresponds to periodic motion. When the LLE is positive, the system corresponds to chaotic motion.

Wolf et al. proposed the Lyapunov exponent, which can be estimated directly based on phase trajectory, phase plane, and phase volume. This kind of method is collectively called the Wolf method, which is widely applied in the study of chaos [5, 33, 34].

The Wolf method is used to calculate LLE, setting time series data for , which correspond to the data of and obtained by solving dynamic equations (Equation (21)), and its corresponding average period, embedding dimension , time delay , and length are all calculated, and then, its corresponding phase space could be express as follows:

Taking initial point as , and then, taking any point near initial point , the distance between the two points is . If the value of is more than the specified value at time point , that is, , is reserved, find a point near point , make , and let the angle between the two points as small as possible, and then, above processes are repeated again and again. When the arrives at the end of time series , total number of iterations during the tracing evolution process is . The LLE could be obtained by the following equation:where is total number of iterations during the tracing evolution process; is the time corresponding to the iteration; is the initial time; and is the distance of the two adjacent points in the phase space. All of these variables are all determined by the time series data.

Lyapunov exponents are shown in Figure 10. From the Figures 10(a) and 10(b), the LLE in direction is 2.7700 and LLE in direction is 1.0742. Because the LLEs in the direction and the direction are greater than zero, so mechanism is in chaos.

The phase diagrams, Poincaré portraits, and LLEs of the mechanism can confirm that the mechanism is in chaos. This is consistent with aperiodic and nonlinear characteristics of dynamic response.

6. The Effect of Main Parameters of the Mechanism on the Chaotic Motions

6.1. Effect of Clearance Value on Chaotic Motion

In order to better study the influence of clearance value on dynamic behavior, the driving speed and friction coefficient should be set as fixed values. It is assumed that driving speeds of crank 1 and crank 4 are and and friction coefficient is 0.01. The effects of different clearance values on chaos phenomena are investigated.

The direction is selected to investigate, and the Poincaré portraits of different clearance values in direction are plotted. When the clearance sizes are 0.01 mm, 0.05 mm, 0.1 mm, and 0.2 mm, corresponding Poincaré portraits are shown in Figures 11(a)11(d). From Figure 11, as the clearance value increases, Poincaré portrait expands, and the chaos phenomenon is more obvious. When the clearance is 0.01 mm, the Poincaré portrait is a point, indicating that mechanism is in periodic motion. When the clearance values are 0.05 mm and 0.01 mm, the Poincaré portrait points get together or form a closed circle, illustrating that mechanism is in quasiperiodic state. When the clearance values is 0.2 mm, the mapping points in the Poincaré portrait are chaotic and irregular, showing that mechanism is in chaos.

When clearance values are 0.01 mm, 0.05 mm, 0.1 mm, and 0.2 mm, corresponding Lyapunov exponents are shown in Figures 12(a)12(d), and the corresponding LLEs are −0.0494, −0.0302, −0.0138, and 0.2580, respectively. When clearance values are 0.01 mm, 0.05 mm, and 0.1 mm, the LLEs are both less than zero, which proves that there is no chaos in the mechanism. While clearance value is 0.2 mm, the LLE is greater than zero, which proves that the mechanism is in chaos. It can be seen that, with increase of clearance size, the LLE increases, and the chaos is more obvious.

Bifurcation is a very important dynamic characteristic of engineering systems. The bifurcation diagram of the system response versus clearance values is drawn to better realize the effect of the clearance value on the system behavior. The clearance values ranged from 0.01 mm to 0.25 mm are selected to research bifurcation phenomena. Bifurcation diagram with varying the clearance size is shown in Figure 13. The axis represents clearance size, and the axis represents relative displacement of shaft inside bearing in direction. From the Figure 13, when clearance size is small, the stability of mechanism is strong, the mechanism performs periodic motions, and with increase of clearance value, mechanism is in chaos, and the chaos phenomenon is more obvious. On the whole, the change of clearance value has great effect on the motion state of the mechanism.

6.2. Effect of Friction Coefficient on Chaotic Motion

In order to better study effect of the friction coefficient on chaos, the driving speed and clearance value should be set as fixed values. It is supposed that the driving speeds of the crank 1 and crank 4 are and , and the clearance value is 1(mm). The effect of different friction coefficients on the nonlinear dynamic behavior has been studied.

The direction of is selected to study, and the Poincaré portraits of different friction coefficients are plotted. When friction coefficients are 0.01, 0.03, 0.06, and 0.1, the corresponding Poincaré portraits are shown in Figures 14(a)14(d). As shown in Figure 14, with the increase of friction coefficient, the chaos phenomenon in the mechanism is weakened.

When friction coefficients are 0.01, 0.03, 0.06, and 0.1, the corresponding Lyapunov exponents are shown in Figures 15(a)15(d), and the corresponding LLEs are 5.4267, 2.8972, 1.2772, and 0.6746, respectively. It can be seen that, with the increase of clearance size, the LLE decreases and chaos phenomenon is more weakened.

Bifurcation diagrams can exhibit the effects of continuous changes in parameters on the stability of the mechanism. The bifurcation diagram with varying the friction coefficient is shown in Figure 16. The axis represents friction coefficient, and the axis represents relative displacement of shaft inside bearing in direction. The variation range of friction coefficient is . From the Figure 16, with the increase of the friction coefficient, the bifurcation diagram tends to converge, chaos is gradually weakened, and the chaotic behavior is dominant in the whole range of the parameter.

6.3. Effect of Driving Speed of Two Cranks on Chaotic Motion

In this section, friction coefficient is set as 0.01, and the clearance value of clearance joint is set as 0.5 mm. The influence of different driving velocities on the nonlinear dynamic behavior of mechanism is studied. The direction is selected to research, and Poincaré portraits are shown in Figures 17(a)17(d). It can be seen from the Figures 17(a) and 17(b), the Poincaré portraits are both a point, so the mechanism is doing periodic motion. The Poincaré portraits of Figures 17(c) and 17(d) are both messy points, it shows that mechanism is in a state of chaos, and the Poincaré portrait of Figure 17(d) is more chaotic than that of the Figure 17(c).

Corresponding Lyapunov exponents are shown in Figures 18(a)18(d) and the corresponding LLEs are −0.0205, −0.0331, 2.0822, and 4.4267, respectively. When driving speed , and , , the LLEs are less than zero, so the mechanism is in the periodic motion state. When driving speed , and , , the LLEs are greater than zero, the mechanism is in the chaotic state, and the LLE of latter is bigger than former, so the chaos phenomena is more obvious when the driving speed is and .

The bifurcation diagram by varying driving speeds is shown in Figure 19. The axis represents driving speed, and the axis represents relative displacement of shaft inside bearing in direction. The variation range of the driving speed of the crank 1 and the crank 4 are and , respectively. When the driving speed of crank 1 ranges from to and the driving speed of crank 4 ranges from to , the mechanism behavior is in period-one state. When driving speed of crank 4 is , the mechanism behavior is in period-two state, since there are only two points on the corresponding bifurcation diagram. When the driving speed of crank 1 ranged from to and the driving speed of crank 4 ranged from to , mechanism shows chaotic motion, and the chaos phenomenon is more and more obvious.

7. Conclusions

This paper studies the nonlinear dynamic behavior of the 2-DOF nine-bar mechanism. The main conclusions are as follows:(i)the nonlinear dynamics model of the 2-DOF nine-bar mechanism with a revolute clearance is established.(ii)Dynamic response of kinematics characteristic of slider, contact force, driving torque, shaft center trajectory, and the penetration depth for 2-DOF nine-bar mechanism are analyzed. It can be seen that clearance has an obvious effect on dynamic response for mechanism, which reduces stability of mechanism and produces vibrations.(iii)Chaos phenomenon existed in the mechanism is identified by the phase diagrams, the Poincaré portraits and LLEs, respectively. The results indicate that the chaos exists in the mechanism.(iv)The influences of the different clearance values, different friction coefficients. and different driving speeds on nonlinear dynamic behavior are studied. Bifurcation diagrams by varying the clearance value, the friction coefficient and the driving speed are drawn. The results show that, with the increase of driving speeds and clearance values, the mechanism gradually changes from periodic state to chaotic state. As the friction coefficient increases, the chaos of the mechanism decreases gradually.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the Natural Science Foundation of Shandong Province (Grant No. ZR2017MEE066), Shandong Young Scientists Award Fund (Grant No. BS2012ZZ008), Taishan Scholarship Project of Shandong Province (No. tshw20130956), and Graduate Innovation Project of Shandong University of Science and Technology (No. SDKDYC180216).