Application of strength reduction method to dynamic anti-sliding stability analysis of high gravity dam with complex dam foundation

Deng-hong CHEN*, Cheng-bin DU

Department of Engineering Mechanics, College of Mechanics and Materials, Hohai University, Nanjing 210098, P. R. China

Application of strength reduction method to dynamic anti-sliding stability analysis of high gravity dam with complex dam foundation

Deng-hong CHEN*, Cheng-bin DU

Department of Engineering Mechanics, College of Mechanics and Materials, Hohai University, Nanjing 210098, P. R. China

Considering that there are some limitations in analyzing the anti-sliding seismic stability of dam-foundation systems with the traditional pseudo-static method and response spectrum method, the dynamic strength reduction method was used to study the deep anti-sliding stability of a high gravity dam with a complex dam foundation in response to strong earthquake-induced ground action. Based on static anti-sliding stability analysis of the dam foundation undertaken by decreasing the shear strength parameters of the rock mass in equal proportion, the seismic time history analysis was carried out. The proposed instability criterion for the dynamic strength reduction method was that the peak values of dynamic displacements and plastic strain energy change suddenly with the increase of the strength reduction factor. The elasto-plastic behavior of the dam foundation was idealized using the Drucker-Prager yield criterion based on the associated flow rule assumption. The result of elasto-plastic time history analysis of an overflow dam monolith based on the dynamic strength reduction method was compared with that of the dynamic linear elastic analysis, and the reliability of elasto-plastic time history analysis was confirmed. The results also show that the safety factors of the dam-foundation system in the static and dynamic cases are 3.25 and 3.0, respectively, and that the F2 fault has a significant influence on the anti-sliding stability of the high gravity dam. It is also concluded that the proposed instability criterion for the dynamic strength reduction method is feasible.

dynamic anti-sliding stability; complex dam foundation; dynamic strength reduction method; instability criteria; elasto-plastic model; dynamic time history analysis; gravity dam

1 Introduction

Complex dam foundations, containing various rock formations, faults, and weak structural layers in the rock mass of the dam foundation, present a challenge to gravity dam design, in which the shallow and deep anti-sliding stability is the main concern. The anti-sliding stability is generally expressed through a safety factor calculated by the rigid body limit equilibrium method according toDesign Specification for Concrete Gravity Dams(SETC 2000), which has significant limitations in practical use. At present, various numerical methods are commonly used in studying the static anti-sliding stability of complex dam foundations, including the strength reduction method (Ge et al. 2003; Cheng et al. 2007; Zhou et al. 2008), the overload method (Yu and Ren 2007; Li et al. 2009), and the rigid body-spring element method (Zhang et al. 2001).

As related to the dynamic parameters of materials, seismic input methods, failure criteria, and other complex problems, nonlinear dynamic analysis of the complex dam-foundation system has been conducted by many scholars. Dai and Li (2007) proposed the strength reduction method for dynamic stability analysis of a complex rock slope under seismic loads with the FLAC software. Wang et al. (2009) analyzed the dynamic strength reduction stability of a plane finite element model of a gravity dam with a two-slip surface using the ADINA software. To investigate landslides induced by the Wenchuan earthquake in China, Zheng et al. (2009b) and Ye et al. (2010) conducted numerical analyses of seismic slope failure mechanisms with the dynamic strength reduction method using the FLAC3D software, in which the tensile shear failure was considered.

In this study, elasto-plastic dynamic analysis of a high roller-compacted concrete (RCC) gravity dam with a complex foundation was conducted and the dynamic strength reduction method was employed to analyze the anti-sliding stability of the dam foundation under dynamic loads. A corresponding dynamic instability criterion was also proposed in this paper.

2 Calculation methods

2.1 Incremental dynamic equation

In dynamic analysis with the nonlinear finite element method, the incremental dynamic equilibrium equation at timetis

The implicit direct integration method of the Hilber-Hughes-Taylor recursive form based on the Newmark method was used to solve Eq. (1) in this study.

2.2 Dynamic strength reduction method

The strength reduction method is a limit equilibrium analysis method based on the definition of the safety factor. It emphasizes the uncertainty and possible weakening effect of the material strength. The shear strength parameters of foundation materials, affected by many factors, are quite difficult to define accurately, and the actual strength is probably lesser than the standard strength that the design requires. In the strength reduction method, the cohesioncand the angle of internal frictionφof materials are decreased continuously, and expressed as

