Published: 30 June 2013

Dynamic characteristics of joint surface considering friction and vibration factors based on fractal theory

Xiaopeng Li1
Yamin Liang2
Guanghui Zhao3
Xing Ju4
Haotian Yang5
1, 2, 3, 4, 5School of Mechanical Engineering &Automation, Northeastern University, China
Views 68
Reads 26
Downloads 1367

Abstract

From the point of view of the micro-geometry, most joint surfaces are composed of the rough surfaces with the self-affine fractal characteristics. And the fractal characteristics have a large impact on the dynamic behaviors of the composite structures notably. In this paper the dynamic characteristics of the joint surface are discussed with the friction and vibration coupling effect considered based on fractal theory. Friction and vibration are caused by the rough surfaces. Firstly the expressions of the stiffness and damping model are obtained based on the fractal contact theory. Secondly the contact parameters of the joint surfaces are studied and analyzed with the contact surfaces replaced by the equivalent damping layer with the stiffness and damping model. Thirdly, based on the fractal contact theory, the contact pairs are established by ANSYS to simulate the effect of the characteristics of joint surfaces by the finite element method. Finally the natural frequency and model analysis obtained from theory and experiment of the composite beams are comparatively analyzed. The results show that the established stiffness and damping models are well suitable for the reality of the joint surface.

1. Introduction

In a view of microscopic geometry, most joint surfaces are composed of rough surfaces with self-affine fractal characteristics. In order to create a more accurate dynamic model of joint surfaces, the dynamic characteristics of joint surfaces should be well understood. Many new theories and contact models of joint surfaces are obtained, although the research of the dynamic characteristics largely depends on the classic theories which contain Hertz, G-W, and M-B [1-3]. The G-W statistical model with characteristic of elastic and plastic contact was put forward by Nayak [1]. Majumdar [2] presented the improved M-B model which contains the affine fractal features. Professor Tian [4] proposed a defect of the fractal model where the contact stiffness is an unconditional equation, and then he put forward the modified mathematical formula based on the M-B model. The elastic, elastic-plastic and plastic fractal models of the rough contact surface were established by Zhu [5], and the influence of the interface friction and the elastic-plastic deformation process of the materials are all considered. The research on the deformation characteristics of joint surfaces begins with the contact mechanics point [6-8]. Then the fractal models with the friction which are seldom studied by the dynamics theory [9-10] by now are established gradually [11-12].

Considering the friction of the joint surfaces, the fractal models of the contact stiffness and damping of joint surface are deduced based on the M-B contact model in this paper. Then the joint surface is replaced by the equivalent damping layer based on the obtained models in the process of simulation. The contents of the dynamic characteristics and the finite element analysis are expanded with the modal strain method, which can improve simulation credibility. Finally the reliability of the equivalent method is verified by comparing the results of the modal experiment with the equivalent damping method.

2. Fractal prediction models of the normal contact parameters of joint surface

The friction between the contact surfaces of the mechanical joint surface makes the contact stress different from the classic Hertz solutions for different value and distribution. The influence of friction is ignored in the classic M-B fractal model. So it’s essential to modify the M-B model with the friction and establish the more accurate fractal model. The most basic theory of the contact surface is the Hertz theory and some fractal models are establied based on it. The section of the static contact deformation of two elastomers in Hertz theory is shown in Figure 1.

According to the Hertz theory, the real contact radius is:

1
r=3QR4E1/3,

where E=1-v12E1+1-v22E2-1 is the equivalent elastic modulus. Q is the contact force. E1, E2 are the elastic modulus of two elastomers. v1 and v2 are the Poisson's ratios of the elastomers. And R1, R2 are the curvature radiuses of the two elastomers, separately. The equivalent curvature radius of single peak is R=1R1+1R2-1.

Based on the M-B fractal model, the mathematical model of the surface outline of the undeformed micro-bulge is expressed as Equation (2) [5]:

2
zx=GD-1a1-0.5Dcosπxa0.5, -0.5a0.5<x<0.5a0.5,

where a is the actual contact area, D is the fractal dimension, G is the fractal scale coefficient.

Fig. 1The section of the static contact deformation of two elastomers

The section of the static contact deformation of two elastomers

When x= 0 in the Equation (2), the contact deformation of micro-bulge can be expressed as Equation (3):

