Response of a Shape Memory Alloy Beam Model under Narrow Band Noise Excitation

Abstract

To describe the hysteretic nonlinear characteristic of the strain-stress relation of shape memory alloy (SMA), a Van-der-Pol hysteretic cycle is applied to simulate the hysteretic loops. Then, the model of a simply supported SMA beam subject to transverse narrow band noise excitation with nonlinear damping was proposed. The deterministic and the stochastic responses are studied, respectively, applying the multiple scale method. The stability of the steady state responses is analyzed by Floquet theory and the moment method. The numerical simulation results quite agree with the theoretical analysis.

1. Introduction

Shape memory alloy (SMA) is a kind of smart materials and is applied in the engineering field widely. There are various mechanical models of SMA proposed in recent decades [1–3]. Most of them are based on thermodynamics theory and micromechanics theory. However, those mechanical models of SMA are mostly shown as equations with subsection function which are hard to be analyzed. Furthermore, SMA beams are widely used in mechanical engineering as sensors and actors. It is of great importance to study the vibrating character of the SMA beam. In the recent decades, some impressive works have been done on the nonlinear vibrating of the SMA beam [4–9]. Many interesting phenomena such as bifurcation, period doubling, jumping, and chaos can be observed. Anyway, the studies mentioned above are all about deterministic excitations. Stochastic excitations widely exist in engineering applications, such as wind loads and seismic waves. By now, few has been done to explore the dynamical responses induced by stochastic excitation. For instance, Yan and Nie [10] studied a SMA superelastic subjected to stationary random excitations.

In this paper, we focus on a SMA beam under narrow band noise excitation. Narrow band noise is a harmonic function with constant amplitude and random frequency, which is a reasonable model for the random excitation or response to engineering systems compared with the harmonic excitations. In Section 2, a Van-der-Pol hysteretic cycle is applied to simulate the hysteretic character of SMA based on the experimental data, and the Galerkin method is used to obtain the one-degree-freedom ordinary differential equation governing the motion of a simply supported SMA beam. Then the multiple scale method is applied to simplify the original function and get the stationary state response in Section 3. Furthermore, in Section 4, the Ito-type stochastic differential equation is obtained by linearization, and the responses and stability of both the noise-free system and the noise included system are studied using Floquet theory and the moment method. Finally, the numerical simulation is carried out to prove the theoretical analysis.

2. Modeling of the SMA Beam

The experimental hysteretic nonlinearity strain-stress curve of SMA [11] is shown in Figure 1. The length of the Ti-Ni SMA specimens is 14 cm, the diameter of the specimens is 1.94 mm, and transformation temperature is −10°C. The experimental equipment is Instron 8501 plus type fatigue testing machine. The testing temperature is set at 25°C, the exciting frequency is set as 1 Hz, and the cycle number is set as 40. According to [11], it is clear that when the cycle number is limited, the shape of the hysteretic loop does not change much, even the exciting frequency is different. As we know, the SMA is usually used as damper in some structure, for example, aseismic structure. In fact, the effective imposing time of aseismic structure is less than 20 seconds, so the number of the hysteretic loops is limited. Then we can assume that the shape of the hysteretic loop remains unchanged.

Figure 1 

The experimental stain-stress curve.