when the structure reaches its limit bearing capacity, whereKis the strength reduction factor of materials at the limit bearing capacity state, which is defined as the safety factor of the dam-foundation system.

The procedure of the dynamic strength reduction method is as follows (Dai and Li 2007): First, the static analysis is carried out by decreasing the shear strength parameters of the dam foundation’s rock mass in equal proportion. Then, the seismic time history analysis is carried out. When the dam foundation reaches the limit bearing capacity under static loads, the dynamic analysis is not processed. Considering the short duration of earthquakes, we assumed the dynamic shear strength parameters of the rock mass to be equal to the static shear strength parameters in this paper. According toSpecifications for Seismic Design of Hydraulic Structures(SETC 2001), the dynamic elastic modulus of concrete is 1.3 times the static elastic modulus.

2.3 Instability criteria

Recently, three major static criteria have been proposed for determining when the dam-foundation system reaches an unstable state based on the strength reduction method (Zhou et al. 2008; Wang et al. 2009): (1) the calculation non-convergence criterion; (2) the displacement mutation criterion; and (3) the plastic yield zone connection criterion, in which the dam foundation is considered unstable when the plastic zones are connected.

Usually, all or some criteria mentioned above are taken as the instability criteria under static loads. When the dynamic strength reduction method is used, the dynamic response of key parts of the dam changes with the decrease of the shear strength parameters of the material. In this paper, the sudden change of the peak values of the dynamic displacements and plastic strain energy with the increase of the strength reduction factor is proposed as an instability criterion for the dynamic strength reduction method. This dynamic criterion was constructed by combining the mutation criterion of dynamic displacements and plastic strain energy and the static criteria (1) and (3).

3 Case study

The RCC gravity dam, with a maximum height of 200.5 m, is located on the middle Beipan River in Guizhou Province in China. The crest length of the dam is 410.0 m, the crest width of the non-overflow section is 12.0 m, the platform width of the overflow dam section is 35.2 m, and the largest bottom width of the dam is 159.05 m. The dam is divided into 20 dam monoliths. Its foundation is situated in the rock layer with complex geological structures. The main geological structure comprises limestone or argillaceous limestone of the Yongningzhen formations T1yn1-1, T1yn1-2, and T1yn1-3, with T1yn1-1occupying the most space. Two faults of the dam, the F1 and F2 faults, extend through the dam foundation to the two banks, with theF1 fault extending through the upstream foundation and mainly developing in T1yn1-1. The F1 fault strike is 20° to 30° with the dam axis and tends toward downstream, the inclination is about 80°, and the effect width is about 30 m. The F2 fault extends through the foundation, about 130 m downstream of the F1 fault, the inclination is about 72°, and the effect width is about 10 m. The two faults occur near the critical site of the dam heel and dam toe and significantly influence the deformation, stress, and anti-sliding stability of the dam.

3.1 Numerical model

A three-dimensional finite element numerical model of an overflow dam section (No. 10 dam monolith with a width of 20.5 m) was established, and the finite element discretization of the dam-foundation system is shown in Fig. 1. Various types of the rock mass were considered in the finite element model, and the F1 and F2 faults were also considered. The foundation meshes were extended a distance of approximately 1.5 times the dam height toward the upstream, downstream, and downward directions. In the global coordinate system of the model, they-axis is horizontal along the positive direction from upstream to downstream, thez-axis is in the vertical direction, and thex-axis is along the dam axis. The dam foundation was discretized into eight-node brick elements with 15 024 elements and 17 839 nodes. The Gauss integration method with 2 ×2×2 Gauss points was adopted to assemble the stiffness matrix.

Fig. 1 Finite element model of overflow dam section

In evaluating the earthquake performance of the dam-foundation system, it is extremely important to establish a reasonable earthquake input model. The earthquake input model most commonly used at present is the massless foundation model proposed by Clough et al. (1985). However, many studies have shown that the influence of an infinite foundation should be considered in numerical analysis of elastic wave propagation. So far, the artificial transmitting boundary (Liao et al. 1984), Higdon boundary (Higdon 1986), and viscous-spring boundary (Lysmer and Kuhlemeyer 1969; Liu et al. 2006; Du et al. 2006) methods have been employed to study the energy dissipation effect of an infinite foundation. The boundary methods mentioned above and the corresponding earthquake input mechanism were deduced based onthe homogenous media and linear elastic wave theory. However, the rock mass examined in this study was regarded as non-homogenous media, so there were still some limitations when the above methods were applied in the practical projects. We think the fixed boundary condition is feasible from an engineering point of view, i.e., the horizontal displacement on the lateral boundary is constrained, and all the displacements on the bottom boundary are constrained completely.

