黏弹形状记忆合金系统的随机响应分析

Stochastic Response Analysis of Viscoelastic Shape Memory Alloy System

摘要

本文研究了窄带噪声激励下黏弹形状记忆合金（Shape Memory Alloy, SMA）弹簧振子系统的动力学响应。根据SMA多项式本构关系，引入黏弹性介质描述弹簧阻尼，建立了随机SMA动力学方程。采用Lindstedt-Poincaré和多尺度相结合方法，得到系统的幅频响应和稳态矩方程，并探讨了幅频响应与稳态矩随黏弹参数、噪声强度等变化的规律。研究结果表明，黏弹项能够影响系统的非线性刚度和结构阻尼；噪声强度和黏弹参数都能够诱导系统发生随机跃迁现象。这些都为黏弹SMA复合结构的稳定性提供了理论参考。

关键词

窄带噪声 / 形状记忆合金 / 多尺度法 / 随机跃迁

Key words

narrow-band noise / shape memory alloy / multi-scale method / random transition

0 引言

形状记忆合金（Shape Memory Alloy， SMA）作为一类特殊的具有功能性和结构强度的新型智能材料，以其独特的形状记忆、超弹性、电阻敏感性和优良的抗腐蚀性等特性被广泛关注^［1-2］。尤其是超弹性特性，也称伪弹性，是指在外力作用下产生的应变可以随着外力的卸载而恢复。SMA这种优异的力学性能主要是由于高温下的马氏体和低温下的奥氏体在温度和压力作用下发生了相变。近年来，SMA在土木工程、航空航天、医疗设备等领域广泛应用^［3-6］。

随着智能材料的快速发展，SMA结构的动力学特性，如混沌、分岔、激变等已成为非线性动力学问题中的一个重要研究方向，主要包括形状记忆合金弹簧振子、形状记忆合金梁、形状记忆合金支架等等。例如，Piccirillo等^［7-8］研究了SMA弹簧振子的主共振和次共振，通过相位图、功率谱图、李雅普诺夫指数和庞加莱图像分析了自由振动和受迫振动时振子表现出的混沌行为，发现了丰富的随机动力学现象。Lacarbonara等^［9］基于热力学模型研究了单自由度SMA振子的非线性响应和分岔现象。Spanos等^［10］通过Preisach方法研究了白噪声激励下SMA振荡器的响应。Dobson等^［11］研究了具有不对称迟滞的单自由度SMA振荡器系统的随机振动，该系统由Bouc-Wen模型描述。Yan和Nie^［12］采用等效线性化技术来预测受到高斯白噪声作用的SMA振荡器的伪弹性响应。Yue等^［13-14］通过数值胞映射方法从全局角度研究了SMA振子系统的动力学行为。Guo等^［15-16］利用多尺度方法和随机平均法分别研究了窄带噪声和高斯色噪声激励下SMA弹簧振子系统的响应问题。Zhang等^［17］研究了SMA层合梁系统的厚度比、激励参数对系统分岔特性的影响。

近年来，黏弹性材料由于既能存储势能又能耗散能量的特性，已被广泛用到力学和工程系统中^［18-20］。Zhu和Cai^［21］利用随机平均法研究了宽带噪声激励下黏弹系统的随机振动，并讨论了黏弹参数对系统的影响。Ling等^［22］等利用随机平均法研究了单自由度强非线性黏弹系统的响应及其稳定性。Deng^［23］应用随机平均法和摄动法计算了矩Lyapunov指数，研究了宽带噪声激励下线性黏弹系统的随机稳定性。Xu等^［24-25］将Lindstedt-Poincaré摄动法与多尺度法相结合，研究了随机激励下含分数阶阻尼的黏弹系统的动力学响应。孙双双等^［26］利用有限元法对黏弹基形状记忆合金复合梁的弯曲问题进行了建模和分析，实验结果表明热黏弹性理论可以用来分析SMA复合材料的响应问题。秦惠增等^［27］研究了具有SMA层的黏弹柔性梁的动力学响应特性，分析了应力诱发马氏体相变的影响。目前，基于黏弹性材料和SMA相结合的复合结构的研究大都是针对确定性系统，随机情形的很少。而在实际工程中，不确定的随机因素会对结构的振动、稳定性造成灾难性破坏。因此，研究随机噪声激励下黏弹SMA系统的动力学问题具有重要的实际应用价值。本文利用SMA多项式本构关系，建立了窄带噪声激励下的黏弹性阻尼SMA弹簧振子振动方程，研究了系统主共振时的稳态响应，分析了黏弹性参数、外部激励幅值、噪声强度对系统稳态响应的影响，为黏弹性SMA复合结构的应用提供了理论基础。

1 SMA弹簧振子模型

考虑如图1具有黏弹性阻尼的单自由度SMA弹簧振子系统，系统由质量块M、黏弹性阻尼器C和SMA弹簧组成。

采用Falk^［28］提出的多项式本构方程描述SMA振子的弹簧恢复力，在Guo等^［15］建立的SMA振动模型的基础上，加入黏弹项 H \left(x,l,t\right)，引入无量纲变换，得到系统振动的无量纲方程为：

\ddot{x} + {\epsilon }^{\mathrm{2}} \alpha \dot{x} + {\epsilon }^{\mathrm{2}} \beta H \left(x,l,t\right) + {\omega }_{\mathrm{0}}^{\mathrm{2}} x + \epsilon \gamma x^{\mathrm{3}} + {\epsilon }^{\mathrm{2}} \delta x^{\mathrm{5}} = {\epsilon }^{\mathrm{2}} f \xi \left(t\right) , （1）