3
δ=GD-1a1-0.5D.

Through Equation (2), curvature radius of elastomer can be obtained as shown in Equation (4):

4
R=a0.5DG1-Dπ.

And also the relationship between a, δ and ac, δc is shown in Equation (5):

5
δcδ=aca1-D.

From the Reference [5] and other references, the relation of the biggest contact area (amax) and the distribution of the contact points (n(a)) is:

6
na=D2ψ1-0.5Damax0.5Da-1-0.5D, 0<a<amax,

where ψ is fractal area expansion coefficient.

When the relative sliding is produced in the joint surfaces, the influence of surface friction on the critical deformation between the elastic and plastic state of micro-bulge must be considered. And it can be expressed as follows:

7
δμc=3πpμ4E2R=3.3πkμσy4E2R,

where σyis the yield strength, kμ is the friction coefficient correction factor, kμ=1-0.228μ, 0μ0.3, kμ=0.932e-1.58μ-0.3, 0.3μ0.9. So the critical area between the elastic and plastic states of micro-bulge can be obtained through Equation (3), Equation (4) and Equation (5) expressed as Equation (8):

8
aμc=3.3π0.5kμσy4E21-DG2.

When the micro-bulge is in the elastic deformation state, according to the Hertz theory the normal contact stiffness of a single micro-bulge is:

9
kn=2rE.

The elastic deformation will happen when the contact area is larger than the critical area, while the influnce of the friction is considered. Considering the distribution of the micro-bulge area between joint surfaces and the elastic deformation, the fractal model of the normal contact stiffness of the whole joint surfaces can be deviated based on the Equation (1), Equation (6) and Equation (9):

10
Kn=2D1-Dπ-0.5ψ1-0.5DE(amax0.5-amax0.5Daμc-0.5D+0.5).

Then the dimensionless expression is:

11
Kn*=2D1-Dπ-0.5ψ1-0.5Da*max0.5-a*max0.5Da*μc-0.5D+0.5.

From Equation (11) it shows that the stiffness of the joint surfaces is combined with the critical area of elastic-plastic deformation in friction state. When the fractal characteristic parameters of joint surfaces, such as D, G and material elastic modulus (E), are obtained, the contact stiffness of joint surfaces can be estimated. The accuracy of estimating the actual contact rigidity will be improved by this method.

Then the influence of the friction coefficient on the normal damping of the joint surfaces can be obtained with the simulation software and the results are shown in Figure 2.

Form the Figure 2, the normal contact stiffness of joint surfaces continuously decreased with the increase of the friction coefficient. When the friction coefficient is less than 0.3, the dimensionless normal contact stiffness appears to have linear attenuation with the increase of the friction. And when the friction coefficient is larger than 0.3, the dimensionless normal contact stiffness appears to have the exponential attenuation with the increase of the friction coefficient, and the attenuation speed is reduced quickly with the increase of the actual combination area.

Based on the fractal theory we can study the expression of the contact stiffness and the influence of the friction coefficient on the stiffness from microcosmic angle. And the method to get the contact stiffness with the fractal theory would be improved.

Fig. 2The relationship between the normal contact stiffness and friction coefficient of joint surfaces

The relationship between the normal contact stiffness and friction coefficient of joint surfaces

Now we will deduce the normal contact damping of the joint surfaces considering the effect of friction. When the elastic strain of single micro-bulge happened, the stored elastic energy can be obtained through integrating the contact load and deformation of single micro-bulge as Equation (12):

12
Ue=0δPe(δ)dδ=0δ43ER0.5δ1.5dδ=815ER0.5δ2.5, 0<δ<δμc.

And when the plastic strain of single micro-bulge happened, the irrecoverable deformation energy can be obtained through integrating the contact load and deformation of single micro-bulge:

13
Up=δμcδPp(a)dδ=δμcδKσyπRδdδ=π2KRσyδ2-δμc2, δ>δμc.

According to Equation (2), and computing the integral of Equation (12) and Equation (13) for area respectively, we can obtain the stored elastic energy and the consumption of the plastic deformation of the whole joint surfaces in a contact process.

Then Ue and Up are changed as Equation (14) and Equation (15):