3.2 Basic parameters of materials

The physical-mechanical parameters of concrete and rock masses are shown in Table 1. As the mass and stiffness proportional damping was applied in the dam-foundation system, the first and second vibration modes were considered. The fundamental frequency and the coefficient of the stiffness matrix of the model are 0.435 4 Hz and 0.005 2, respectively, when the damping ratio of concrete is assumed to be 5%.

Table 1 Physical-mechanical parameters of concrete and rock masses

3.3 Load analysis

The static loads under normal operating conditions included the self-weight, hydrostatic pressure acting on the upstream and downstream faces of the dam, silt pressure, and uplift pressure. The dynamic loads included earthquake loads in the horizontal and vertical directions and the hydrodynamic pressure. The acceleration time histories of the design earthquake produced from the design response spectra are shown in Fig. 2. The peak earthquake accelerations are 0.123gin the horizontal direction and 0.082gin the vertical direction, respectively, wheregis the gravitational acceleration. The earthquake load was input by exerting a uniform inertia force on the dam body (Zhang et al. 2009). The effects of hydrodynamic pressure acting on the upstream face of the dam were approximated with theadded mass assumption (Westergaard 1933):

wherePwis the hydrodynamic pressure acting on the upstream face of the dam,ahis the design earthquake acceleration in the horizontal direction,ρwis the density of water,His the water depth, andhis the height above the base of the dam. Although Eq. (3) gives a conservative estimate of hydrodynamic pressures (Saleh and Madabhushi 2010), it is appropriate from an engineering point of view.

Fig. 2 Acceleration time histories of design earthquake

All static and dynamic loads mentioned above were applied according to the following steps: (1) The self-weight of the foundation was first applied, and the self-weight stress of the foundation was taken as the initial stress of the system. (2) The self-weight of the dam was applied. (3) The reservoir water pressure acting on the upstream and downstream dam faces, silt pressure, and uplift pressure were applied. (4) The hydrodynamic pressure and the earthquake load were applied at a fixed time step Δt= 0.02 s.

3.4 Constitutive model and yield criterion of rock mass

The isotropic, elasto-plastic model with linear softening characteristics was used to examine the constitutive relationship of the rock mass materials of the dam foundation. The most typical Drucker-Prager (D-P) criterion in geotechnical engineering was used as the yield criterion:

whereFis the yield function, andI1andJ2are the first invariant of the stress tensor and the second invariant of the partial stress tensor, respectively. Bothαandkare positive constants that depend oncandφof the materials, andcandφdepend on the correlation between the Mises conical surface and the Mohr-Coulomb hex-pyramid surface. Different yield criteria can be realized in finite element calculations by modifying the expressions ofαandk. In this study,αandkwere given as follows (Zheng et al. 2009a):

4 Result analysis

4.1 Static strength reduction analysis

The horizontal displacement and vertical displacement at the dam crest and the plastic strain energy of the dam-foundation system varying with the strength reduction factorKsare shown in Figs. 3 and 4. A coefficientr, which represents the change rate of the displacements or plastic strain energy, is defined as follows:

where Δuis the displacement or plastic strain energy increment corresponding to the increment of the strength reduction factorΔKs.

Fig. 3 Relationship between static displacements at dam crest and strength reduction factorKs

Fig. 4 Relationship between plastic strain energyEsand strength reduction factorKs

Figs. 3 and 4 and Eq. (6) show that the change rates of the static displacements and the plastic strain energy are –12.44 cm, –8.56 cm, and 289.92 MJ, respectively, when thestrength reduction factor is increased toKs=3.25, which are greater than those forKs＜3.25. In other words, the displacements and plastic strain energy of the dam-foundation system change dramatically whenKsof the rock mass is increased to 3.25. The displacements and plastic strain energy have a similar mutation pattern as the strength reduction factor changes.