其中黏弹项 H \left(x,l,t\right)为

H \left(x,l,t\right) = {\int }_{\mathrm{0}}^t \frac{\mathrm{1}}{l}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{t-s}{l}\right) x ( s ) \mathrm{d} s , （2）

其中 x为无量纲位移，t是无量纲时间， \alpha， \beta， l及 {\omega }_{\mathrm{0}}分别为系统阻尼、黏弹项系数、黏弹参数及系统固有频率， \gamma和 \delta是非线性项系数，f为参激的激励强度， \epsilon为小参数。系统所受外部激励 \xi \left(t\right)为窄带噪声，具体形式为

\xi \left(t\right) = \mathrm{c} \mathrm{o} \mathrm{s} \left[\Omega t+hW\left(t\right)\right] , （3）

其中 \Omega是随机激励的频率，h是标准维纳过程 W \left(t\right)的强度。

2 窄带噪声激励下的理论分析

在这一部分我们把Lindstedt-Poincaré方法和多尺度方法相结合，研究系统的近似解析解。

首先，根据Lindstedt-Poincaré方法，引入时间变换，将系统的振幅和频率联系起来，即

\tau = \omega t，（4）

则黏弹项变形为

\begin{array}{c}H\left(x,l,t\right)={\int }_{\mathrm{0}}^t\frac{\mathrm{1}}{l}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{t-s}{l}\right)x(s)\mathrm{d}s=\\ {\int }_{\mathrm{0}}^{\tau }\frac{\mathrm{1}}{l}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{\mathrm{1}}{l}\left(\frac{\tau }{\omega }-\frac{\eta }{\omega }\right)\right]x(\eta )\mathrm{d}\frac{\eta }{\omega }={\int }_{\mathrm{0}}^{\tau }\frac{\mathrm{1}}{l\omega }\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\tau -\eta }{l\omega }\right)x(\eta )\mathrm{d}\eta =H(\overline{l},x,\tau ),\end{array} （5）

其中 \overline{l} = l \omega。因此把式（3） -式（5）代入式（1），得到新的系统方程：

{\omega }^{\mathrm{2}} x ″ + {\epsilon }^{\mathrm{2}} \omega \alpha x \text{'} + {\epsilon }^{\mathrm{2}} \beta H ( \overline{l} , x , \tau ) + {\omega }_{\mathrm{0}}^{\mathrm{2}} x + \epsilon \gamma x^{\mathrm{3}} + {\epsilon }^{\mathrm{2}} \delta x^{\mathrm{5}} = {\epsilon }^{\mathrm{2}} f \mathrm{c} \mathrm{o} \mathrm{s} \left[\frac{\Omega \tau }{\omega }+hW\left(\frac{\tau }{\omega }\right)\right] , （6）

其中 x ″代表对时间变量 \tau的二阶导数。

其次，引入快与慢时间尺度， T_{\mathrm{0}} = \tau， T_{\mathrm{1}} = \epsilon \tau ， T_{\mathrm{1}} = {\epsilon }^{\mathrm{2}} \tau，令

D_{\mathrm{0}} = \frac{\partial }{\partial T_{\mathrm{0}}} , D_{\mathrm{1}} = \frac{\partial }{\partial T_{\mathrm{1}}} , D_{\mathrm{2}} = \frac{\partial }{\partial T_{\mathrm{2}}} ,

则常微分算子可以展开为下列偏微分序列

\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}\tau }=D_{\mathrm{0}}+\epsilon D_{\mathrm{1}}+{\epsilon }^{\mathrm{2}}{D}_{\mathrm{2}}+\cdots ,\\ \frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{\tau }^{\mathrm{2}}}=D_{\mathrm{0}}^{\mathrm{2}}+\mathrm{2}\epsilon D_{\mathrm{0}}{D}_{\mathrm{1}}+{\epsilon }^{\mathrm{2}}\left(D_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{D}_{\mathrm{0}}{D}_{\mathrm{2}}\right)+\cdots \mathrm{。}\end{array} （7）

进而，方程（6）的近似解析解可表示为如下形式

x = x_{\mathrm{0}} \left(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}}\right) + \epsilon x_{\mathrm{1}} \left(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}}\right) + {\epsilon }^{\mathrm{2}} x_{\mathrm{2}} \left(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}}\right) \mathrm{。} （8）

同时将频率以 {\omega }_{\mathrm{0}}^{\mathrm{2}} = {\omega }^{\mathrm{2}} - \epsilon {\omega }_{\mathrm{1}} - {\epsilon }^{\mathrm{2}} {\omega }_{\mathrm{2}}的形式进行摄动，将式（7）和式（8）代入式（6）中，重新整理，根据小参数 \epsilon的阶次，得到一系列的摄动方程

O \left({\epsilon }^{\mathrm{0}}\right) : {\omega }^{\mathrm{2}} D_{\mathrm{0}}^{\mathrm{2}} x_{\mathrm{0}} + {\omega }^{\mathrm{2}} x_{\mathrm{0}} = \mathrm{0} , （9a）

O \left({\epsilon }^{\mathrm{1}}\right) : {\omega }^{\mathrm{2}} D_{\mathrm{0}}^{\mathrm{2}} x_{\mathrm{1}} + {\omega }^{\mathrm{2}} x_{\mathrm{1}} = - \mathrm{2} {\omega }^{\mathrm{2}} D_{\mathrm{0}} D_{\mathrm{1}} x_{\mathrm{0}} + {\omega }_{\mathrm{1}} x_{\mathrm{0}} - \gamma x_{\mathrm{0}}^{\mathrm{3}} , （9b）