Now, Van-der-Pol hysteretic cycle model was introduced to describe the hysteretic characteristic of SMA. Before the numerical simulation, some assumptions must be given. Firstly, because the difference of the curves is small when the cycle number is limited, we take the average stress value as the final experimental data (see in Figure 2). Secondly, the residual strain is obviously small, so it could be ignored for convenience. Thirdly, considering that the bending deformation is small in most working conditions, the tensile-compression asymmetry of SMA behavior is also ignored. In short, in this paper we focus on the method of simulating the hysteretic loop. The initial Van-der-Pol hysteretic cycle model describes hysteretic cycle which is symmetrical about the initial point (0, 0). Since the curve is obtained by isostrain rate tensile test, we define a parameter 𝑥̇ as loading rate, which is chosen as 1 while loading (−1 as unloading). This parameter just has mathematical meaning. Then the model can be shown as follows: 𝑦=𝑓(𝑥)=𝑓₀(𝑥)𝑥+𝑎1−𝑏²𝑥̇=𝑏₁𝑥+𝑏₂𝑥³𝑥+𝑎1−𝑏²𝑥̇, (1) where 𝑓₀(𝑥) is skeleton curve of hysteretic cycle and 𝑎 and 𝑏 are coefficients. The parameters 𝑎 and 𝑏 are used to change the shape of the hysteretic loop. Supposing the strain-stress curve of SMA is symmetrical about the point 𝐺 (𝜀₀,𝜎₀), one may substitute 𝑦=𝜎−𝜎₀ and 𝑥=𝜀−𝜀₀ into (1). Considering the residual strain of SMA should be avoided in simulation,that is, when 𝜀=0, there is 𝜎=0. Then the final expression of the stress-stain relation of SMA is 𝜎=𝑎₁𝜀+𝑎₂𝜀²+𝑎₃𝜀³+𝑎₄𝜀−𝑎₅𝜀²𝜀,̇ (2) where 𝑎₁=𝑏₁+3𝑏₂𝜀²₀, 𝑎₂=−3𝑏₂𝜀₀, 𝑎₃=𝑏₂, 𝑎₄=2𝑎/𝜀₀, and 𝑎₅=𝑎/𝜀²₀.

Figure 2 

The fitting diagram of stain-stress curve.

Then the strain-stress curve of SMA obtained by numerical fitting can be shown in Figure 2.

One hinged-hinged SMA beam is shown in Figure 3 with thickness ℎ, length 𝐿, and width 𝑏. The axial force may be expressed in the form 𝑁=𝑝₀. The transverse narrow band noise excitation is in the form 𝐹𝜂(𝑡)sin((𝜋/𝐿)𝑥). Here, only the first vibration mode is considered, as we think it will have enough accuracy for most engineering applications. 𝜂(𝑡) is in the form of cos(Ω₀𝑡+𝛿^∗𝑊(𝑡)), where Ω₀ denotes the center frequency, 𝑊(𝑡) denotes the standard Wiener process, and 𝛿^∗ denotes the small strength of the 𝑊(𝑡).

Figure 3 

The hinged-hinged SMA beam.

Considering the boundary conditions at 𝑥=0, 𝑤=0, 𝑤_𝑥𝑥=0 and at 𝑥=𝐿, 𝑤=0, 𝑤_𝑥𝑥=0, the transverse deflection 𝑤 is expressed as 𝑤(𝑡,𝑥)=𝑢(𝑡)sin(𝜋𝑥/𝐿), where 𝑢(𝑡) is the amplitude. Consider the small geometrical deformation condition 𝜀=−𝑦𝜕²𝑤/𝜕𝑥².

The bending moment 𝑀 is presented as 𝑀=𝑏^ℎ/2_−ℎ/2−𝑦𝑎₁𝜀+𝑎₂𝜀²+𝑎₃𝜀³+𝑎₄𝜀−𝑎₅𝜀²𝜀̇𝑑𝑦=−𝐼₁𝑎₁𝜋²𝜋𝑢sin𝐿𝑥−𝐼₃𝑎₃𝜋⁶𝑢³sin³𝜋𝐿𝑥+𝐼₃𝑎₅𝜋⁶𝑢²sin³𝜋𝐿𝑥𝑢̇, (3) where 𝐼₁=𝑏ℎ³/(12𝐿²) and 𝐼₃=𝑏ℎ⁵/80𝐿⁶.

The dynamical motion equation [8] is expressed as 𝜕²𝑀𝜕𝑥²+𝑁𝜕²𝑤𝜕𝑥²+𝑐𝜕𝑤+𝜕𝑡𝜌𝐴𝜕²𝑤𝜕𝑡²=𝐹𝜂(𝑡), (4) where 𝑐 is the linear damping coefficient, 𝜌 is the density of the SMA, and 𝐴 is the area of the cross-section of the beam. Substituting (2) and (3) into (4) and applying the Galerkin approach, one obtains a differential equation of motion governing the deflection of a beam as follows: ̈𝑎𝑢+¹𝐼₁𝜋⁴−𝑝₀𝜋²𝜌𝐴𝐿²𝑐𝑢+−𝜌𝐴3𝑎₅𝐼₃𝜋⁸4𝜌𝐴𝐿²𝑢²𝑢̇+3𝑎₃𝐼₃𝜋⁸4𝜌𝐴𝐿²𝑢³=4Ω𝜋𝜌𝐴𝐹cos₀𝑡+𝛿^∗𝑊(𝑡). (5) Considering 𝜔₀=((𝑎₁𝐼₁𝜋⁴−𝑝₀𝜋²)/𝜌𝐴𝑙²), 𝑡^∗=𝜔₀𝑡, 𝑥=𝑢/𝐿, Ω=Ω₀/𝜔₀, and 𝛿=√𝜔₀𝛿^∗, it is easy to know that 𝑥̇𝑡^∗=𝑑𝑥𝑑𝑡^∗=1𝜔₀𝐿𝑑𝑢,𝑥̈𝑡𝑑𝑡^∗=1𝜔₀²𝐿𝑑𝑢.𝑑𝑡 (6)