Fig. 3(a) shows that the displacements of the dam are directly related to the shear strength parameters of the rock masses of the dam foundation. The foundation rock masses of the project (Fig. 1), from upstream to downstream, are the T1f2-1, T1f2-2, T1f2-3, and T1yn1-1layers; the F1 fault; the T1yn1-2layer; the F2 fault; and the T1yn1-3, T1yn2, and T1yn3-1layers. The rock masses from the dam heel to the upstream truncated boundary are T1f2-3, T1f2-2, and T1f2-1, respectively, and the strengths of the three types of rock masses are generally low (Table 1), especially the T1f2-2rock mass. With the increase of the strength reduction factor, the horizontal displacement components produced by the dam weight gradually offset the horizontal displacements generated by the hydrostatic pressure, resulting in dam deformation from downstream to upstream, which is different from other projects (Zhou et al. 2008; Wang et al. 2009), but the mutation patterns of the displacements and plastic strain energy are significant.

The plastic yield zones of the dam foundation under the static loads atKs=2.0 andKs=3.0are shown in Fig. 5. From the progressive failure process, it can be seen that, with the decrease of the shear strength of the rock mass, the tension shear failure occurs first in the dam heel region, and the failure area is very small. Then, the compression shear failure occurs in the dam toe region, and the plastic failure area gradually enlarges, occurring in the F2 fault, the T1f2-3and T1yn1-2rock layers, and the F1 fault, successively. WhenKsis increased to 3.25, the plastic yield zones from upstream to downstream are connected completely, and the maximum equivalent plastic strain value is 4.685 × 10-3. WhenKsis increased to 3.45, the nonlinear computation has not yet converged, and the dam-foundation system nearly achieves its limit bearing capacity.

Fig. 5 Plastic yield zones under static loads

According to the mutation criterion of static displacements and the plastic strain energy and the plastic yield zone connection criteria, the safety factor of the dam-foundation system is 3.25 under static loads, which means that the anti-sliding stability is guaranteed whenKs≤3.25.

4.2 Dynamic strength reduction analysis

The elasto-plastic time history analysis of the dam was carried out according to the dynamic computational method presented in section 2.1. In order to confirm the accuracy of the elasto-plastic time history analysis, the results of dynamic linear elastic analysis and elasto-plastic analysis (Ks=1.0) were compared, as shown in Table 2 and Fig. 6. The maximum horizontal displacementUymax, maximum vertical displacementUzmax, maximum principal tensile stress at the dam heelσ1max, and maximum principal compressive stress at the dam toeσ3maxobtained by the two analysis methods are very close. Because the nonlinear characteristic of the rock mass of the dam foundation was considered, residual deformation was produced, and a slight drift in the crest displacements of the dam was observed in elasto-plastic analysis, which led to a little difference between the calculation results of the two analysis methods. Therefore, the elasto-plastic time history analysis method presented in this paper is feasible.

Table 2 Comparison of liner elastic and elasto-plastic time history analyses

Fig. 6 Results of liner elastic and elasto-plastic time history analyses

The peak values of the horizontal displacement and vertical displacement at the dam crest and plastic strain energy of the dam-foundation system varying with the strength reduction factor are shown in Figs. 7 and 8, respectively. It is shown that the displacements and plastic strain energy of the dam-foundation system change suddenly whenKsis increased to 3.0. Thechange rates of the peak values of the dynamic displacements and plastic strain energy are–7.51 cm, –6.54 cm, and 435.32 MJ, respectively, atKs= 3.0, which are greater than those forKs< 3.0. The change patterns of the displacements and plastic strain energy with the strength reduction factor are basically the same.

Fig. 7 Relationship between peak values of dynamic displacements at dam crest andKs

The horizontal displacement at the dam crest under static loads varying with the strength reduction factor is shown in Fig. 3(a). The deformation range of the dam crest in the horizontal direction under earthquake loads is –3 to 3 cm (Ren et al. 2009). Under the combined action of the static and dynamic loads, the horizontal hydrodynamic pressure and earthquake loads are added, so the maximum horizontal displacement at the dam crest is towards downstream, and it also has a descending trend. The time history of the horizontal displacement at the dam crest with the strength reduction factor is shown in Fig. 9.

Fig. 8 Relationship between peak value of plastic strain energyEsmaxof dam-foundation system andKs