\begin{array}{c}O\left({\epsilon }^{\mathrm{2}}\right):{\omega }^{\mathrm{2}}{D}_{\mathrm{0}}^{\mathrm{2}}{x}_{\mathrm{2}}+{\omega }^{\mathrm{2}}{x}_{\mathrm{2}}=-\mathrm{2}{\omega }^{\mathrm{2}}{D}_{\mathrm{0}}{D}_{\mathrm{1}}{x}_{\mathrm{1}}-{\omega }^{\mathrm{2}}\left(D_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{D}_{\mathrm{0}}{D}_{\mathrm{2}}\right)x_{\mathrm{0}}-\omega \alpha D_{\mathrm{0}}{x}_{\mathrm{0}}-\beta H(\overline{l},x,\tau )+\\ {\omega }_{\mathrm{1}}{x}_{\mathrm{1}}+{\omega }_{\mathrm{2}}{x}_{\mathrm{0}}-\mathrm{3}\gamma x_{\mathrm{0}}^{\mathrm{2}}{x}_{\mathrm{1}}-\delta x_{\mathrm{0}}^{\mathrm{5}}+f\mathrm{c}\mathrm{o}\mathrm{s}\left[\frac{\Omega \tau }{\omega }+\overline{h}W\left(\tau \right)\right],\end{array} （9c）

其中 \overline{h} = \frac{h}{\sqrt[]{\omega }}。

由微分方程知识可设（9a）的解为 x_{\mathrm{0}} = a \mathrm{c} \mathrm{o} \mathrm{s} \left(T_{\mathrm{0}}+\phi \right)，等价形式为

x_{\mathrm{0}} = a \left(T_{\mathrm{1}},T_{\mathrm{2}}\right) {\mathrm{e}}^{\mathrm{i}{T}_{\mathrm{0}}} + \overline{A} \left(T_{\mathrm{1}},T_{\mathrm{2}}\right) {\mathrm{e}}^{-\mathrm{i}{T}_{\mathrm{0}}}，（10）

其中 \overline{A}为 A的共轭，显然我们可以得到 A = \frac{\mathrm{1}}{\mathrm{2}} a {\mathrm{e}}^{\mathrm{i}\phi }。

将式（10）代入方程（9b）中，消除久期项，得到

- \mathrm{2} \mathrm{i} {\omega }^{\mathrm{2}} D_{\mathrm{1}} A + {\omega }_{\mathrm{1}} A - \mathrm{3} \gamma A^{\mathrm{2}} \overline{A} = \mathrm{0} \mathrm{。} （11）

对于方程（11）来说，我们假设 D_{\mathrm{1}} A = \mathrm{0}，可得 {\omega }_{\mathrm{1}} = \mathrm{3} \gamma A \overline{A}。显然 {\omega }_{\mathrm{1}}是一个实数，也就是说我们的假设是合理的。将式（10）和（11）代入方程（9b）中，消除久期项，求解方程可得

x_{\mathrm{1}} = \frac{\gamma }{\mathrm{8}{\omega }^{\mathrm{2}}} A^{\mathrm{3}} {\mathrm{e}}^{\mathrm{3}\mathrm{i}{T}_{\mathrm{0}}} + c . c .，（12）

其中 c . c .表示其前一项的共轭复数。因此我们得到了系统的一阶近似解，为了进一步研究系统的主共振响应，将参激频率展开为如下形式

\Omega = \omega \left(\mathrm{1}+{\epsilon }^{\mathrm{2}}\sigma \right)，（13）

根据黏弹项表达式（5），忽略瞬态部分可得

