1 Introduction

Shape Memory Alloy (SMA) is a special material. Since Greninger and Mooradian^[1] observed the shape memory effect on SMA, this material has attracted considerable attention in many scientific fields, such as aerospace^[2], aviation^[3], and medical application^[4]. The characteristic ^[5-6] and constitutive model ^[7-8] of SMA were developed in the early 20^th century. Subsequently, voluminous literature has focused on the analysis of the solution to SMA model in determinative cases. For an SMA actuator, the analysis results have been proposed^[9], which were applied to the field of mechanical design^[10]. In view of the two-point distribution technique, the exact solution to SMA's heat exchanger was received^[11].

In stochastic cases, the research of probability response is an important topic for nonlinear dynamical systems. In recent years, many methods^[12-14] have been developed to solve the problem of probability response, including stochastic averaging method, equivalent nonlinear method^[15], and Monte Carlo method. In Ref. [16], super-elastic system's response under white noise excitation was predicted by the Monte Carlo method. However, the Monte Carlo method involves huge amount of computation and complex algorithms. The stochastic averaging method is a powerful approximate technique for the prediction of probability response. Feng^[17-18] applied stochastic averaging method to explore the response of Duffing oscillator under white noise. Bifurcation phenomenon is one of the focus topics in dynamical system. In the nonlinear vibro-impact system^[19] and SMA system^[20], the stochastic averaging method is applied to analyze bifurcation.

Non-smooth factors are widely existent in engineering fields. In previous literature, the methods of non-smooth transformation and stochastic averaging method were rarely studied for the SMA beam model. Thus, this paper focuses on discussing the probability response and bifurcation for the SMA beam model using stochastic averaging method. Section 2 introduces the SMA beam model with rigid constraint. In view of the non-smooth transformation, the stochastic averaging method is carried out to obtain the average Fokker Planck Kolmogorov (FPK) equation in Section 3. Section 4 shows the numerical results to verify the analytical results. Moreover, the critical values of stochastic bifurcation are investigated in detail.

2 Shape Memory Alloy Beam Model

Consider the model of SMA beam^[21]with constraint in Fig. 1. l, h, b represent the length, thickness, and width of the beam, respectively. The model is excited by

$ F = {s_0} + s\zeta (t) $

(1)

where s₀ represents generalized displacement at the initial moment, ζ(t) is the standard Gaussian white noise, and t is time. ζ(t) satisfies

$ \begin{array}{*{20}{c}} {E[\zeta (t)] = 0}\\ {E[\zeta (t + \tau )\zeta (t)] = \delta (\tau )} \end{array} $

(2)

where τ is the time difference, E is the mean, and δ(τ) represents Dirac delta function.

The motion equation for the beam is as follows:

$ \frac{{{\partial ^2}M}}{{\partial {x^2}}} + F\frac{{{\partial ^2}\omega }}{{\partial {x^2}}} + c\frac{{\partial \omega }}{{\partial t}} + \rho A\frac{{{\partial ^2}\omega }}{{\partial {t^2}}} = 0 $

(3)

where M is bending moment of the beam, c represents linear damping, ω is lateral displacement of the beam, and ρA is the mass of the beam per unit length.

Then, the stochastic nonlinear motion equation of SMA beam can be expressed as

$ \ddot x - kx + \alpha {x^3} + \left( {\tilde \mu - \tilde \gamma {x^2} + \tilde \beta {{\dot x}^2}} \right)\dot x = \tilde f\zeta (t) $

(4)

where

$ k = {{\rm{ \mathsf{ π} }}^4} - \frac{{{s_0}{{\rm{ \mathsf{ π} }}^2}}}{{{a_1}{I_1}}},\alpha = \frac{{3{a_3}{I_3}{{\rm{ \mathsf{ π} }}^8}}}{{4{a_1}{I_1}}},\tilde \mu = \frac{{cl}}{{{{\left[ {\rho A{a_1}{I_1}} \right]}^{1/2}}}} $