Fig. 9 Time history of horizontal displacement at dam crest at differentKs

The plastic yield zones of the dam foundation under the dynamic loads atKs=2.0 andKs=3.0are shown in Fig. 10. It also can be seen from the progressive failure process that,with the decrease of the shear strength of the rock mass, tension shear failure occurs first in the dam heel region, and the failure area is very small. When the strength reduction factor is increased toKs=3.0, the plastic yield zones from upstream to downstream are connected completely, and the residual deformation appears at the same time, which is reflected in Fig. 9. When the strength reduction factor is increased toKs=3.35, the nonlinear computation does not converge in the dynamic steps.

Fig. 10 Plastic yield zones under dynamic loads

According to the dynamic displacement and plastic strain energy mutation criterion and the plastic yield zone connection criterion, the safety factor of the dam-foundation system is 3.0 under the dynamic loads.

Figs. 5 and 10 also show that the F2 fault almost reaches the maximum value of the equivalent plastic strain, which means that it has a significant influence on the anti-sliding stability of the overflow dam monolith. The F1 fault has a limited influence because of the cutoff trench near the dam heel. Engineering measures could be taken to deal with the F2 fault to further improve the safety factor of the dam-foundation system.

5 Conclusions

(1) Through the static and dynamic strength reduction analysis and dynamic linear elastic analysis of an overflow dam monolith of a high gravity dam, the reliability of elasto-plastic time history analysis based on the dynamic strength reduction method is confirmed. The results show that the safety factors of the dam-foundation system in the static and dynamic cases are 3.25 and 3.0, respectively, and that the F2 fault has a significant influence on the anti-sliding stability of the high gravity dam.

(2) The instability criterion of the dynamic strength reduction method proposed in this paper is appropriate for the dynamic analysis, and the dynamic calculation method will provide some meaningful suggestions for the dynamic analysis of similar projects.

(3) For dynamic strength reduction analysis, there is a lack of rational criteria forselecting the dynamic parameters and determining their ranges. Thus, we considered the static shear strength parameters as the dynamic shear strength parameters in the present study. The experimental study (Li et al. 2008) indicates that the mechanical characteristic and failure criterion of the rock mass under the static and dynamic loads are different. Thus, further study is needed for rational selection of dynamic parameters.

Cheng, Y. M., Lansivaara, T., and Wei, W. B. 2007. Two-dimensional slope stability analysis by limit equilibrium and strength reduction methods.Computers and Geotechnics, 34(3), 137-150. [doi:10.1016/ j.compgeo.2006.10.011]

Clough, R. W., Chang, K. T., Chen, H. Q., and Ghanaat, Y. 1985.Dynamic Interaction Effects in Arch Dams,Report No. UCB/EERC-85/11. Berkeley: Earthquake Engineering Research Center, University of California.

Dai, M. L., and Li, T. C. 2007. Analysis of dynamic stability safety evaluation for complex rock slopes by strength reduction numerical method.Chinese Journal of Rock Mechanics and Engineering, 26(s1), 2749-2754. (in Chinese)

Du, X. L., Zhao, M., and Wang, J. T. 2006. A stress artificial boundary in FEA for near-field wave problem.Chinese Journal of Theoretical and Applied Mechanics, 38(1), 49-56. (in Chinese)

Ge, X. R., Ren, J. X., Li, C. G., and Zheng, H. 2003. 3D-FEM analysis of deep sliding stability of 3# dam foundation of left power house of the Three Gorges Project.Chinese Journal of Geotechnical Engineering, 25(4), 389-394. (in Chinese)

Higdon, R. L. 1986. Absorbing boundary conditions for difference approximations to multi-dimensional wave equation.Mathematics and Computation, 47, 437-459.

Li, L. C., Tang, C. A., Zhu, W. C., and Liang, Z. Z. 2009. Numerical analysis of slope stability based on the gravity increase method.Computers and Geotechnics, 36(7), 1246-1258. [doi:10.1016/j.compgeo.2009. 06.004]

Li, X. B., Zhou, Z. L., Lok, T. S., Hong, L., and Yin, T. B. 2008. Innovative testing technique of rock subjected to coupled static and dynamic loads.International Journal of Rock Mechanics and Mining Sciences, 45(5), 739-748. [doi:10.1016/j.ijrmms.2007.08.013]