\begin{array}{c}H\left(\overline{l},x_{\mathrm{0}},T_{\mathrm{0}}\right)=\frac{\mathrm{1}}{\overline{l}}{\int }_{\mathrm{0}}^{T_{\mathrm{0}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{T_{\mathrm{0}}-s}{\overline{l}}\right)x_{\mathrm{0}}(s)\mathrm{d}s=\frac{A}{\overline{l}}{\mathrm{e}}^{-\frac{T_{\mathrm{0}}}{\overline{l}}}{\int }_{\mathrm{0}}^{T_{\mathrm{0}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{1}}{\overline{l}}+\mathrm{i}\right)s\mathrm{d}s+c.c.=\\ A{\mathrm{e}}^{-\frac{T_{\mathrm{0}}}{\overline{l}}}\frac{\mathrm{1}-\overline{l}\mathrm{i}}{\mathrm{1}+{\overline{l}}^{\mathrm{2}}}\left({\mathrm{e}}^{\left(\frac{\mathrm{1}}{\overline{l}}+\mathrm{i}\right)T_{\mathrm{0}}}-\mathrm{1}\right)+c.c.=\frac{\mathrm{1}-\overline{l}\mathrm{i}}{\mathrm{1}+{\overline{l}}^{\mathrm{2}}}{x}_{\mathrm{0}}-\frac{\mathrm{1}-\overline{l}\mathrm{i}}{\mathrm{1}+{\overline{l}}^{\mathrm{2}}}A{\mathrm{e}}^{-\frac{T_{\mathrm{0}}}{\overline{l}}}+c.c.\mathrm{。}\end{array} （14）

将式（12） -（14）代入方程（9c）中，消除久期项可得

- \mathrm{2} {\omega }^{\mathrm{2}} \mathrm{i} D_{\mathrm{2}} A - \omega \alpha \mathrm{i} A - \frac{\mathrm{1}-\overline{l}\mathrm{i}}{\mathrm{1}+{\overline{l}}^{\mathrm{2}}} \beta A + {\omega }_{\mathrm{2}} A - \frac{\mathrm{3}{\gamma }^{\mathrm{2}}}{\mathrm{8}{\omega }^{\mathrm{2}}} A^{\mathrm{3}} {\overline{A}}^{\mathrm{2}} - \mathrm{10} \delta A^{\mathrm{3}} {\overline{A}}^{\mathrm{2}} + \frac{\mathrm{1}}{\mathrm{2}} f {\mathrm{e}}^{\mathrm{i}(\sigma T_{\mathrm{2}}+\tilde{h}W(T_{\mathrm{2}}))} = \mathrm{0} , （15）

其中

\overline{h} W ( T_{\mathrm{0}} ) = \overline{h} W \left(\frac{T_{\mathrm{2}}}{{\epsilon }^{\mathrm{2}}}\right) = \left(\frac{\overline{h}}{\epsilon }\right) W \left(T_{\mathrm{2}}\right) = \tilde{h} W ( T_{\mathrm{2}} ) \mathrm{。}

对于式（15），如果假设 D_{\mathrm{2}} A = \mathrm{0}，得到 {\omega }_{\mathrm{2}}是一个复数，说明假设不合理，所以令 {\omega }_{\mathrm{2}} = \mathrm{0}，那么方程（15）变为

- \mathrm{i} {\omega }^{\mathrm{2}} \left(D_{\mathrm{2}}a+\mathrm{i}a{D}_{\mathrm{2}}\phi \right) - \frac{\mathrm{1}}{\mathrm{2}} \omega \alpha \mathrm{i} a - \frac{\mathrm{1}-\overline{l}\mathrm{i}}{\mathrm{2}(\mathrm{1}+{\overline{l}}_{\mathrm{2}})} \beta a - \frac{\mathrm{3}{\gamma }^{\mathrm{2}}}{\mathrm{256}{\omega }^{\mathrm{2}}} a^{\mathrm{5}} - \frac{\mathrm{5}\delta }{\mathrm{16}} a^{\mathrm{5}} - \frac{\mathrm{1}}{\mathrm{2}} f {\mathrm{e}}^{\mathrm{i}\psi } = \mathrm{0} , （16）

其中 \psi = \sigma T_{\mathrm{2}} + \tilde{h} W \left(T_{\mathrm{2}}\right) - \phi。

根据以上分析可得系统频率 \omega近似为： \omega = \sqrt[]{{\omega }_{\mathrm{0}}^{\mathrm{2}}+\frac{\mathrm{3}}{\mathrm{4}}\epsilon \gamma a^{\mathrm{2}}}。

分离方程（16）的实部和虚部得到

\left\{\begin{array}{l}{D}_{\mathrm{2}}a=-\frac{\mathrm{1}}{\mathrm{2}\omega }\alpha a+\frac{\overline{l}}{\mathrm{2}{\omega }^{\mathrm{2}}\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta a+\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}}f\mathrm{s}\mathrm{i}\mathrm{n}\psi ,\\ a{D}_{\mathrm{2}}\psi =a\left(\sigma +\tilde{h}{D}_{\mathrm{2}}W\right)-\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta a-\frac{\mathrm{3}{\gamma }^{\mathrm{2}}}{\mathrm{256}{\omega }^{\mathrm{4}}}{a}^{\mathrm{5}}-\frac{\mathrm{5}\delta }{\mathrm{16}{\omega }^{\mathrm{2}}}{a}^{\mathrm{5}}+\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}}f\mathrm{c}\mathrm{o}\mathrm{s}\psi \mathrm{。}\end{array}\right. （17）

通过求解公式（17）可以得到系统的振幅 a和相位 \psi。因此，结合公式（10），公式（12）和公式（17），我们得到了系统（6）的二阶近似解析解

x = a \mathrm{c} \mathrm{o} \mathrm{s} \left(T_{\mathrm{0}}+\phi \right) + \epsilon \frac{\delta }{\mathrm{32}{\omega }^{\mathrm{2}}} a^{\mathrm{3}} \mathrm{c} \mathrm{o} \mathrm{s} \mathrm{3} \left(T_{\mathrm{0}}+\phi \right) + O \left({\epsilon }^{\mathrm{2}}\right) , （18）

其中 O \left(\epsilon \right)代表小参数 \epsilon的高阶项。

首先考虑确定情形下 h = \mathrm{0}的SMA弹簧振子系统。当 D_{\mathrm{2}} a = D_{\mathrm{2}} \psi = \mathrm{0}时，公式（17）可化简为

\left\{\begin{array}{l}-\frac{\mathrm{1}}{\omega }\alpha a_{\mathrm{0}}+\frac{\overline{l}}{{\omega }^{\mathrm{2}}\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta a_{\mathrm{0}}+\frac{\mathrm{1}}{{\omega }^{\mathrm{2}}}f\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_{\mathrm{0}}=\mathrm{0},\\ \sigma -\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta -\frac{\mathrm{3}{\gamma }^{\mathrm{2}}}{\mathrm{256}{\omega }^{\mathrm{4}}}{a}_{\mathrm{0}}^{\mathrm{4}}-\frac{\mathrm{5}\delta }{\mathrm{16}{\omega }^{\mathrm{2}}}{a}_{\mathrm{0}}^{\mathrm{4}}+\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}{a}_{\mathrm{0}}}f\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_{\mathrm{0}}=\mathrm{0}\mathrm{。}\end{array}\right. （19）