$ \tilde \gamma = \frac{{3{a_5}{I_3}{{\rm{ \mathsf{ π} }}^8}}}{{4{{\left[ {\rho A{a_1}{I_1}} \right]}^{1/2}}l}},\tilde \beta = \frac{{{c^2}{l^2}}}{{\rho A{a_1}{I_1}}} $

$ \tilde f = \frac{{4{l^2}}}{{{a_1}{I_1}{\rm{ \mathsf{ π} }}}},{I_1} = \frac{{b{h^3}}}{{12{l^2}}},{I_3} = \frac{{b{h^5}}}{{80{l^6}}} $

$ {a_1} = {b_1} + 3{b_2}\varepsilon _0^2,{a_3} = {b_2},{a_5} = \frac{a}{{\varepsilon _0^2}} $

ε₀ is the symmetric center of the ring, and s represents generalized displacement.

Fig. 2 shows the sectional view of the model. The dimensionless equation of motion is

$ \left\{ {\begin{array}{*{20}{l}} {\ddot x - kx + \alpha {x^3} + \left( {\tilde \mu - \tilde \gamma {x^2} + \tilde \beta {{\dot x}^2}} \right)\dot x = \tilde f\zeta (t),x > 0}\\ {{{\dot x}_ + } = - R{{\dot x}_ - },x = 0} \end{array}} \right. $

(5)

where k, α, $ \tilde{\mu}, \quad \tilde{\gamma}, \quad \tilde{\beta}$ are constants, f is noise amplitude, and R is restitution coefficient, which is usually decided by the material of the collision surface and satisfies 0 < R≤1. Generally, a collision may cause energy loss for the system when 0 < R < 1.

3 Probability Response

Non-smooth transformation is introduced^[22-25]:

$ \begin{array}{*{20}{c}} {x = {x_1} = |y| + \Delta ,\dot x = {x_2} = \dot y\text{sgn}(y)}\\ {\ddot x = \ddot y\text{sgn}(y),\text{sgn}(y) = \left\{ {\begin{array}{*{20}{r}} {1,y > 0}\\ { - 1,y < 0} \end{array}} \right.} \end{array} $

(6)

With Δ=0 and Eq.(6), Eq.(5) can be rewritten as

$ \begin{array}{*{20}{c}} {\ddot y - ky + \alpha {y^3} + \left( {\tilde \mu - \tilde \gamma {y^2} + \tilde \beta {{\dot y}^2}} \right)\dot y + }\\ {(1 - R)\dot y|\dot y|\delta (y) = \tilde f \text{sgn} (y)\zeta (t)} \end{array} $

(7)

where δ(y) is Dirac delta function.

Compare the energy difference between the stochastic excitation and the damp. If system energy is powerful within a period, Eq.(7) can be regarded as a quasi-conservative system. Assume $\tilde{\mu}=\varepsilon \mu, \tilde{\gamma}=\varepsilon \gamma $, $\tilde \beta = \varepsilon \beta , \tilde f = \sqrt \varepsilon f $, R=1－εr. ε is a small parameter scale. Let y₁=y, y₂=$ \dot y$, Eq.(7) can be rewritten as

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot y}_1} = {y_2}}\\ {{{\dot y}_2} = k{y_1} - \alpha y_1^3 - \varepsilon \left( {\mu - \gamma y_1^2 + \beta y_2^2} \right){y_2} - }\\ {\;\;\;\;\;(1 - R){y_2}\left| {{y_2}} \right|\delta \left( {{y_1}} \right) + \sqrt \varepsilon f \text{sgn} \left( {{y_1}} \right)\zeta (t)} \end{array}} \right. $

(8)