Substituting (6) and no dimensional variables into (5) and dropping the over star yield the dimensionless function as follows: 𝑥̈+𝑥+𝜀2𝜇−𝛽𝑥²𝑥̇+𝜀𝛼𝑥³=𝜀𝑓cosΩ𝑡+𝛿𝑊(𝑡), (2) where 𝜇=(𝑐/2𝜌𝐴𝜔₀), 𝛽=(3𝑎₅𝐼₃𝜋⁸/4𝜌𝐴𝜔₀), 𝛼=(3𝑎₃𝐼₃𝜋⁸/4𝜌𝐴𝜔²₀), 𝑓=(4/𝜋𝜌𝐴𝜔²₀𝐿)𝐹, and 𝜀 is the mark for small perturbation terms.

3. Multiple Scale Method

A uniformly approximate solution of (2) is sought in the form 𝑥(𝑡,𝜀)=𝑥₀𝑇₀,𝑇₁+𝜀𝑥₁𝑇₀,𝑇₁+⋯, (8) where 𝑇₀=𝑡, 𝑇₁=𝜀𝑡 are fast and slow scale, respectively.

In this paper only the first-order uniform expansion of the solution 𝑥₀(𝑇₀,𝑇₁) is discussed. By denoting 𝐷₀=𝜕/𝜕𝑇₀ and 𝐷₁=𝜕/𝜕𝑇₁, the ordinary-time derivatives can be transformed into partial derivatives as 𝑑𝑑𝑡=𝐷₀+𝜀𝐷₁𝑑+⋯²𝑑𝑡²=𝐷²₀+𝜀2𝐷₀𝐷₁+⋯. (9) Substituting (8) and (9) into (2) and comparing coefficients of 𝜀 with equal powers, one obtains the following equations: 𝐷²₀𝑥₀+𝑥₀𝐷=0,₂₀𝑥₁+𝑥₁=−2𝐷₀𝐷₁𝑥₀−𝛼𝑥³₀−2𝜇−𝛽𝑥²₀𝐷₀𝑥₀+𝑓cosΩ𝑇₀+𝛿𝑊(𝑡). (10) The general solution of (10) can be expressed as 𝑥₀𝑇₀,𝑇₁=𝐴𝑇₁exp𝑖𝑇₀+cc, (11) where the cc is the complex conjugate of its preceding terms and 𝐴(𝑇₁) is the slowly varying amplitude of the response. Substituting (9) and (11) into (10), one obtains 𝐷²₀𝑥₁+𝑥₁=−2𝐷₁𝐴𝑖e^𝑖𝑇0−𝛼𝐴³e^𝑖3𝑇0+3𝐴²𝐴e^𝑖𝑇0−2𝜇𝑖𝐴e^𝑖𝑇0+𝛽𝑖𝐴³e^𝑖3𝑇0+𝑖𝐴²𝐴e^𝑖𝑇0+12𝑓e^{𝑖(𝑇0+𝛿𝑊(𝑇1))}+cc, (12) where the over bar stands for the complex conjugate, 𝛿=𝛿/√𝜀. For the Wiener process 𝑊(𝑡), 𝐸[𝑊(𝑡)]=0, 𝐸[𝑊²(𝑡)]=𝑡, where 𝐸[⋅] denotes the expectation, one has 𝛿𝑊(𝑡)=𝛿√𝜀𝑊(𝜀𝑡)=𝛿𝑊𝑇₁. (13) In this paper we focus on the primary resonances of the system (2). Introducing the detuning parameter 𝜎₀ as Ω=1+𝜀𝜎₀, one has Ω𝑡=𝑇₀+𝜎₀𝑇₁. (14) Using (14) and eliminating the secular terms of (12) yield −2𝐷₁𝐴𝑖−3𝛼𝐴²𝐴−2𝜇𝑖𝐴+𝑖𝛽𝐴²1𝐴+2𝑓e^{𝑖(𝜎0𝑇1+𝛿𝑊(𝑇1))}=0. (15)Expressing 𝐴(𝑇₁) into the polar form𝐴𝑇₁=𝑎𝑇₁exp𝑖𝜑𝑇₁. (16) Substituting (16) into (15) and separating the real and imaginary parts of (15), one obtains 𝑑𝑎𝑑𝑇₁1=−𝜇𝑎+8𝛽𝑎³−12𝑎𝑓sin𝜂,𝑑𝜂𝑑𝑇₁=−𝜎₀3𝑎+8𝛼𝑎³−12𝑓cos𝜂−𝛿𝑎𝑊^′𝑇₁, (17) where 𝜂=𝜑−𝜎₀𝑇₁−𝛿𝑊^′(𝑇₁) and (17) is the first-order equation governing the modulation of the amplitude and phase.