消除式（19）中的参数 {\psi }_{\mathrm{0}}，得到黏弹性阻尼SMA模型的幅频响应方程为

{\left[\frac{\omega \alpha }{f}{a}_{\mathrm{0}}-\frac{\overline{l}}{f\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta a_{\mathrm{0}}\right]}^{\mathrm{2}} + \left[-\frac{\mathrm{2}{\omega }^{\mathrm{2}}{a}_{\mathrm{0}}}{f}\sigma +\frac{\mathrm{1}}{f\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta a_{\mathrm{0}}+\right. {\left.\frac{\mathrm{3}{\gamma }^{\mathrm{2}}}{\mathrm{128}{\omega }^{\mathrm{2}}f}{a}_{\mathrm{0}}^{\mathrm{5}}+\frac{\mathrm{5}\delta }{\mathrm{8}f}{a}_{\mathrm{0}}^{\mathrm{5}}\right]}^{\mathrm{2}} = \mathrm{1} \mathrm{。} （20）

进一步考虑噪声 h \ne \mathrm{0}对SMA弹簧振子系统的影响。基于摄动思想，引入如下展开式

a = a_{\mathrm{0}} + a_{\mathrm{1}} , \psi = {\psi }_{\mathrm{0}} + {\psi }_{\mathrm{1}} , （21）

其中a_{\mathrm{1}}和{\psi }_{\mathrm{1}}分别是确定性结果a_{\mathrm{0}}和{\psi }_{\mathrm{0}}

的小摄动，可以被视为小量。将式（21）代入式（17）

，\left\{\begin{array}{l}{D}_{\mathrm{2}}{a}_{\mathrm{1}}=-\frac{\mathrm{1}}{\mathrm{2}\omega }\alpha a_{\mathrm{1}}+\frac{\overline{l}}{\mathrm{2}{\omega }^{\mathrm{2}}\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta a_{\mathrm{1}}+\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}}f{\psi }_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_{\mathrm{0}},\\ D_{\mathrm{2}}{\psi }_{\mathrm{1}}=\tilde{h}{D}_{\mathrm{2}}W\left(T_{\mathrm{1}}\right)-\frac{\mathrm{3}{\gamma }^{\mathrm{2}}}{\mathrm{64}{\omega }^{\mathrm{4}}}{a}_{\mathrm{0}}^{\mathrm{3}}{a}_{\mathrm{1}}-\frac{\mathrm{5}\delta }{\mathrm{4}{\omega }^{\mathrm{2}}}{a}_{\mathrm{0}}^{\mathrm{3}}{a}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}{a}_{\mathrm{0}}}f{\psi }_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_{\mathrm{0}}-\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}{a}_{\mathrm{0}}^{\mathrm{2}}}f{a}_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_{\mathrm{0}}\mathrm{。}\end{array}\right. （22）

结合式（21）消除高阶小量可以得到如下It \widehat{\mathrm{o}}方程

\left\{\begin{array}{l}\mathrm{d}{a}_{\mathrm{1}}=\left(-\frac{\mathrm{1}}{\mathrm{2}\omega }\alpha a_{\mathrm{1}}\right.+\frac{\overline{l}}{\mathrm{2}{\omega }^{\mathrm{2}}\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta a_{\mathrm{1}}+\left.\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}}f{\psi }_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_{\mathrm{0}}\right)\mathrm{d}{T}_{\mathrm{2}},\\ \mathrm{d}{\psi }_{\mathrm{1}}=\left(-\frac{\mathrm{3}{\gamma }^{\mathrm{2}}}{\mathrm{64}{\omega }^{\mathrm{4}}}\right.a_{\mathrm{0}}^{\mathrm{3}}{a}_{\mathrm{1}}-\frac{\mathrm{5}\delta }{\mathrm{4}{\omega }^{\mathrm{2}}}{a}_{\mathrm{0}}^{\mathrm{3}}{a}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}{a}_{\mathrm{0}}}f{\psi }_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_{\mathrm{0}}-\left.\frac{\mathrm{1}}{\mathrm{2}{\omega }^{\mathrm{2}}{a}_{\mathrm{0}}^{\mathrm{2}}}f{a}_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_{\mathrm{0}}\right)\mathrm{d}{T}_{\mathrm{2}}+\tilde{h}\mathrm{d}W\left(T_{\mathrm{2}}\right)\mathrm{。}\end{array}\right. （23）

在稳态情形下，根据It \widehat{\mathrm{o}}微积分的性质有

\frac{\mathrm{d}E\left(a_{\mathrm{1}}\right)}{\mathrm{d}{T}_{\mathrm{2}}} = \frac{\mathrm{d}E\left({\psi }_{\mathrm{1}}\right)}{\mathrm{d}{T}_{\mathrm{2}}} = \frac{\mathrm{d}E\left(a_{\mathrm{1}}^{\mathrm{2}}\right)}{\mathrm{d}{T}_{\mathrm{2}}} = \frac{\mathrm{d}E\left({\psi }_{\mathrm{1}}^{\mathrm{2}}\right)}{\mathrm{d}{T}_{\mathrm{2}}} = \frac{\mathrm{d}E\left(a_{\mathrm{1}}{\psi }_{\mathrm{1}}\right)}{\mathrm{d}{T}_{\mathrm{2}}} = \mathrm{0} \mathrm{。} （24）