Using the Itô formula, Eq.(8) is reduced to the following Itô differential equation:

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{d}}{y_1} = {y_2}{\rm{d}}t}\\ {{\rm{d}}{y_2} = \left[ {k{y_1} - \alpha y_1^3 - \varepsilon \left( {\mu - \gamma y_1^2 + \beta y_2^2} \right){y_2} - } \right.}\\ {\left. {\left( {1 - R} \right){y_2}\left| {{y_2}} \right|\delta \left( {{y_1}} \right) + \frac{1}{2}\varepsilon {f^2}y_1^2} \right]{\rm{d}}t + }\\ {\sqrt \varepsilon f \text{sgn} \left( {{y_1}} \right){\rm{d}}W(t)} \end{array}} \right. $

(9)

in which W(t) is the Wiener process.

When ε=0 in Eq.(8), the unperturbed system is obtained. The corresponding energy function and potential function are

$ H = \frac{{y_2^2}}{2} - \frac{k}{2}y_1^2 + \alpha \frac{{y_1^4}}{4} $

(10)

$ U\left( {{y_1}} \right) = - \frac{k}{2}y_1^2 + \alpha \frac{{y_1^4}}{4} $

(11)

The derivative of the energy Eq.(10) on the time t is

$ \begin{array}{l} \dot H = {y_2}\left[ { - \varepsilon \left( {\mu - \gamma y_1^2 + \beta y_2^2} \right){y_2} - (1 - } \right.\\ \;\;\;\;\;\;\left. {R){y_2}\left| {{y_2}} \right|\delta \left( {{y_1}} \right) + \sqrt \varepsilon f \text{sgn} \left( {{y_1}} \right)\zeta (t)} \right] \end{array} $

(12)

According to the Itô formula, the Itô stochastic differential equation of Eq. (12) can be expressed as

$ \begin{array}{l} {\rm{d}}H = \left[ {{y_2}\left[ { - \varepsilon \left( {\mu - \gamma y_1^2 + \beta y_2^2} \right){y_2} - } \right.} \right.\\ \;\;\;\;\;\;\;\;\left. {\left. {(1 - R){y_2}\left| {{y_2}} \right|\delta \left( {{y_1}} \right)} \right] + \frac{1}{2}\varepsilon {f^2}} \right]{\rm{d}}t + \\ \;\;\;\;\;\;\;\;\sqrt \varepsilon f{y_2} \text{sgn} \left( {{y_1}} \right){\rm{d}}W(t) \end{array} $

(13)

It is known from the Khassminskii theorem^[26], y₁ and y₂ are two fast mutative processes, and H(t) is a slow mutative process, which is approach to a Markov process as ε→0. Therefore, the averaging Itô stochastic diffusion process of the Markov process can be described as

$ {\rm{d}}H = \langle m(H)\rangle {\rm{d}}t + \langle \sigma (H)\rangle {\rm{d}}W(t) $

(14)

where W(t) is the Wiener process, 〈m(H)〉 and 〈σ(H)〉 are the mean drift and mean diffusion coefficients:

$ \begin{array}{l} \langle m(H)\rangle = \frac{1}{{{T_{\frac{1}{4}}}(H)}}[ - \varepsilon \mu B(H) + \varepsilon \gamma S(H) - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\varepsilon \beta Q(H) + (R - 1)H + \frac{1}{2}\varepsilon {f^2}{T_{\frac{1}{4}}}(H)} \right] \end{array} $

(15)

$ \left\langle {{\sigma ^2}(H)} \right\rangle = \frac{{\varepsilon {f^2}}}{{{T_{\frac{1}{4}}}(H)}}B(H) $

(16)