Liao, Z. P., Wong, H. L., Yang, B. P., and Yuan, Y. F. 1984. A transmitting boundary for transient wave analyses.Scientia Sincia, Ser. A, 27(10), 1063-1076.

Liu, J. B., Du, Y. X., Du, X. L., Wang, Z. Y., and Wu, J. 2006. 3D viscous-spring artificial boundary in time domain.Earthquake Engineering and Engineering Vibration, 5(1), 93-102. [doi:10.1007/s11803-006-0585-2]

Lysmer, J., and Kuhlemeyer, R. L. 1969. Finite dynamic model for infinite media.Journal of Engineering Mechanics Division, 95(4), 859-878.

Ren, Q. W., Du, C. B., Chen, D. H., and Yuan, J. T. 2009.Dynamic Analysis and Seismic Safety Evaluation of RCC Gravity Dam of Guangzhao Hydropower Project. Nanjing: Department of Engineering Mechanics, Hohai University. (in Chinese)

Saleh, S., and Madabhushi, S. P. G. 2010. Response of concrete dams on rigid and soil foundations under earthquake loading.Journal of Earthquake and Tsunami, 4(3), 251-268. [doi:10.1142/S179343111 0000820]

The State Economic and Trade Commission, P. R. China (SETC). 2000.Design Specification for Concrete Gravity Dams (DL5108-1999). Beijing: China Electric Power Press. (in Chinese)

The State Economic and Trade Commission, P. R. China (SETC). 2001.Specifications for Seismic Design of Hydraulic Structures (DL5073-2000). Beijing: China Electric Power Press. (in Chinese)

Wang, D. S., Zhang, L. J., Xu, J. J., He, M. J., and Zhang, W. D. 2009. Seismic stability safety evaluation ofgravity dam with shear strength reduction method.Water Science and Engineering, 2(2), 52-60. [doi:10.3882/j.issn.1674- 2370.2009.02.006]

Westergaard, H. M. 1933. Water pressures on dams during earthquakes.Transactions of ASCE, 98, 418-433.

Ye, H. L., Zheng, Y. R., Huang, R. Q., Li, A. H., and Du, X. L. 2010. Application of dynamic analysis of strength reduction in the slope engineering under earthquake.Engineering Sciences, 12(3), 41-48. (in Chinese)

Yu, T. T., and Ren, Q. W. 2007. Evaluation of global safety degree of Jinping high arch dam.Chinese Journal of Rock Mechanics and Engineering, 26(4), 787-794. (in Chinese)

Zhang, C. H., Pan, J. W., and Wang, J. T. 2009. Influence of seismic input mechanisms and radiation damping on arch dam response.Soil Dynamics and Earthquake Engineering, 29(9), 1282-1293. [doi:10.1016/j.soil dyn.2009.03.003]

Zhang, J. H., He, J. D., and Fan, J. W. 2001. Static and dynamic stability assessment of slopes or dam foundations using a rigid body-spring element method.International Journal of Rock Mechanics and Mining Sciences, 38(8), 1081-1090. [doi:10.1016/S1365-1609(01)00072-7]

Zheng, Y. R., Tang, X. S., Zhao, S. Y., Deng, C. J., and Lei, W. J. 2009a. Strength reduction and step-loading finite element approaches in geotechnical engineering.Journal of Rock Mechanics and Geotechnical Engineering, 1(1), 21-30. [doi:10.3724/SP.J.1235.2009.00021]

Zheng, Y. R., Ye, H. L., and Huang, R. Q. 2009b. Analysis and discussion of failure mechanism and fracture surface of slope under earthquake.Chinese Journal of Rock Mechanics and Engineering, 28(8), 1714-1723. (in Chinese)

Zhou, W., Chang, X. L., Zhou, C. B., and Liu, X. H. 2008. Failure analysis of high-concrete gravity dam based on strength reserve factor method.Computers and Geotechnics, 35(4), 627-636. [doi:10.1016/j. compgeo.2007.10.005]

This work was supported by the National Basic Research Program of China (973 Program, Grant No. 2007CB714104), the National Natural Science Foundation of China (Grant No. 50779011), and the Innovative Project for Graduate Students of Jiangsu Province (Grant No. CX09B_155Z).

*Corresponding author (e-mail:chdenghong@163.com)

Received Oct. 29, 2010; accepted Jan. 5, 2011