其中 E \left(\cdot \right)表示数学期望算子。则可以得到系统振幅 a的一阶和二阶稳态矩

E \left(a\right) = a_{\mathrm{0}} + E \left(a_{\mathrm{1}}\right) = a_{\mathrm{0}} , （25a）

E \left(a^{\mathrm{2}}\right) = a_{\mathrm{0}}^{\mathrm{2}} + E \left(a_{\mathrm{1}}^{\mathrm{2}}\right) = a_{\mathrm{0}}^{\mathrm{2}} + \frac{{\tilde{h}}^{\mathrm{2}}{P}_{\mathrm{2}}^{\mathrm{2}}}{P_{\mathrm{1}}{P}_{\mathrm{2}}{Q}_{\mathrm{1}}+P_{\mathrm{1}}{Q}_{\mathrm{2}}^{\mathrm{2}}+P_{\mathrm{1}}^{\mathrm{2}}{Q}_{\mathrm{2}}+Q_{\mathrm{1}}{Q}_{\mathrm{2}}{P}_{\mathrm{2}}} , （25b）

其中

\begin{array}{c}{P}_{\mathrm{1}}=\frac{\mathrm{1}}{\omega }\alpha -\frac{\overline{l}}{{\omega }^{\mathrm{2}}\left(\mathrm{1}+{\overline{l}}^{\mathrm{2}}\right)}\beta ,P_{\mathrm{2}}=\frac{\mathrm{1}}{{\omega }^{\mathrm{2}}}f\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_{\mathrm{0}},\\ Q_{\mathrm{1}}=\frac{\mathrm{3}{\gamma }^{\mathrm{2}}}{\mathrm{32}{\omega }^{\mathrm{4}}}{a}_{\mathrm{2}}^{\mathrm{3}}+\frac{\mathrm{5}\delta }{\mathrm{2}{\omega }^{\mathrm{2}}}{a}_{\mathrm{0}}^{\mathrm{3}}+\frac{\mathrm{1}}{{\omega }^{\mathrm{2}}{a}_{\mathrm{0}}^{\mathrm{2}}}f\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_{\mathrm{0}},\\ Q_{\mathrm{2}}=\frac{\mathrm{1}}{{\omega }^{\mathrm{2}}{a}_{\mathrm{0}}^{\mathrm{2}}}f\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_{\mathrm{0}}\mathrm{。}\end{array}

3 随机响应分析

这一部分，我们将通过数值模拟验证所得近似解析解的正确性，进而论证理论分析方法的有效性。在所有数值模拟中，主要系统参数根据文献［26］给出： \alpha = \mathrm{0.2}， \epsilon = \sqrt[]{\mathrm{0.1}} , \beta = \mathrm{0.5}， {\omega }_{\mathrm{0}} = \mathrm{2.0}， \gamma = \mathrm{87.4}， \delta = \mathrm{2.8} \times {\mathrm{10}}^{\mathrm{4}}，其他参数的取值会在对应位置具体说明。

如图2所示，关于确定性系统的幅频响应，通过将经典的林德斯泰特-庞加莱（Lindstedt-Poincaré，L-P）与多尺度方法结合，获得的近似解析解结果和数值结果具有很好的一致性。从图2（a）与图2（b）可以看出，当外部激励的频率 \Omega或幅值 f在一定的范围内时，系统呈现出多值响应现象，如图中蓝色区域所示，在这个范围内，黏弹性SMA系统的响应包含高振幅振荡和低振幅振荡，即系统出现了双稳态行为，这很大程度上依赖于系统的初值条件。这种现象通过图3就可以很直观地说明。如图3（a），当初始条件为 x \left(\mathrm{0}\right) = \mathrm{0.01} , \dot{x} \left(\mathrm{0}\right) = \mathrm{0.01}时，系统处于低振幅振荡状态，对应于双稳态区域幅值响应的下分支。图3（b）中，当初始条件为 x \left(\mathrm{0}\right) = \mathrm{0.1} , \dot{x} \left(\mathrm{0}\right) = \mathrm{0.1}时，系统处于高振幅振荡状态，对应于双稳态区域振幅响应的上分支。此外，在幅频响应中，随着外部激励频率 \Omega或幅值 f的增大或减小，可以观察到系统响应从上半支/下半支到下半支/上半支的随机跃迁现象。

接下来，我们研究在确定情形下不同外部激励幅值 f及黏弹性参数 l对SMA弹簧系统稳态响应，结果如图4所示。从图4（a）可以看到，随着外部激励幅值 f的增大，双稳区域向右移动，对应区间的面积也增大。同样，随着黏弹性参数 l的增大，双稳区域面积也逐渐增大，如图4（b）所示。这主要是因为增大黏弹性参数 l降低了系统的结构阻尼，提高了SMA系统的结构刚度。

现在，我们对具有非零噪声强度 h \ne \mathrm{0}的黏弹性SMA系统的稳态响应进行详细分析。前面利用L-P与多尺度相结合的方法，得到了随机激励下系统的解析一阶矩（23a）和二阶矩（23b）。利用四阶龙格库塔算法，从图5（a）和图5（b）可以看出理论解和数值解吻合得很好，验证了本方法的有效性。