The steady state response (𝑑𝑎/𝑑𝑇₁)=(𝑑𝜂/𝑑𝑇₁)=0 is decided by −𝜇1𝑎+8𝛽𝑎³=12𝑓sin𝜂,−𝜎₀3𝑎+8𝛼𝑎³−𝛿𝑎𝑊^′(𝑇)=12𝑓cos𝜂, (18) where 𝑎 and 𝜂 denote the steady amplitude and phase.

4. Responses and Stability

Then the Floquet theory is applied to investigate the stability of the steady state response by linearization. Assuming 𝑎=𝑎+Δ𝑎𝜂=𝜂+Δ𝜂, where Δ𝑎,Δ𝜂 are the small perturbation terms, substituting the above equations into (17), and ignoring the high order terms, the linearization equation (20) is expressed as 𝑑Δ𝑎𝑑𝑇₁3=−𝜇Δ𝑎+8𝛽𝑎²1Δ𝑎−2𝑓cos𝜂Δ𝜂,𝑑Δ𝜂𝑑𝑇₁=34𝛼1𝑎+𝑎²12𝑓cos𝜂−𝛿𝑊^′(𝑡)+1Δ𝑎𝑎12𝑓sin𝜂Δ𝜂. (19) Substituting (18) into (19) yields the Ito-type stochastic equation: 3𝑑Δ𝑎=−𝜇+8𝛽𝑎²Δ𝑎+𝜎₀3𝑎−8𝛼𝑎³Δ𝜂𝑑𝑇₁,𝑑Δ𝜂=9𝛼8𝜎𝑎−𝑎1Δ𝑎+−𝜇+8𝛽𝑎²Δ𝜂𝑑𝑇₁−𝛿𝑑𝑊(𝑡), (20) where 𝑑𝑊(𝑡) denotes the standard Gauss white noise.

4.1. Noise-Free Condition

Assuming 𝛿 is sufficiently small, we determine the response of system (2) for the noise-free condition (𝛿=0). In this case, (18) can lead to the following frequency response formula: 𝐺=−𝜇1𝑎+8𝛽𝑎³²+−𝜎₀3𝑎+8𝛼𝑎³²−𝑓²4=0. (21)

Obviously, in some proper parametrical conditions, there are multiple solutions coexisting of (21). Some of them are stable and the others are unstable. It is necessary to know how the system's parameters influence the numbers of the solution. The singularity theory [12] is applied to study the condition for multisolutions of (21). Choosing (21) as the bifurcation function and detuning parameter 𝜎₀ amplitude 𝑓 as universal unfolding parameters, the Hysteresis point set is given by 𝐺𝑎,𝜎₀,𝑓=𝐺_𝑎𝑎,𝜎₀,𝑓=𝐺_𝑎𝑎𝑎,𝜎₀,𝑓=0 (22) which yields 3𝑓16²𝜇𝛽+3𝛼𝜎₀−𝜇²+𝜎₀²²=0. (23) When the parameter 𝜎₀ is increased across the line expressed by (23) for a set of given other parameters, the number of the solutions of (21) bifurcates from one to three.