14
Ue=aμcamaxUenada=4D152.5-1.5DEπ-0.5G2D-1ψ1-0.5Damax0.5Damax2.5-1.5D-aμc2.5-1.5D, amax>aμc,
15
Up=0amaxUpnada=D42-DKσyGD-1ψ1-0.5Damax2-0.5D+14KσyGD-1ψ1-0.5Damax-0.5Daμc2, (amaxaμc).

Damping factor (η) is a kind of parameter that describes the damping characteristics of the structure or material and it is equal to the ratio of the dissipative energy and the storage of strain in a movement process:

16
η=Upamax=aμcUeamax>aμc=15×2.5-1.5DK16DEπ-0.5GD-1×σyaμc2-0.5DD2-D+1amax0.5Damax2.5-1.5D-aμc2.5-1.5D.

When the structure is in resonance state, the relationship of the damping factor (η), the normal contact damping (Cn) and damping ratio (ξ) is:

17
η=2ξ=2CnCc,

where Cc is the critical contact damping of joint surfaces. The Cn can be obtained by Equation (17). And then with the simulation method the influence of the friction coefficient on the normal damping of joint surfaces can be studied, the results are shown in Figure 3.

Fig. 3The relationship between the damping factor and friction coefficient of joint surfaces

The relationship between the damping factor and friction coefficient of joint surfaces

Form the Figure 3 the normal contact damping of the joint surfaces continuously increased with the increase of the friction coefficient. While the real contact area of the joint surfaces is becoming bigger, the friction coefficient is becoming bigger too. While the real contact area of the joint surfaces is smaller, the better damping performance will be obtained with changing the friction coefficient.

3. The fractal prediction model of the tangential contact parameters of joint surface

The former researches show that when the effect of dry friction is taken into consideration, the normal stress on the boundary of the contact circle is low, while the tangential stress is close to the infinite. So the boundary of the contact area will slip even under any size of the tangential load. The incomplete sliding will occur on joint surfaces when tangential load is smaller than the maximum static friction stress. In this condition the single point contact area can be divided into no slip region (adhesion circle with the radius c) and partly slip annular region of boundary. They are shown in Figure 4.

Fig. 4Tangential contact of micro-bulge with the effect of friction considered

Tangential contact of micro-bulge with the effect of friction considered

The distribution of shear stress of adhesion region is assumed as follows:

18
τstickr'=-μp0cr1-r'c20.5, 0<r'<c.

The distribution of the shear stress of partly slip region is:

19
τslipr'=up01-r'r20.5, c<r'<r,

where p0 is the maximum contact stress of the Hertz contact center, p0=1.5pm, pm is the average contact pressure.

The displacement of the adhesion region δstick is:

20
δstick=14Gsrπμp01-cr2,

where Gs is the equivalent shear modulus, μ is the friction coefficient.

The total tangential stress of joint surfaces is:

21
Qt'=23πτ2r2-τ1c2=2πr23μp01-cr3,
22
cr=1-Qt'2μπr2p0313=1-Qt'μQ13.

The shear stiffness of the single contact point can be obtained with differential of δstick and taking Equation (21) into Equation (22). It is expressed as Equation (23):

23
kt'=8Gsr1-cr331-cr2=431-cr31-cr2kt.

Only the plastic deformation of the contact point happened under the normal loading, the plastic flowing will happen along the tangential direction under the influence of the tangential forces. In this condition the ability to resist the shear deformation of the contact point is nearly lost. Hence, while obtaining the shear stiffness of joint surfaces, only the shear stiffness of the contact point where elastic deformation happened under the effect of the normal loading is taken into consideration. According to the distribution of the contact area of the fractal joint surfaces, the shear stiffness under the tangential force can be obtained as follows:

24
Kt*'=34π134D1-DQt*'1-1-Qt*'23Qn*1/3ψ1-0.5D×G*1-D6a*max0.5Da*max1-D3-a*μc1-D3,

where Kt'*=Kt'/GsAa0.5 is the dimensionless shear stiffness, Qn*=Qn/EAr is the dimensionless normal loading, Qt*'=Qt'/μQ is the dimensionless tangential loading.

And then with the simulation method the influence of the friction coefficient on the tangential stiffness of joint surfaces can be studied, the results are shown in Figure 5.