此外，如图5所示，随着外部激励频率 \Omega在一定范围内增大或减小时，黏弹性SMA弹簧振子系统出现了多值区域，同确定性情形一样。图6（a）， 6（b）研究了不同黏弹性参数 l下系统的稳态一阶矩和稳态二阶矩，我们发现随着黏弹性参数 l的逐渐增大，双稳区域有向右移动，对应的区域面积增大。根据式（14）可知，随着黏弹性参数 l \to \mathrm{0}，黏弹性阻尼退化为经典的纯黏性阻尼，在这种情况下双稳现象消失，说明了在随机情形下，黏弹性参数 l可影响系统的双稳特性。

接下来，我们采用Monte Carlo数值模拟计算了响应幅值的稳态概率密度。如图7所示，在随机扰动作用下，我们发现稳态概率密度的拓扑结构发生了显著变化。当噪声强度 h = \mathrm{0.01}时，观察到一个低振幅振荡的峰。随着噪声强度 h的逐渐增大，稳态概率密度出现双峰。当噪声强度 h足够大时，稳态概率密度又一次变为单峰。从稳态密度的变化可以猜想到，黏弹性SMA系统的稳态响应包含两种可能的运动形式，即高幅值和低幅值的随机振动，这意味着此刻系统发生了随机跃迁现象。

最后，为了更直观地展示图5所示的跃迁现象，绘制了一些时间历程图来说明。在图8（a）中，当 h = \mathrm{0.09}，即噪声强度较小时，系统响应仅受到轻微扰动，且没有观察到随机跃迁现象。随着噪声强度的增大，如图8（b） -8（d）所示，系统响应出现了从高振幅振荡到低振幅振荡的明显跃迁现象，并且这种现象变得越来越频繁。

以上分析表明，窄带噪声会诱导黏弹性阻尼SMA弹簧振子系统出现随机跃迁现象，这种稳态概率密度拓扑结构变化的现象也称为随机P分岔。需要注意的是，频繁的随机跃迁会造成SMA系统内部的损伤，进而造成安全隐患。因此，在实际工程应用中，设计人员应采取相应的措施来避免这种情况出现。

4 结论

本文考虑材料的黏弹性，建立了随机激励下黏弹SMA弹簧振子模型。利用理论方法得到系统的二阶近似解析解，并通过数值模拟验证了该方法的有效性。当外部激励频率或幅值取值在一定范围内时，系统出现双稳行为；随着激励幅值、黏弹参数及黏弹系数的逐渐增大，振幅的稳态矩和双稳区域的频率区间也随之增大；窄带噪声能够诱导SMA系统出现高幅值和低幅值振荡状态之间的随机跃迁现象。随机跃迁的发生意味着系统稳定性的改变，为了提升系统的稳定性与安全性，可以从SMA模型响应的早期预警方面进一步研究。

参考文献

[1] KUMBHAR S B, CHAVAN S P, GAWADE S S. Adaptive Tuned Vibration Absorber Based on Magnetorheological Elastomer-shape Memory Alloy Composite[J]. Mech Syst Signal Process, 2018, 100: 208-223. DOI: 10.1016/j.ymssp.2017.07.027 .

[2] COSTA D D A, SAVI M A, DE PAULA A S, et al. Chaos Control of a Shape Memory Alloy Structure Using Thermal Constrained Actuation[J]. Int J Non Linear Mech, 2019, 111: 106-118. DOI: 10.1016/j.ijnonlinmec.2019.02.006 .

[3] CHO H K, RHEE J. Nonlinear Finite Element Analysis of Shape Memory Alloy (SMA) Wire Reinforced Hybrid Laminate Composite Shells[J]. Int J Non Linear Mech, 2012, 47(6): 672-678. DOI: 10.1016/j.ijnonlinmec.2011.11.002 .

[4] FONSECA L M, RODRIGUES G V, SAVI M A, et al. Nonlinear Dynamics of an Origami Wheel with Shape Memory Alloy Actuators[J]. Chaos Solitons Fractals, 2019, 122: 245-261. DOI: 10.1016/j.chaos.2019.03.033 .

[5] KONH B, HONARVAR M, HUTAPEA P. Design Optimization Study of a Shape Memory Alloy Active Needle for Biomedical Applications[J]. Med Eng Phys, 2015, 37(5): 469-477. DOI: 10.1016/j.medengphy.2015.02.013 .

[6] COSTA D D A, SAVI M A. Nonlinear Dynamics of an SMA-pendulum System[J]. Nonlinear Dyn, 2017, 87(3): 1617-1627. DOI: 10.1007/s11071-016-3137-y .

[7] PICCIRILLO V, BALTHAZAR J M, PONTES B R. Analytical Study of the Nonlinear Behavior of a Shape Memory Oscillator: Part I: Primary Resonance and Free Response at Low Temperatures[J]. Nonlinear Dyn, 2010, 59(4): 733-746. DOI: 10.1007/s11071-009-9573-1 .

[8] PICCIRILLO V, BALTHAZAR J M, PONTES B R. Analytical Study of the Nonlinear Behavior of a Shape Memory Oscillator: Part II: Resonance Secondary[J]. Nonlinear Dyn, 2010, 60(4): 513-524. DOI: 10.1007/s11071-009-9611-z .

[9] LACARBONARA W, BERNARDINI D, VESTRONI F. Nonlinear Thermomechanical Oscillations of Shape-memory Devices[J]. Int J Solids Struct, 2004, 41(5/6): 1209-1234. DOI: 10.1016/j.ijsolstr.2003.10.015 .

[10] SPANOS P D, CACCIOLA P, REDHORSE J. Random Vibration of SMA Systems via Preisach Formalism[J]. Nonlinear Dyn, 2004, 36(2): 405-419. DOI: 10.1023/B: NODY.0000045514.54248.fa .