$ \begin{array}{l} {T_{\frac{1}{4}}}(H) = \frac{{T(H)}}{4} = \int_0^A {\frac{1}{{\sqrt {2H - 2U\left( {{y_1}} \right)} }}{\rm{d}}{y_1}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{{\sqrt {\alpha {A^2} - k} }}{E_1}(K) = \frac{1}{{\sqrt M }}{E_1}(K) \end{array} $

$ B(H) = \frac{{2\sqrt M }}{{3\alpha }}\left[ {(M - N){E_1}(K) + (2N - M){E_2}(K)} \right] $

$ \begin{array}{l} S(H) = \frac{{4\sqrt M }}{{15{\alpha ^2}}}\left[ {\left( {3MN - 2{M^2} - {N^2}} \right){E_1}(K) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {2\left( {{M^2} + {N^2} - 2MN} \right){E_2}(K)} \right] \end{array} $

$ \begin{array}{l} Q(H) = \frac{{4\sqrt M }}{{35{\alpha ^2}}}\left[ {\left( {8{N^3} - 13M{N^2} + 3{M^2}N + 2{M^3}} \right){E_1}(K) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {2\left( { - 8{N^3} + 12M{N^2} - 2{M^2}N - {M^3}} \right){E_2}(K)} \right] \end{array} $

$ M = \alpha {A^2} - k,N = \frac{{\alpha {A^2}}}{2},K = \frac{N}{M} $

where E₁(K), E₂(K) are the first and second kind of complete elliptic integrals, respectively, and

$ {E_1}(K) = \int_0^{\frac{{\rm{ \mathsf{ π} }}}{2}} {\frac{1}{{\sqrt {1 - K{{\sin }^2}\theta } }}{\rm{d}}\theta } $

$ {E_2}(K) = \int_0^{\frac{\pi }{2}} {\sqrt {1 - K{{\sin }^2}\theta } {\rm{d}}\theta } $

where A is the positive root of equation U(A) =H. Thus, combination with Eq.(11) obtains the following equation:

$ A = \sqrt {\frac{{\sqrt {{k^2} + 4\alpha H} + k}}{\alpha }} $

(17)

The average FPK equation of Eq.(14) is

$ \frac{{\partial p}}{{\partial t}} = - \frac{\partial }{{\partial H}}\left[ {\langle m(H)\rangle p} \right] + \frac{1}{2}\frac{{{\partial ^2}}}{{\partial H}}\left[ {\left\langle {{\sigma ^2}(H)} \right\rangle p} \right] $

(18)

where p represents the probability density function (PDF) of energy function H.

The steady state solution implies $ \frac{{\partial p}}{{\partial t}} = 0$ in Eq.(18), and the corresponding boundary conditions is taken as

$ \left\{ {\begin{array}{*{20}{l}} {0 \le p < \infty ,\quad H = 0}\\ {p \to 0,\frac{{{\rm{d}}p}}{{{\rm{d}}H}} \to 0,H \to \infty } \end{array}} \right. $

(19)

Therefore, the steady state solution of slow variable H can be described as

$ p(H) = {C_0}\exp [ - \lambda (H)] $

(20)

where

$ \begin{array}{l} \lambda (H) = - 2\int_0^H {\frac{{ - \mu B(x) + \gamma C(x) - \beta Q(x) + (R - 1)x}}{{B(x)}}{\rm{d}}x} - \\ \;\;\;\;\;\;\;\;\;\;\ln \frac{{{T_{\frac{1}{4}}}(H)}}{{{T_{\frac{1}{4}}}(0)}} \end{array} $

(21)

and C₀ is a normalized constant. The joint PDF of Eq.(5) can be obtained as below:

$ \begin{array}{l} p\left( {{x_1},{x_2}} \right) = {p_{{y_1}{y_2}}}\left( {{y_1},{y_2}} \right)\left| {\frac{{\partial \left( {{y_1},{y_2}} \right)}}{{\partial \left( {{x_1},{x_2}} \right)}}} \right| = \\ \;\;\;\;\;\;\;{{\tilde p}_{{y_1}{y_2}}}\left( {{x_1},{x_2}} \right) + {{\tilde p}_{{y_1}{y_2}}}\left( { - {x_1}, - {x_2}} \right) = \\ \;\;\;\;\;\;\;2{{\tilde p}_{{y_1}{y_2}}}\left( {{x_1},{x_2}} \right) = {\left. {2\frac{{p(H)}}{{T(H)}}} \right|_{H = \frac{{y_2^2}}{2} - \frac{k}{2}y_1^2 + \alpha \frac{{y_1^4}}{4}}} \cdot \\ \;\;\;\;\;\;\;\left| {\frac{{\partial \left( {{y_1},{y_2}} \right)}}{{\partial \left( {{x_1},{x_2}} \right)}}} \right|,{x_1} \ge 0 \end{array} $

(22)

The edge PDFs are received from Eq.(22):

$ p\left( {{x_1}} \right) = \int_{ - \infty }^{ + \infty } p \left( {{x_1},{x_2}} \right){\rm{d}}{x_2} $

(23)

$ p\left( {{x_2}} \right) = \int_{ - \infty }^{ + \infty } p \left( {{x_1},{x_2}} \right){\rm{d}}{x_1} $

(24)

4 Analysis of Stochastic Response

In this section, the system parameters are taken as k=－1.0, α=0.01, γ=0.01, β=0.01. Based on the original Eq.(5), the numerical results can be obtained by Monte Carlo method. Moreover, stochastic bifurcations is discussed.

4.1 Numerical Simulation

The influence of restitution coefficient R and noise amplitude f on the probability responses of Eq.(5) is considered. Let the parameter μ be 0.01 in this subsection.

As illustrated in Fig. 3, it is obvious that analytical results agree with numerical results. When noise amplitude is fixed, the smaller restitution coefficient R could lead to higher peaks of PDFs p(H), displacement p(x₁), and velocity p(x₂). In addition, Fig. 3(a), (b), and (c) show that the influence of noise amplitude f is similar to restitution coefficient R. These results indicate that the smaller restitution value can make the larger loss of energy, and stochastic perturbation can result in the system's deviation from the equilibrium position.

4.2 Stochastic Bifurcations

In this subsection, the noise amplitude f is fixed as 0.1. The effect of linear damping coefficient μ and restitution coefficient R on PDFs of Eq.(5) were explored.

It can be seen from Figs. 4-5, analytical results agree well with numerical results when μ=－0.01. When R=0.98, a single peak appears in the joint PDFs as shown in Fig. 4. As the restitution coefficient R increases to 1.0, the join PDFs transits to the shape of crater. This change of PDFs implies that Eq.(5) undergoes a P-bifurcation.

To discover the detailed process of P-bifurcation, the curve of the edge PDFs for different restitution coefficient R is shown in Fig. 6.

Obviously, there is only a singular peak for R=0.98. As restitution coefficient R increases, the curve becomes flatter and eventually turns into a crater for R=1.0. It is can be seen in Fig. 6 that the critical value is about R≈0.993 during transition. On the other hand, the evolution of edge PDFs implies that a smaller R can cause a larger loss of energy.

Comparison of Fig. 7 and Fig. 8 shows that analytical results coincide with numerical results well. As μ decreases from -0.01 to -0.02, the curve changes from single peak to a crater. The changes indicate the occurrence of P-bifurcation.

For different linear damping coefficients μ, Fig. 9 describes the procedure of P-bifurcation in detail. It is easy to see one peak occurs for μ=－0.01. As μ decreases, the curve gets increasingly flatter. Finally, when μ=－0.02, the form of crater appears. From the evolution of edge PDFs, it was found that μ≈－0.018 was the critical value that made stochastic bifurcation happen.

5 Conclusions

In this paper, the probability response of SMA impact model under external Gaussian white noise excitation is investigated. First, the non-smooth transformation was employed to deal with the discontinuous position. Then the original system turned into an approximate system associated with the Dirac function. Lastly, the stochastic response of the approximate system was obtained after analysis. Meanwhile, comparison of analytical results with numerical results demonstrates the effectiveness of the stochastic averaging method. It reveals that the stationary probability response of system is affected by the increase of noise amplitude and restitution force.

Furthermore, this paper explores the P-bifurcation of SMA beam system. It is worth noting that the restitution coefficient R and linear damping coefficient μ can cause P-bifurcation. However, the P-bifurcation usually results in transition of dynamical system, in which unstable vibration is generated and the stability of the SMA beam is destroyed. Therefore, by adjusting the system parameters, the security of the SMA beam could be enhanced.