The relationships between the tangential stiffness of joint surfaces and the friction factor with the fractal dimension D= 1.25 and D= 1.62 are shown in Figure 5. From the Figure 5 the tendence can be found, the tangential stiffness decreased with the friction factor increasing. But when D= 1.25 the influence of friction factor on the tangential stiffness of joint surfaces will be reduced. At the same time there is no use in improving the tangential stiffness of the joint surfaces by only changing the friction factor.

Besides the friction generating from the slip of the joint surfaces will consume the energy of the system, which shows the dry friction damping characteristics. On the basis of this characteristic, the characteristics of tangential stiffness of the whole joint surfaces under the tangential cyclic load can be obtained.

Assume that the normal load on joint surfaces is Qn, the amplitude is Qt and the period of motion is ϖs. Any two asperities are not only affected by the normal load Q=aQn/Ar, but also by the tangential load Qt'=aQt/Ar. When Qt'<μQ the two asperities are in the state of partial slip, the relationship between the tangential load and relative displacement is shown as ellipse.

Fig. 5The relationships of tangential stiffness and friction factor with different actual contact areas

The relationships of tangential stiffness and friction factor with different actual contact areas

a) Fractal dimension D= 1.25

The relationships of tangential stiffness and friction factor with different actual contact areas

b) Fractal dimension D= 1.62

The energy consumed in one period of the whole joint surfaces is:

25
Ed=0a maxedn(a)da=11.5-0.5D9(μQn)2π0.5D5GsAr2Λψ1-0.5Damax1.5.

Under the normal and tangential loading, the micro sliding area consuming the energy will be formed nearby to the contact point before the macro slipping of joint surfaces occurres. So the dry friction damping under the periodic tangential load can be replaced by the equivalent viscosity damping coefficient.

The equivalent viscosity damping coefficient of the dry friction damping is [6]:

26
Cte=Edπϖsδstick2,

where Ed is the work that the friction does in a whole period.

Substituting Equation (25) and Equation (20) into Equation (26), the expression of Cte is:

27
Cte=64Gsπ0.5-73Dϖs7.5-2.5D34EQnAr2/3×Λψ1-0.5Damax1.5a1+0.5DG1-D2/3.

If the contact radius is changed with the contact area changed, the dimensionless equivalent tangential damping factor of joint surfaces is expressed as follows:

28
Cte*=64π0.5-73D342/37.5-2.5DG*1-D3×Λψ1-0.5DQn*2/3amax*2D-56,

where G*=G2/Aa is the nondimensional parameter of the fractal dimension, Qn*=3Qn/4EAr is the nondimensional parameter of the normal load, amax*=amax/Aa is the biggest single contact area of the nondimensional parameters, Λ is the function about the friction coefficient of joint surfaces.

The factors in Equation (28) have a great influence on the contact characteristics of joint surfaces. And then with the simulation method the influence of the friction coefficient on the tangential damping ratio of joint surfaces can be studied, the results are shown in Figure 6.

Fig. 6The relationship between the tangential damping factor and the friction factor with different contact area

The relationship between the tangential damping factor  and the friction factor with different contact area

From the Figure 6 we can find that the tangential damping factor of joint surfaces becomes more stable with the friction coefficient increasing. The location of the break of Cte*-μ curve will be affected by the size of the actual contact area. The critical friction factor μc will move backwards as the actual contact area increases. The application of the theory of tangential damping factor of joint surfaces obtained in this paper is only applied to the state when micro partial slip occurred between joint surfaces.

4. Application and validation of the fractal contact models

When the machine is working, macro or micro displacement of joint surface will convert part of the mechanical energy into heat and noise consumption through friction, it essentially shows the damping characteristics. Based on the previous researches the joint surface can be equivalent to damping layer. Material constants of the damping layer can be determined based on the concept that the fractal contact stiffness is equal to strain energy. That’s to say, tangential parameters and normal paramenters obtained from the fractal model can be converted into the material properties of the damping layer. The damping layer is shown in Figure 7. The dynamic characteristics of joint surface are simulated with the finite element software by introducing the influence of the friction coefficient for the material properties of the layer.

By material mechanics, the power (W) generated by the normal loading (Q) in the elastic range is shown as follows:

29
W=Q2h2EAa=Q22Kn.

In the linear elastic range the strain energy of the equivalent damping layer of joint surfaces is:

30
Vε=Q2h2EiAa.

Due to the fact that the power (W) generated by the external loading is equal to strain energy stored by the equivalent damping layer, equation is shown as follows:

31
Ei=KnhAa.

In the same way the shear modulus (Gi) of the material of the equivalent damping layer is shown as follows:

32
Gi=KthAa.

The relationship between the shear modulus (Gi), elastic modulus (Ei) and Poisson's ratio (vi) of the material of the equivalent damping layer is shown as follows:

33
νi=Ei2Gi-1.

Because the normal stiffness (Kn) and the tangential stiffness (Kt) of joint surfaces are the functions of the friction coefficient, Equations (30-32) show the relationship between the material coefficient of the layer and the friction coefficient. They also consider the influence of the friction coefficient, so they better match the engineering practice.

By using the stylus profile meter to measure the contact surface of the composite beams in Figure 8, the fractal dimension (D) of its profilogram is 1.4781 and the fractal scale parameters (G) are 0.00016.

Fig. 7The equivalent damping layer of joint surfaces

The equivalent damping layer of joint surfaces

Fig. 8The composite beams for experiment

The composite beams for experiment

When the tightening torque of bolt is 4 Nm, by substituting the fractal parameters of steel plate (D), (G), elastic modulus (E), Poisson's ratio (v), density (ρ) and friction coefficient (μ) into Equation (10) and Equation (23), the normal contact stiffness (Kn) and the tangential contact stiffness (Kt) of joint surfaces are obtained. The parameters are summarized in Table 1.

Table 1The material constants and fractal parameters

μ
E (Pa)
ν
ρ (kg/m3)
D
G
0.12
2.1e11
0.3
7850
1.4781
0.00016

By putting these parameters into Equations (30), (31) and (32), the material constant of the equivalent damping layer is obtained. Considering the surface morphology of joint surfaces and the oxide layer and the plastic deformation layer by loading, the thickness of the equivalent damping layer is taken 8 mm. The material constants are summarized in Table 2.

Table 2The material constants of the equivalent damping layer

Ei (Pa)
νi
ρ (kg/m3)
h (m)
5.5e8
0.35
7850
0.0008

According to the structure of the composite beams in order to accurately mesh the finite element model, which contains the top and the bottom steel plates and the equivalent damping layer, by ignoring the bolts the model is built. In the model the thickness of each steel plate is 6 mm, the equivalent damping layer is 0.8 mm, and each layer material is divided as 40×10×2 by the Solid-45 element with the simulation software, the total number of mesh points is 1200, it is shown in Figure 9. And each model shape and natural frequency are obtained through the method of Block Lanczos [13, 14] of modal analysis. The natural frequencies and loss factors of the composite beams processed by the equivalent damping layer are as summarized in Table 3.

Fig. 9The finite element model of the equivalent damping layer of the composite beams

The finite element model of the equivalent damping layer of the composite beams

Table 3The natural frequency and loss factor of composite beams processed by equivalent damping layer

Modal order
1
2
3
4
5
Natural frequency (Hz)
374.5
912.2
1590
2384
3309
Loss factor
0.086
0.211
0.271
0.283
0.267

The relationship between the modal loss factor and the damping ratio of the system is linear under each resonance frequency. The damping characteristics of the composite beams can be demonstrated by the modal loss factor. The purpose of this experiment is to obtain the natural frequency and the damping ratio of the composite beams. The magnitude-frequency response curve of the composite beams is shown as in Figure 10. And the comparisons of the results of the equivalent damping layers and experiment are listed in Table 4.

Fig. 10The magnitude-frequency curve of the composite beams

The magnitude-frequency curve of the composite beams

From the Table 4 modal parameters of each order of the composite beams with modal test are in agreement with the results, which are obtained by the finite element calculation method of equivalent damping layers with the influence of joint surfaces contact characteristics considered. The error of two kinds of results is about 5 %, and few may be close to 10 %. This method can provide the simulation value of modal loss factor related to structural damping. Above all the finite element calculation method of equivalent damping layers is used for the dynamics analysis of mechanical structure which considers the influence of joint surfaces.

Table 4Comparisons of the results of the equivalent damping layers and experiment

Modal
order
Experimental
natural frequency
Equivalent
damping layers
Error
(%)
Modal experiment
damping ratio (%)
Equivalent
damping layers
1
365
375
2.6
0.4
0.09
2
897
912
1.7
1.2
0.21
3
1691
1590
6.0
1.1
0.27
4
2430
2384
1.9
1.4
0.28
5
3030
3309
9.2
1.5
0.27