[11] DOBSON S, NOORI M, HOU Z, et al. Modeling and Random Vibration Analysis of SDOF Systems with Asymmetric Hysteresis[J]. Int J Non Linear Mech, 1997, 32(4): 669-680. DOI: 10.1016/s0020-7462(96)00090-x .

[12] YAN X, NIE J. Response of Sma Superelastic Systems under Random Excitation[J]. J Sound Vib, 2000, 238(5): 893-901. DOI: 10.1006/jsvi.2000.3020 .

[13] YUE X L, XIANG Y L, XU Y, et al. Global Dynamics of the Dry Friction Oscillator with Shape Memory Alloy[J]. Arch Appl Mech, 2020, 90(12): 2681-2692. DOI: 10.1007/s00419-020-01741-z .

[14] YUE X L, XIANG Y L, ZHANG Y, et al. Global Analysis of Stochastic Bifurcation in Shape Memory Alloy Supporter with the Extended Composite Cell Coordinate System Method[J]. Chaos Interdiscip J Nonlinear Sci, 2021, 31(1): 013133. DOI: 10.1063/5.0024992 .

[15] GUO R, LIU Q, LI J L, et al. Response Statistics of a Shape Memory Alloy Oscillator with Random Excitation[J]. Appl Sci, 2021, 11(21): 10175. DOI: 10.3390/app112110175 .

[16] BU M, LI J, GUO R, et al. Dynamical Responses of a Gaussian Colored Noise-driven Shape Memory Alloy Oscillator with a Periodic Force[J]. Chaos Solitons Fractals, 2023, 166: 112966. DOI: 10.1016/j.chaos.2022.112966 .

[17] ZHANG X H, Wu Z Q. Bifurcation Analysis of Shape Memory Alloy Laminated Beam[J]. Chin J Appl Mech, 2017, 34: 397-403. DOI: 10.11776/cjam.34.03.C023 .

[18] DRESSLER M. Inertial Theories for Dilute Viscoelastic Polymer Blends with a Volume Preserving Microstructure[J]. J Non Newton Fluid Mech, 2012, 173/174: 40-48. DOI: 10.1016/j.jnnfm.2012.02.002 .

[19] REY A D, HERRERA-VALENCIA E E. Liquid Crystal Models of Biological Materials and Silk Spinning[J]. Biopolymers, 2012, 97(6): 374-396. DOI: 10.1002/bip.21723 .

[20] WANG X, SCHOEN J A, RENTSCHLER M E. A Quantitative Comparison of Soft Tissue Compressive Viscoelastic Model Accuracy[J]. J Mech Behav Biomed Mater, 2013, 20: 126-136. DOI: 10.1016/j.jmbbm.2013.01.007 .[PubMed]

[21] ZHU W Q, CAI G Q. Random Vibration of Viscoelastic System Under Broad-band Excitations[J]. Int J Non Linear Mech, 2011, 46(5): 720-726. DOI: 10.1016/j.ijnonlinmec.2011.02.004 .

[22] LING Q, JIN X L, HUANG Z L. Response and Stability of SDOF Viscoelastic System Under Wideband Noise Excitations[J]. J Frankl Inst, 2011, 348(8): 2026-2043. DOI: 10.1016/j.jfranklin.2011.05.019 .

[23] DENG J. Higher-order Stochastic Averaging for a SDOF Fractional Viscoelastic System Under Bounded Noise Excitation[J]. J Frankl Inst, 2017, 354(17): 7917-7945. DOI: 10.1016/j.jfranklin.2017.09.019 .

[24] XU Y, LI Y G, LIU D. Response of Fractional Oscillators with Viscoelastic Term Under Random Excitation[J]. J Comput Nonlinear Dyn, 2014, 9(3): 031015. DOI: 10.1115/1.4026068 .

[25] LIU Q, XU Y, KURTHS J. Bistability and Stochastic Jumps in an Airfoil System with Viscoelastic Material Property and Random Fluctuations[J]. Commun Nonlinear Sci Numer Simul, 2020, 84: 105184. DOI: 10.1016/j.cnsns.2020.105184 .

[26] 孙双双, 孙国钧, 杜晓伟, 黏弹基形状记忆合金复合梁的弯曲[J]. 上海交通大学学报, 2002, 36(11): 1663-1666. DOI: 10.3321/j.issn: 1006-2467.2002.11.030 .

[27] SUN S S, SUN G J, DU X W, et al. Bending of Viscoelastic Beam with Embedded Shape Memory Alloy Wires[J]. J Shanghai Jiaotong Univ, 2002, 36(11): 1663-1666. DOI: 10.3321/j.issn: 1006-2467.2002.11.030 .

[28] 秦惠增, 任勇生. 含形状记忆合金(SMA)层柔性梁的振动频响特性[J]. 机械强度, 2002, 24(1): 45-48. DOI: 10.3321/j.issn:1001-9669.2002.01.013 .

[29] QIN H Z, REN Y S. Frequency Response Characteristics of Vibration of Flexible Beam Covered with Shape Memory Alloy(Sma) Layers[J]. J Mech Strength, 2002, 24(1): 45-48. DOI: 10.3321/j.issn:1001-9669.2002.01.013 .

[30] FALK F. Model Free Energy, Mechanics, and Thermodynamics of Shape Memory Alloys[J]. Acta Metall, 1980, 28(12): 1773-1780. DOI: 10.1016/0001-6160(80)90030-9 .

基金资助

山西省高等学校科技创新项目(2022L316)