For the triple-solution condition, it is necessary to know the stable conditions for the steady state solutions. The Jacoby matrix is a usually used way to judge the stability, which is expressed as follows: 3det−𝜇+8𝛽𝑎²−𝜆𝜎₀3𝑎−8𝛼𝑎³9𝛼𝑎8−𝜎₀𝑎𝛽8𝑎²−𝜇−𝜆=0. (24)

The characteristic equation can be expressed as 𝜆²−𝑃𝜆+𝑄=0, (25) where 1𝑃=2𝛽𝑎²𝜇−2𝜇𝑄=²+𝜎²₀−𝜇𝛽+3𝛼𝜎₀2𝑎²+27𝛼²+3𝛽²64𝑎⁴. (26) The roots of (25) are 𝜆_1,2=𝑃±𝑃²−4𝑄2. (27) Obviously, the stable condition of the system (19) without noise is 𝑃<0,𝑄>0 (condition 1), which means that the solution is realizable by numerical simulation. And the condition for existing complex roots is 𝑃²−4𝑄<0 (condition 2), which means that there exist periodical vibrations.

In this paper, the parameters are chosen as 𝜇=0.1, 𝛽=0.1, 𝛼=0.5, 𝑓=0.2, and 𝜀=0.1. By function (23), one knows that if the detuning parameter 𝜎₀ is larger than 0.193, the number of steady state solutions bifurcates from one to three. In Figure 4(a), one can observe clearly that when 𝜎₀ is chosen as 0.22, there are three steady state solutions coexisting. The largest one denoted as 𝑎_𝐿 and the smallest one denoted as 𝑎_𝑆 are stable sinks, while the middle one denoted as 𝑎_𝑢 is unstable saddle. In Figure 4(b), when 𝜎₀ is chosen as 0.18, there is only steady state solution which is a stable sink in the phase plane.

(a)

(b)

Figure 4 

The steady state solutions of the noise-free equation (17) (a) 𝜎0=0.22 (b) 𝜎0=0.18.

The numerical simulation is carried out to prove the theoretical analysis. In Figure 5, the black dots are numerical solutions and the black line is the theoretical solution given by (21). The red shadow region is obtained by 𝑄<0 (condition 1), which means that, in this region, the unstable solution cannot be obtained by simulation. While the red line decided by 𝑃²−4𝑄=0 (condition 2) illustrates that there are no periodical solutions between the two red lines, one can find that the theoretical analysis quite agrees with the numerical solutions.

Figure 5 

Frequency response (black dots: numerical simulations, black line: theoretical solution, red shadow: unstable region, and red lines: no periodical solution region).

When there are two stable steady state solutions coexisting, one of them has bigger amplitude than the other. Which kind of amplitude the system will act is decided by the initial conditions. In Figure 6, the 𝑦-axes mean the nondimensional displacement of the beam, the 𝑦̇-axes mean the nondimensional velocity, then the domain of attraction is obtained by simulation, and the black region is the set of initial points leading to large amplitude, while the white region is the domain of attraction for small amplitude.

Figure 6 

Domain of attraction to the large amplitude.

4.2. Noise Included Condition

For the noise included condition, that is, 𝛿≠0, the steady state response given by (18) is needed to study its stochastic response. The moment method [13–15] is an effective way to obtain the steady state moments 𝐸[Δ𝑎],𝐸[Δ𝑎²]. It is obvious that (𝑑𝐸[Δ𝑎]/𝑑𝑇₁)=(𝑑𝐸[Δ𝜂]/𝑑𝑇₁)=0.