5. Conclusions

Based on the M-B model, the fractal model of the normal contact stiffness of joint surfaces is studied considering the friction coefficient between joint surfaces. The fractal model of the normal contact damping and the tangential contact stiffness is established based on the concept of modal loss factor. The fractal model of the tangential contact damping based on the equivalent treatment method of the friction damping is obtained when the surface is in the sliding state. After obtaining a fractal prediction model of the contact parameters of joint surfaces, this paper studies the influence of the friction coefficient on the contact stiffness and damping with simulation method. The simulation results show that the fractal prediction model considering the friction coefficient is certainly right.

The joint surface of the composite beams is equivalent to a damping layer and through this way the material constant of the equivalent damping layer is gained. The modal analyses of the composite beams are made with the modal strain energy method and through this analysis each of natural frequencies and modal loss factors are also gained.

Finally the modal experiment of the composite beams is done, and the modal parameters of joint surfaces and the composite beam are gained. The error of the results between experiment and equivalent damping is very small, so it is feasible to simulate the joint surface with the equivalent method.

References

  • Ni Jun. A perspective review of CNC machine accuracy enhancement through real-time error compensation. China Mechanical Engineering, Vol. 8, Issue 1, 1997, p. 29-34.
  • Jianguo Yang, Yongqiang Ren, Weibin Zhu, et al. Research on on-line modeling method of thermal error compensation model for CNC machines. Chinese Journal of Mechanical Engineering, Vol. 39, Issue 3, 2003, p. 81-84.
  • Shiping Su. Study on the Methods of Precision Modeling and Error Compensation for Multi-Axis CNC Machine Tools. National University of Defense Technology, Beijing, 2002.
  • Xiaoli Li. Modeling and Compensation for the Comprehensive Geometric Errors of NC Machine Tools. Huazhong University of Science and Technology, Wuhan, 2006.
  • Huanlao Liu. Research on the Geometric Error Measurement and Error Compensation of the Numerical Control Machine Tools. Huazhong University of Science and Technology, Wuhan, 2005.
  • R. Hocken, A. J. Simpon, et al. Three dimensional metrology. Annals of the CIRP, Vol. 26, Issue 2, 1997, p. 403-408.
  • P. M. Ferreira, C. R. Liu. An analytical quadratic model for the geometric error of a machine tool. Journal of Manufacturing Systems, Vol. 5, Issue 1, 1986, p. 51-63.
  • S. Yang, J. Yuan, J. Ni. The improvement of thermal error modeling and compensation on machine tools by CMAC neural network. International Journal of Machine Tools and Manufacture, Vol. 36, Issue 4, 1996, p. 527-537.
  • G. Q. Chen, J. X. Yuan, J. Ni. A displacement measurement approach for machine geometric error assessment. International Journal of Machine Tools & Manufacture, Vol. 41, Issue 1, 2001, p. 149-161.
  • Zhiyong Qu, Weishan Chen, Yao Yu. Research on volumetric error identification of guide. Chinese Journal of Mechanical Engineering, Vol. 42, Issue 4, 2006, p. 201-205.
  • Youwu Liu, Libing Liu, Xiaosong Zhao, et al. Investigation of error compensation technology for NC machine tool. China Mechanical Engineering, Vol. 9, Issue 12, 1998, p. 48-52.
  • G. Zhang, R. Yang Ou, B. Lu. A displacement method for machine geometry calibration. Annals of the CIRP, Vol. 37, Issue 1, 1988, p. 515-518.
  • Grimes R. G., Lewis J. G., Simon H. D. A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems. SIAM Journal Matrix Analysis Applications, Vol. 15, Issue 1, 1994, p. 228-272.
  • Hongsong Zhang, Renxi Hu, Shiting Kang. ANSYS 13.0 Finite Element Analysis from Entry to the Master. China Machine Press, 2011, p. 296-318.

About this article

Received
05 March 2013
Accepted
03 June 2013
Published
30 June 2013
Keywords
joint surfaces
contact theory
fractal theory
mechanical dynamics
Acknowledgements

This paper is supported by the National Natural Science Foundation (51275079), Program for New Century Excellent Talents in University (NCET-10-0301) and Fundamental Research Funds for the Central Universities (N110403009).