Abstract
This study employs experiments and numerical simulation to analyze the dynamic response of steel beams under huge-mass impact. Results show that lateral torsional buckling (LTB) occurs for a narrow rectangular cross-section steel beam under transverse impact. The experiments were simulated using LS-DYNA. The numerical simulation is in good agreement with experimental results, thus indicating that the LTB phenomenon is the real tendency of steel beams under impact. Meanwhile, the study shows that LS-DYNA can readily predict the LTB of steel beams. A numerical simulation on the dynamic response of H-shaped cross-section steel beams under huge-mass impact is conducted to determine the LTB behavior. The phenomenon of dynamic LTB is illustrated by displacement, strain, and deformation of H-shaped steel beams. Thereafter, a parametric study is conducted to investigate the effects of initial impact velocity and momentum on LTB. The LTB of H-shaped cross-section steel beams under transverse impact is primarily dependent on the level of impact kinetic energy, whereas impact momentum has a minor effect on LTB mode.
1. Introduction
Beams and columns of special steel frame structures should withstand not only normal design load but also explosion, impact, collision, fire, and other extreme loads during the period of service. Hence, the dynamic response and failure behavior of steel members under these extreme loads need to be studied. In fact, explosion and shock have often been associated with fire. Therefore, the effect of the interaction between the structure under impact load and fire needs to be investigated.
Stability is a potent issue in the design of steel structures. Bad stability of steel components may cause serious structural failure. Lateral torsional buckling of a component is a very common type of instability that needs further investigation [1]. Lateral torsional buckling (overall instability) of steel beams under static loads can be described as follows [2]: when steel beams, without lateral support in the net span, are subjected to a certain value of transverse load or moment about the major-axis, a large lateral displacement and twist angle occur as shown in Figure 1. Then, the beam loses its carrying capacity. Obviously, out-of-plane buckling analysis is more difficult than analysis of in-plane bending buckling. However, a number of studies have been conducted based on elastic buckling theory, particularly experimental studies and theoretical analysis [3] on lateral torsional buckling of steel beams under static load. Yang et al. carried out experimental tests and numerical simulations on lateral torsional buckling behavior of singly symmetric I-beams fabricated from Q460GJ steel [4, 5]. Their results showed that steel beams developed lateral torsional buckling under concentrated point loads at the mid-span. Hence, local buckling was not observed. Karmazínová et al. [6] studied the lateral flexural–torsional buckling of steel sigma-cross-section beams with web holes. Their study conclusions became the background of the supplements to specified provisions for the design of steel structures. Kala and Valeš [7] examined a hot-rolled steel I-beam subjected to lateral torsional buckling (LTB) due to bending moment and obtained the stochastic effects of initial imperfection and residual stress on the resistance of the investigated steel beam. The development of numerical methods, such as finite element (FE), has facilitated research on lateral torsional buckling of steel beams from experimental studies to numerical analysis [8–10]. Moreover, research on the fire resistance of steel structures has made great progress [11]. Lateral torsional buckling problem of steel beams under fire has been given considerable attention [12–14]. This problem is due to the substantial number of the torsional buckling modes of steel beam during fire accidents. A steel structure can endure fire for a very long time. Hence, analysis of structures under fire can be regarded as a static problem. When structures are subjected to fire, elevated temperature causes degradation of material properties and triggers plastic deformation of the structure. Thus, elastic buckling theory is not suitable anymore. Therefore, numerical simulation is an effective research method to deal with inelastic buckling problems. For example, Vila Real et al. [15, 16] performed a numerical study and proposed a design curve for the LTB of steel beams. Yin and Wang [17] have numerically investigated the effects of several design factors on lateral torsional buckling bending moment resistance of steel I-beams submitted to nonuniform temperature distribution. Couto et al. [18] studied the lateral torsional buckling of beams with slender cross-sections in the case of fire using a numerical method.
Whether LTB behavior occurs when steel beam is subjected to the combined action of impact and fire loading is difficult to determine because of the complexity of two-load interaction effects. In addition, LTB occurrence on the steel beam when it is only subjected to impact is difficult to investigate because impact load affects the plastic buckling of beams whereas the inertial effect may delay or avoid plastic buckling. On one hand, steel is a rate sensitive material meaning that the yield strength of steel can be improved because of the strain rate effect. Hence, the strain rate effect of steel may lead to the occurrence of elastic LTB. These complex factors hinder the full understanding of LTB behavior of steel beams under impact. In recent years, several studies were conducted on the dynamic response of steel members under impact. The drop weight impact test of clamped steel beams struck transversely at the center by a mass with a rectangular indenter was conducted by Villavicencio and Guedes Soares [19] to define the proper boundary conditions of impacted beams. Wang et al. [20] studied the dynamic behavior of steel beams with fin-plate connections subjected to falling floor impact loads. The results show that steel beams show dynamic behaviors under a certain impact energy with varying impact masses and velocities. A higher impact velocity can slightly improve the energy dissipation rate. Steel beams develop different impact resistances while impact loads are applied at various locations. Al-Thairy and Wang [21] examined the behavior and failure modes of axially compressed steel columns subject to transverse impact by a rigid mass at varying impact speeds and locations. Ning and Zhao [22, 23] completed the experimental study on the stability of aluminum alloy cantilever high beams under cylinder head impact. However, studies on LTB of steel beams under impact load are lacking.
Above all, systemic research on LTB of steel beams under impact load is necessary. This article is a preliminary attempt of related studies. This paper presents the experimental results of the LTB of rectangular section steel beams under drop hammer impact. Afterwards, the LTB behavior of rectangular section and H-shaped steel beams under impact is discussed.
2. Experimental Test Introduction
The drop weight test of steel beams was carried out by setting up a DHR9401 drop hammer test machine. The overall height of the DHR9401 drop hammer test machine is about 13.47 m above the floor and its impact velocity is up to 15.70 m/s. Figure 2 shows the components of the drop hammer test machine. The impact hammer comprises several main components. The first piece is a toughened flat head which is 80 mm long and 30 mm thick and made from high-strength chromium 15 (64HRC) as shown in Figure 3. The second component of the striker is the weight element. This component contains the main mass of the striker and delivers a variety of striker weights from 2 to 250 kg. The third part of the striker is load cells. The overall weight of the impact hammer is 57.8 kg and can be raised to the required height to produce varied impact velocities (up to 15.7 m/s) and energies [24].
Rolled steel plates (Q235) were selected to make the specimens. The specimens were classified according to the thickness of the beam: #A: mm3, #B: mm3, and #C: mm3. The effective spans of steel beams were 500, 750, and 1000 mm, as shown in Figure 4.
Table 1 shows the dimensions of steel beams. The constraints of steel beams were implemented by the length of the reserved. Two supports were fixed on both ends of the experimental platforms. Bolts flank the support side to implement lateral fixation. The rotation constraint was fixed by a steel cover plate as shown in Figure 5. The material properties of the specimens are determined from standard static tensile coupon tests which have been carried out by a CM5105A computer-controlled electronic omnipotence test machine. The experiment was repeated 3 times and the data averaged. Figure 6 shows the stress–strain curve of the steel. The modulus of elasticity and yield stress of the steel beams are 185 GPa and 315 MPa, respectively.
3. Experimental Test Results
3.1. General Results
In these tests, the ultimate deformation modes and lateral and vertical deflections on mid-span of beams were recorded. Then, the influence of section properties, height–width ratio, length–height ratio, and impact velocity on the ultimate deformation modes was studied. Table 2 shows the results of the test specimens, where is the impact height, is the vertical displacement, is the lateral displacement, and is the initial impact velocity.
In Table 2, it is noticed that the lateral and vertical deflections increase with the increase in impact velocity. It is also observed that the deformation of beams is influenced by the height–width and length–height ratio. The length–height ratio varied at a range of 10 and 20 of specimens A1-3 and A3-1, while the impact velocity is 3.13 m/s and height–width ratio is 10. The lateral displacements of specimens A1-3 and A3-1 are 0.9 and 21.78 mm, respectively. The vertical displacements of specimens A1-3 and A3-1 are 9.82 and 20.4 mm. Thus, the deformation of the beams is proportional to the length–height ratio. The height–width ratio varies at a range of 10, 14, and 20 of specimens A3-3, B3-2, and C3-2 whereas the impact velocity is 5.42 m/s and length–height ratio is 20. The lateral displacements of specimens A3-3, B3-2, and C3-2 are 75.1, 51.3, and 17.6 mm, respectively. The vertical displacements of specimens A3-3, B3-2, and C3-2 are 94.4, 46.36, and 18.14 mm, respectively. Thus, the deformation of the beams is inversely proportional to the height–width ratio. Moreover, the ultimate deformation modes of beams are influenced by section properties, height–width and length–height ratios, and impact velocity.
3.2. Ultimate Deformation Modes
Figure 7 shows the ultimate deformation modes of specimen #A1 for various velocities. Specimen A1-1 has noticeable in-plane deformation without the lateral displacements on twists. The local deformation of specimens A1-2 and A1-3 has been found without the LTB. When impact velocity reaches 4.43 m/s, in specimens A1-4 and A1-5, the LTB occurs. At the same time, specimen A1-5 has an risen deformation. The experimental results show that LTB occurs for rectangular section steel beams under transverse impact. This LTB phenomenon is mainly influenced by the impact velocity. It can be seen that the LTB critical velocity of A1 must be between 3.13 and 4.43 m/s.
Figure 8 shows the ultimate deformation modes of specimen #A2 with different velocities of 4.43, 5.42, and 6.26 m/s. Specimens A2-1, A2-2, and A2-3 all show LTB behavior. In the meantime, the deformation of #A2 increases with the increase of the impact velocity. Figure 9 shows the ultimate deformation modes of specimen #A3. Specimen #A3 also has the LTB behavior. The LTB deformation of #A3-1 is mainly the out-of-plane bending deformation. The torsion deformation of #A3-1 is relatively small. From the ultimate deformation modes of A1-4, A2-1, and A3-2, it can be observed that the length–height ratio has a considerable influence on LTB modes.
Figure 10 shows the ultimate deformation modes of specimen #B1 for various velocities. When the impact velocity is 3.13 m/s, specimen B1-1 has an in-plane deformation without the lateral displacements on twists. Specimen B1-2 has an out-plane deformation with an impact velocity of 4.43 m/s. When the impact velocity reaches 5.42 and 6.26 m/s, specimens B1-3 and B1-4 have the behavior of LTB with flexural buckling as the main deformation. When the impact velocity reaches 7 and 7.67 m/s, specimens B1-5 and B1-6, respectively, have different LTB modes with torsional deformation. Once again, it can be observed that the impact velocity has a considerable influence on the deformation of beams. Figure 11 shows the ultimate deformation modes of specimen #B3 with different velocities. It can be seen that specimen B3-1 presents the flexural–torsional buckling behavior with flexural buckling as the main deformation. When the impact velocity reaches 5.42 and 6.26 m/s, the beam presents various LTB modes.
Figures 12, 13, and 14 show the ultimate deformation modes of specimens #C1, #C2, and #C3 with different velocities. Flexural–torsional buckling noticeably occurs for specimens #C1, #C2, and #C3 at different impact velocities. At the same impact kinetic energy, specimens B3 and C3 have different LTB modes. This difference is caused by the variation of section heights between B3 and C3. The LTB of the steel beams with smaller section height is similar to the first-order static buckling mode. With the increase of the section height, the torsional stiffness of the steel beam increases and the possibility of first-order static buckling-like LTB mode decreases. As a result, the steel beams with larger section height are more likely to behave the second-order buckling mode.
The abovementioned results show that rectangular section steel beams under transverse impact can lead to LTB behavior. The LTB mode of rectangular section steel beams under impact is not a single pattern and depends on the geometry of the specimen.
4. Numerical Simulation
4.1. Finite Element Model
Numerical simulation is carried out by finite element code LS-DYNA and element SOLID164 was used to simulate the solid beam and drop hammer. SOLID164 [25] is used for the 3D modeling of solid structures. The element is defined by eight nodes having the following degrees of freedom at each node: translations, velocities, and accelerations in the nodal , , and directions. This element is used in explicit dynamic analyses. In this paper, the elastic-perfectly plastic model is adopted for the numerical analysis using the Material Type 3 [26] of LS-DYNA. This model is suited for modeling isotropic and kinematic hardening plasticity with the option of including rate effects. As a type of strain rate sensitive material, steel has a considerably high yield stress under high strain rates. When structures are subjected to impact loads, strain rate effect has a considerable influence on the dynamic behavior and damage [20]. The strain rate effects are described by employing the Cowper–Symonds [27] equation with material parameters of , . The relationship between yield stress and strain rate is expressed in the following equation:where and are the dynamic and static yield stress of the material, respectively, and is the equivalent plastic strain rate. The values for and have been chosen to be consistent with those given in the previous section. The density of steel is 7,850 kg/m3 while the Poisson ratio is 0.3. The values of elastic modulus, yield strength, and tangent modulus are similar to those in Figure 3 which features the stress–strain curve. Since the current study aims to investigate the behavior of steel beams under transverse impact load, the emphasis is on the beam rather than the drop hammer. Therefore, only a brief description is given on the drop hammer which is assumed to be a rigid cuboid with dimensions of mm3. To increase the calculation precision and efficiency, the impact location and nearby supports are simulated in a fine mesh whereas a coarse mesh is used elsewhere. The drop hammer is modeled as a nondeformable solid object. The number of elements has no effect on the calculation precision. Hence, the drop hammer is simulated in a coarse mesh and divided into 1,536 eight-node hexahedrons. A steel beam is divided into 8,216 eight-node hexahedrons. The maximum mesh size of the solid elements is 15 mm as shown in Figure 15. An initial imperfection is generally applied to a structure model in buckling analysis. The buckling mode under impact load should be the superposition of static buckling modes whereas the percentage of the low order buckling mode is comparatively high. For the low-speed impact problem in this paper, buckling should be close to the first-order mode. Therefore, in the analysis, the initial imperfection has the same shape as the first-order mode for bending vibration with amplitude of , where is the beam span.
In this paper, the boundary conditions of beam are simulated by restraining the , , and direction node displacement in both ends of the steel beam constrained area. The impact hammer moves only in the vertical direction while other displacement and angle directions are restrained. In the numerical simulation analysis, the impact energies are given by setting the different initial velocities of the drop hammer. During the impact experiment, the contact interface is formed using the impact hammer and steel surface. When two solid bodies touch, contact stresses are transmitted across their common interface. This contact interface is able to pass the force but is not allowed to penetrate. To model the contact behavior between the front surface of the impact hammer and top surface of beams, the automatic surface-to-surface constraint algorithm with friction is used, where the friction coefficients are introduced to simulate the friction forces transmitted across the contact interface.
4.2. Comparison of Simulation and Experiment Results
The numerical simulation on the experiment of steel beams #A1, #B3, and #C2 under the drop hammer was carried out using the previously mentioned FE technique. Then, the simulation and experiment results are compared. Figure 16 shows both results of deformation mode of #A1. It is noticed that, under different impact velocities, the deformation behavior of the steel beams is different. When m/s, the experimental and numerical results show that the beams do not produce LTB, but only the in-plane bending deformation. When m/s, local deformation occurs in the mid-span and the steel beams do not produce overall buckling yet. However, at the initial impact velocity of 4.43 m/s, the experiment and numerical results all show that the beam produces out-of-plane buckling and mainly lateral bending deformation. When impact velocity is up to 6.26 m/s, the LTB behavior of beams is presented by experimental and numerical results and mainly torsional deformation. Figures 17 and 18 show the experimental and numerical results of the deformation mode of #B3 and #C2, respectively. It is noticed that the numerical simulation results are basically consistent with the experiment results under varied impact velocities. Table 3 shows the lateral and vertical deflections at mid-span of beams under different impact velocities. It is noticed that the results of the numerical simulation are in good agreement with the experimental results. Thus, the validity and reliability of the simulation method used in the paper are verified. At the same time, it is observed that the numerical simulation method can predict the LTB of steel beams.
5. Numerical Simulation on the LTB of H-Shaped Steel Beams under Transverse Impact
Nowadays, H-shaped steel beams are commonly used as important bearing component of steel structures. Steel beams and columns may fail under a lateral torsional buckling mode under impact and lead to the collapse of the structure. Therefore, studying the buckling tendency of H-shaped steel members under impact is increasingly important [28]. In this section, the numerical simulation analysis on the dynamic response of H-shaped steel beams is carried out using the numerical simulation method which has been verified. The LTB behavior of H-shaped steel beams is studied based on the characteristics of lateral displacement and strain. In this research, the following important parameters have been identified for parametric study: impact velocity and impact momentum.
5.1. Numerical Model
Figure 19 shows the dimension of an H-shaped cross-section; the span of the beam is 1.5 m. The beam is fixed on both ends. Steel is assumed to be grade Q235 while the steel modulus of elasticity and yield strength are assumed to be 206 GPa and 235 MPa, respectively. The impact energies are given by setting different impact heights. The mass of the drop hammer is assumed to be 200 kg in this research.
Explicit dynamical analysis is used to simulate the dynamic response of the beams under impact. The simulation model comprises the following mesh and material properties.
Strain rate effect: Cowper–Symonds equation with s−1, , and mesh: steel beam is divided into 22,716 eight-node hexahedrons (SOLID164) and the maximum mesh size of the solid elements is 15 mm as shown in Figure 20. Since the current study aims to investigate the behavior and failure modes of H-shaped steel beams under transverse impact load, the emphasis is on the beam rather than the drop hammer. Hence, only one drop hammer geometry with a dimension of mm3 is considered. In the numerical simulation analysis, the density of the drop hammer is adjusted to obtain the required mass. The drop hammer is divided into 1,536 eight-node hexahedrons. The initial distance between the striking tip of the drop hammer and top surface of H-shaped steel beam is set at 5 mm. In this section, the initial imperfection has the same shape as the first mode for bending vibration with an amplitude of , where is the beam span.
5.2. Numerical Results
In the current section, the numerical results of the H-shaped steel beam under transverse impact, such as the entire process of the impact, residual deformation of beam, and time history curves of strain, will be discussed. Output data curves of displacement during the simulation are at a 3,750 Hz sampling rate.
5.2.1. Plastic Deformation of Beam
Figure 21 shows the residual deformation of the H-shaped steel beam under transverse impact when the initial impact velocity is 6 m/s. It is noticed from Figure 21 that an LTB behavior is observed while the maximum of the lateral deformation is produced mid-span. At the same time, a local deformation of the beam flange is observed at the location of the steel beam top flange in contact with the impact hammer.
5.2.2. Impact Process
Figure 22 illustrates the whole dynamic response of H-shaped steel beams under transverse impact. It is observed that the lateral displacement of the beam increases with the increase in contact time of the drop hammer and beam. The maximum lateral displacement occurs at an impact time of 13 ms, which includes elastic and plastic deformation. With the movement of the drop hammer and beam, the elastic deformation of beam gradually recovers and the whole displacement of the beam is reduced until the elastic deformation of beam fully recovers. The plastic deformation of the H-shaped steel beam occurs at 35 ms.
5.2.3. Plastic Deformation of Beam
In the current section, the lateral and vertical displacements–time history curves at different nodes of the steel beam are studied. Figure 23 shows the location of the nodes and all points belong to the beam web.
The lateral displacement–time history curves at different nodes of the steel beam are given as shown in Figure 24. It is observed that the lateral displacement of nodes from various locations of beam is different. The residual lateral displacements of nodes A and E are 15.4 mm whereas node C is at 36.7 mm. The difference of the residual lateral displacements on nodes A and C indicates that the beam produces lateral bending deformation. The residual lateral displacements of nodes A, C, and E are greater than zero and corresponding displacements of nodes B, D, and F are less than zero. This shows that the out-of-plane torsion deformation of the steel beam is produced during movement. From the abovementioned discussion, the H-shaped steel beam under transverse impact produces the LTB phenomenon.
Figure 25 shows that the residual vertical displacements of nodes C and D are 41.2 and 28.3 mm, respectively. From the difference of the residual vertical displacement, it can be seen that the H-shaped steel beam under transverse impact occurs in not only in-plane bending deformation but torsional deformation. Figure 25 also demonstrates that the residual vertical displacements of nodes A and B agree with nodes E and F, respectively. This agreement is brought by the constraint and symmetrical load of specimen. Hence, the residual vertical displacement–time history curves of the beam are also symmetrical.
5.2.4. Strain–Time History Curves
Figure 27 presents the strain–time history curves for the two groups of symmetrical elements from the mid-span of the H-shaped steel beam web and flange, respectively. Figure 26 shows the locations of symmetrical elements. Figure 27 shows that the strain of the two symmetrical elements of the beam web is coincident during 0 to 4 ms. This behavior indicates that the H-shaped steel beam web moves in a plane at the beginning and then produces an out-of-plane movement. Figure 27 also shows that the strain–time curves of the two symmetrical elements on the beam flange are different at the beginning of the movement. This result indicates that the H-shaped steel beam flange produces out-of-plane movement first. This is because the drop hammer comes into contact with the beam flange directly and leads the top flange to produce the buckling first.
Based on the above discussion, it can be found that the H-shaped steel beams under transverse impact can cause the LTB behavior.
5.3. Effects of Different Drop Hammers
In this research, the following two important parameters have been identified for parametric study: impact velocity and impact momentum. This section will present numerical simulation results to help understand the effects of these design parameters.
5.3.1. Effects of Impact Velocity
In this section, the dynamic responses of H-shaped steel beams are studied at different impact velocities where the initial impact velocity is 2, 4, 6, and 8 m/s.
Figure 28 presents the plastic deformation of the H-shaped steel beam under different impact velocities. When m/s, the H-shaped steel beam only occurs in in-plane bending deformation whereas LTB behavior is not observed. When = 4, 6, and 8 m/s, it is observed that the H-shaped steel beam produces the LTB accompanied by the local deformation of the beam flange. Figure 29 shows the impact response for the lateral displacement of mid-span of the H-shaped steel beam under different impact velocities. It can be seen that the plastic deformation of H-shaped steel beam increases with the increase of initial impact velocity. Hence, with the increase in initial impact velocity, the H-shaped steel beams are more sensitive to the LTB.
5.3.2. Impact Momentum Effects on the LTB of H-Shaped Steel Beams
In this section, the dynamic response analyses on H-shaped steel beam subjected to transverse impact are carried out using three cases. The three cases have different masses of drop hammer, impact velocities, and same kinetic energy of 3600 J. Table 4 shows the three sets of conditions and impact momentum effects on the LTB of H-shaped steel beams, where is the impact kinetic energy of the drop hammer, is the impact momentum of the drop hammer, and is the contact time.
Table 4 shows that if the impact kinetic energy remains constant, then the residual vertical deformation of the H-shaped steel beam and contact time decrease with the decrease in impact momentum. The transfer energy is proportional to the mass ratio whereas the mass ratio decreases with the decrease in impact momentum. Thus, the residual vertical deformation of the H-shaped steel beam is small under low impact momentum. Figure 30 also shows that the change of the impact momentum has negligible effect in the LTB forms of the H-shaped steel beam. Figure 31 presents the time history curves of the lateral displacement on the beam mid-span. It can be seen that the lateral displacement decreases with the decrease in impact momentum. Under constant impact kinetic energy and varying impact momentum, the maximum vertical displacement difference is only 1.6 mm but the maximum lateral displacement difference is 7.7 mm as shown in Table 4 and Figure 31. This is because the LTB deformation is governed by the lateral displacement compared to vertical displacements. Based on the discussion above, it is conducted that the impact momentum is not a deciding factor for the LTB of the H-shaped steel beam under transverse impact.
6. Conclusion
Through drop weight impact test, it is found that the narrow rectangular cross-section steel beam under transverse impact leads to LTB. The numerical simulation results are consistent with experimental results, which means that the LTB behavior from the experimental test is not an accidental consequence due to structure or material defects, but the real response mode of the structure. The consistency of the numerical simulation and experiment results also shows that the FE code ANSYS/LS-DYNA can effectively simulate LTB with material and geometrical nonlinearities. On this basis, the abovementioned numerical method is used to conduct a study on LTB of H-shaped steel beams subject to transverse impact. From the results, the following conclusions are drawn:(1)After studying the impact process, the deformation modes, lateral displacement, and strain of the H-shaped steel beam, it can be seen that not only can the rectangular section beam under transverse impact lead to LTB, but the H-shaped steel beam under transverse impact also easily causes LTB accompanied by the local deformation of the top flange.(2)The LTB of the H-shaped steel beam under transverse impact is primarily dependent on the level of impact kinetic energy. With the increase in impact kinetic energy, the H-shaped steel beam subject to transverse impact is more prone to LTB and local deformation.(3)If the impact kinetic energy remains constant, then the plastic deformation of the H-shaped steel beam and contact time increase with the increase in impact momentum whereas different impact momentum has a minor effect on LTB of the H-shaped steel beam under transverse impact.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (Grant nos. 10872117, 11272189, 11672165, and 51208289).