Tracking expectation on both sides of (20), one obtains 𝐸[Δ𝑎]=𝐸Δ𝜂=0. (28) Similarly, we have (𝑑𝐸[Δ𝑎²]/𝑑𝑇₁)=(𝑑𝐸[Δ𝜂²]/𝑑𝑇₁)=(𝑑𝐸[Δ𝑎Δ𝜂]/𝑑𝑇₁)=0. Applying Ito's rule of stochastic differentiation, one obtains 34𝛽𝑎²−2𝜇𝐸Δ𝑎²+23𝑎𝜎−4𝛼𝑎³𝐸9Δ𝑎Δ𝜂=0,𝑎𝛼4−2𝜎𝑎𝐸+𝛽Δ𝑎Δ𝜂𝑎²4−2𝜇𝐸Δ𝜂²=−𝛿²,9𝛼8𝜎𝑎−𝑎𝐸Δ𝑎²+𝛽2𝑎²𝐸+−−2𝜇Δ𝑎Δ𝜂3𝛼8𝑎³+𝑎𝜎𝐸Δ𝜂²=0. (29) The solutions of (29) are 𝐸Δ𝑎²=𝛿²3𝑎³𝛼−8𝑎𝜎²×4𝜇−𝑎²𝛽3𝑎⁴9𝛼²+𝛽²−32𝑎²𝜇𝛽𝜇+3𝛼𝜎+64²+𝜎²⁻¹,𝐸Δ𝑎Δ𝜂=𝛿²3𝑎²𝛽−8𝜇3𝑎³𝛼−8×𝑎𝜎4𝜇−𝑎²𝛽3𝑎⁴9𝛼²+𝛽²−32𝑎²𝜇𝛽𝜇+3𝛼𝜎+64²+𝜎²⁻¹,𝐸Δ𝜂²=𝛿²3𝑎⁴9𝛼²+4𝛽²−16𝑎²5𝛽𝜇+6𝛼𝜎+642𝜇²+𝜎²×4𝜇−𝑎²𝛽3𝑎⁴9𝛼²+𝛽²−32𝑎²𝜇𝛽𝜇+3𝛼𝜎+64²+𝜎²⁻¹. (30) Combining (28) and (30), one obtains 𝐸[𝑎]=𝑎,𝐸𝑎²=𝑎²+𝐸Δ𝑎²𝐸𝜂=𝜂,𝐸𝜂²=𝜂²+𝐸Δ𝜂². (31) Form (30) it is easy to see that, for 𝐸[Δ𝑎²]≥0, 𝐸[Δ𝜂²]≥0, there must be 4𝜇−𝑎²𝛽≥0 and 3𝑎²(9𝛼²+𝛽²)−32(𝛽𝜇+3𝛼𝜎₀)+64(𝜇²+𝜎₀²). One finds that the existing condition for second-order moment is the same as the stability condition for the first-order moment. By the method proposed by Zhu [16], the numerical simulation is carried out with small noise density 𝛿=0.05. When 𝜎₀=0.22, the results are shown in Figures 7 and 8. For comparison, the noise-free condition is also given in Figure 7. The initial points chosen as (1,2,.5) can lead to the phase portraits as Figures 7(a) and 7(b). For the noise-free condition, the phase portrait is a limit cycle (see Figure 7(a)), while, for the noise included system, the limit cycle is a diffused one (see Figure 7(b)).

(a)

(b)

Figure 7 

Phase plot for the noise-free and noise included system.

Figure 8 

The frequency response with noise (black dots: numerical simulations and line: theoretical results).

Figure 8 shows the frequency response with noise of this system. Obviously, with small enough 𝛿, the theoretical analysis also approximately agrees with the numerical test. It is found when the exciting frequency is low, the numerical solution agrees with the theoretical solution and when the exciting frequency becomes larger, the numerical solution seems to be not so good compared with the theoretical solution.

5. Summary and Conclusion

In this paper, a Van-der-Pol hysteretic cycle is used to simulate the hysteretic character of SMA, which is suitable and convenient for modeling. Then, a nonlinear vibration model of simply supported SMA beam under narrow band noise excitation is built and approximately solved by the multiple scale method. The steady state responses are obtained for both noise-free condition and noise included condition, respectively. After that, the stability of the steady state is analyzed. Some complex dynamical phenomena can be observed, which is helpful for industrial application.

When the strength of the disturbing noise is small enough, the original system can be simplified to an Ito-type differential equation. The singularity theory is applied to study how the parameters influence the number of stationary solutions. It is found that, with some proper parameters there exist three steady state solutions; two of them are stable and realizable, while the other one is unstable. The initial conditions decide which kind of amplitude can appear. The moment method is used to obtain the first-order and second-order moments of the random steady state responses. Similar results are found compared with the deterministic condition. But the difference is that the limit cycle is diffused by the noise. The results of numerical simulation approximately agree with the analysis.

Conflict of Interests

The author declares that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 11272229 and Grant no. 11302144) and the Science Foundation of Tianjin education committee (Grant no. 20120902).