一类随机边界刚性约束悬臂梁系统周期运动的稳定性分析

Stability Analysis of Periodic Motion for a Class of Cantilever Beam System with Rigid and Random Constraints

摘要: 本文研究了一类二阶可微随机约束下的碰撞悬臂梁系统周期解的稳定性。通过推导含参数的随机零时间不连续映射给出相应的跳跃矩阵，结合跳跃矩阵和光滑流映射的基解矩阵得到了随机线性化矩阵。基于随机线性化矩阵探讨了随机约束对系统稳定性的影响，进一步调查了周期解失稳之后的倍化分岔现象，数值仿真验证了理论的有效性。

Abstract: This paper studies the stability of periodic solution for a class of second-order differentiable cantilever beam system with rigid and random constraints. The saltatory matrix is presented by deducing the random zero-time discontinuous mapping with parameters. The random linearization matrix is obtained by combining the saltatory matrix and the fundamental solution matrix of smooth flow mapping. Based on the random linearization matrix, the influence of stochastic constraints on the system stability is discussed, and the doubling bifurcation phenomenon after the instability of the periodic solution is further investigated. The validity of the theory is verified by numerical simulations.

1. 引言

悬臂梁结构广泛存在于机械工程、民用建筑、航空航天等领域，工程中许多结构都可以简化成附加集中质量块的悬臂梁模型。在实际工程零件中，零件内部或零件之间往往有较大的间隙，间隙的存在会导致复杂的动力学行为，研究这类振动问题对于认识系统的故障机理及特征、减小系统振动和冲击，以及保障安全运行都有重要的意义。因此，如何优化机械动力系统的设计，减少带间隙系统的部件损耗是实际工程中的一个重要课题。

在以往的研究中人们通常考虑的是理想模型，然而自然界和生活中很多现象都是随机的，系统受到随机噪声的扰动是难以避免的，在模型求解中加入随机因素则更符合客观规律。近年来，关于这方面的研究逐渐成为一个新兴的领域。碰撞振动系统在不同随机噪声激励下的随机响应近年来受到了广泛的关注并取得了很大的进展。目前的文献大都是利用一些随机的方法将随机激励下的碰撞振动系统转化为确定性系统来研究其随机响应 [1] - [6]，其中应用最广泛的是随机平均法 [3] [4] [5] [6] [7] 和正交多项式逼近法 [8]。Simpson等 [9] 在分段光滑系统中增加了一类小振幅的高斯噪声，导出了该噪声的Nordmark映射，研究了其中参数对系统擦边分岔的影响，并探讨了随机动力学和滑动近周期轨道的随机平均法 [10]。此外，他根据振动冲击系统的噪声源推导出了三种随机Nordmark映射，并研究了它们对擦边分岔的影响 [11]。Staunton等考虑了间隙中的随机性，将随机过程引入分段光滑系统和混合系统的不连续边界，分别推导出相应的随机跳跃矩阵和随机不连续映射，并讨论了理论的实用性 [12] [13]。

本文以Ornstein-Uhlenbeck (以下简称OU)过程为基本过程，该随机过程已经广泛应用在生物医学 [14]、金融领域 [15]、风速模型、电路、机器学习和控制论等领域。根据文献 [12] [13] 的方法，推导出2阶可微的随机过程，以及相应的随机速度和随机加速度，其优势是在理论推导时随机因素也相应的可以求导。本文以文献 [16] 中描述的两自由度悬臂梁冲击振动系统在随机约束下的力学模型为研究对象，研究了当系统轨迹与碰撞面横截相交时系统周期运动的稳定性。首先，建立了系统在随机约束下的无量纲动力学方程，通过解耦得到系统的周期解，并建立横截周期条件，得到了周期解的初始相位角和积分常数。其次，利用定相位的庞加莱截面，导出了含参数随机零时间不连续映射、跳跃矩阵和线性化矩阵。最后，利用推导的理论研究了随机约束系统的周期运动稳定性及分岔问题。

2. 随机边界约束两自由度悬臂梁系统的动力学方程及横截周期运动

2.1. 力学模型和动力学方程

图1所示的力学模型是一类随机约束下悬臂梁系统的两自由度简化模型，两个集中质量块质量分别为 $m_1$ 和 $m_2$，质块 $m_1$ 到固定端的距离为l，质块 $m_2$ 到固定端的距离为 $l_2$，抗弯刚度为 $EI$。约束边界 $\tilde{\Sigma }$ 是在原确定约束面 $\tilde{\Sigma }$ 基础上添加一类随机过程而构成的一个随机约束面。作用在质块 $m_1$ 上的周期激振力为 $e{\eta }^2\sin \left(\eta \tilde{t}+{\tau }_0\right)$。

Figure 1. Stochastic constrained two-degree-of-freedom cantilever beam system

图1. 随机边界刚性约束两自由度悬臂梁系统模型

在任意相邻两次碰撞 ${\tilde{x}}_1-\tilde{d}<{\tilde{\chi }}_2\left(\tilde{t}\right)$ 之间，系统的动力学方程为：

$\left[\begin{array}{cc}{m}_1& 0\\ 0& m_2\end{array}\right]\{\begin{array}{c}{\ddot{\stackrel{˜}{x}}}_1\\ {\ddot{\stackrel{˜}{x}}}_2\end{array}\}+\upsilon \left[\begin{array}{cc}{m}_1& 0\\ 0& m_2\end{array}\right]\{\begin{array}{c}{\stackrel{\dot{}}{\tilde{x}}}_1\\ {\stackrel{\dot{}}{\tilde{x}}}_2\end{array}\}+\left[\begin{array}{cc}K{\tilde{B}}_{11}& K{\tilde{B}}_{12}\\ K{\tilde{B}}_{21}& K{\tilde{B}}_{22}\end{array}\right]\{\begin{array}{c}{\tilde{x}}_1\\ {\tilde{x}}_2\end{array}\}=\{\begin{array}{c}e{\eta }^2\\ 0\end{array}\}\sin \left(\eta \tilde{t}+{\tau }_0\right)$ (1)

这里， ${\tilde{\chi }}_2\left(\tilde{t}\right)$ 表示一类二阶可微的随机过程，其推导将在第3节中给出。

发生碰撞 ${\tilde{x}}_1-\tilde{d}={\tilde{\chi }}_2\left(\tilde{t}\right)$ 时，根据碰撞法则，质块 $m_1$ 的碰撞方程为：

${\stackrel{\dot{}}{\tilde{x}}}_{1+}=-r{\stackrel{\dot{}}{\tilde{x}}}_{1-}$ (2)

其中， ${\stackrel{\dot{}}{\tilde{x}}}_{1-}$ 和 ${\stackrel{\dot{}}{\tilde{x}}}_{1+}$ 表示质块与约束碰撞前后的瞬时速度， $\eta$ 是激振频率， $\upsilon$ 是阻尼系数，r为碰撞恢复系数。

刚度矩阵形式为 $\left[\begin{array}{cc}K{\tilde{B}}_{11}& K{\tilde{B}}_{12}\\ K{\tilde{B}}_{21}& K{\tilde{B}}_{22}\end{array}\right]=\frac{1}{W}\left[\begin{array}{cc}12EI{l}_2^3& -6EI{l}_2^2\left(3l-l_2\right)\\ -6EI{l}_2^2\left(3l-l_2\right)& 12EI{l}^3\end{array}\right]$，其中 $W=4l^3{l}_2^3-l_2^4{\left(3l-l_2\right)}^2$。

引入无量纲量：

$m=\frac{m_2}{m_1}$， $x_i={\tilde{x}}_i\cdot \frac{m_1}{e}$， $\xi =\frac{\upsilon }{2}\sqrt{\frac{m_1}{K{\tilde{B}}_{11}}}$， $t=\tilde{t}\cdot \sqrt{\frac{K{\tilde{B}}_{11}}{m_1}}$， $K_{ij}=\frac{K{\tilde{B}}_{ij}}{K{\tilde{B}}_{11}}$， $\omega =\eta \cdot \sqrt{\frac{m_1}{K{\tilde{B}}_{11}}}$， $d=\tilde{d}\frac{m_1}{e}$， $\left(i,j=1,2\right)$， ${\chi }_2\left(t\right)={\tilde{\chi }}_2\left(\tilde{t}\right)\frac{m_1}{e}$。

使用上面的无量纲量可将系统(1)和(2)转化为如下的无量纲形式：

$\left[\begin{array}{cc}1& 0\\ 0& m\end{array}\right]\{\begin{array}{c}{\ddot{x}}_1\\ {\ddot{x}}_2\end{array}\}+2\xi \left[\begin{array}{cc}1& 0\\ 0& m\end{array}\right]\{\begin{array}{c}{\stackrel{\dot{}}{x}}_1\\ {\stackrel{\dot{}}{x}}_2\end{array}\}+\left[\begin{array}{cc}1& K_{12}\\ K_{21}& K_{22}\end{array}\right]\{\begin{array}{c}{x}_1\\ x_2\end{array}\}=\{\begin{array}{c}{\omega }^2\\ 0\end{array}\}\sin \left(\omega t+{\tau }_0\right)$， $x_1-d<{\chi }_2\left(t\right)$ (3)

${\stackrel{\dot{}}{x}}_{1+}=-r{\stackrel{\dot{}}{x}}_{1-}$， $x_1-d={\chi }_2\left(t\right)$ (4)

用 $\Psi$ 表示(3)式的正则模态矩阵， ${\omega }_1$ 和 ${\omega }_2$ 是系统未发生碰撞时的固有频率。取 $\Psi$ 为转换矩阵，令 $\stackrel{⌢}{X}={\left(x_1,x_2\right)}^{\text{T}}$ 和 $Y={\left(y_1,y_2\right)}^{\text{T}}$，在变换 $\stackrel{⌢}{X}=\Psi Y$ 下，方程(3)可变为

$I\ddot{Y}+C\stackrel{\dot{}}{Y}+\Lambda Y=P\sin \left(\omega t+{\tau }_0\right)$ (5)

其中，I表示 $2\times 2$ 的单位矩阵， $C=2\xi I$， $\Lambda =diag\left[{\omega }_1^2,{\omega }_2^2\right]$， $P={\Psi }^{\text{T}}{\left({\omega }^2,0\right)}^{\text{T}}$。

根据方程(5)可得方程(3)的通解为：

$x_i\left(t\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}{\psi }_{ij}}\left({\text{e}}^{-{\eta }_{j}t}\left(a_j\cos \left({\omega }_{dj}t\right)+b_j\sin \left({\omega }_{dj}t\right)\right)+A_j\sin \left(\omega t+{\tau }_0\right)+B_j\cos \left(\omega t+{\tau }_0\right)\right)$， $i=1,2$ (6)

${\stackrel{\dot{}}{x}}_i\left(t\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}{\psi }_{ij}}\left({\text{e}}^{-{\eta }_{j}t}\left(\left(-{\eta }_j{b}_j-{\omega }_{dj}{a}_j\right)\sin \left({\omega }_{dj}t\right)+\left(-{\eta }_j{a}_j+{\omega }_{dj}{b}_j\right)\cos \left({\omega }_{dj}t\right)\right)+A_j\omega \cos \left(\omega t+{\tau }_0\right)-B_j\omega \sin \left(\omega t+{\tau }_0\right)\right)$， $i=1,2$ (7)

其中 ${\psi }_{ij}$ 为正则模态矩阵 $\Psi$ 中对应的元素， ${\eta }_j=\xi$， ${\omega }_{dj}=\sqrt{{\omega }_j^2-{\xi }^2}$， $a_j$ 和 $b_j$ 是积分常数 $A_j$ 和 $B_j$ 是振幅常数，表达式如下：

$A_j=\frac{{\omega }_j^2-{\omega }^2}{{\left({\omega }_j^2-{\omega }^2\right)}^2+{\left(2\xi \omega \right)}^2}{P}_j$， $B_j=-\frac{2\xi \omega }{{\left({\omega }_j^2-{\omega }^2\right)}^2+{\left(2\xi \omega \right)}^2}{P}_j$。

2.2. 横截周期运动

如图1所示，以质块 $m_1$ 的静平衡位置为坐标原点，垂直向下为正方向， $t=0$ 时刻从约束面出发。本文考虑的随机边界的变化在 ${10}^{-3}~{10}^{-4}$，极其微小，因此随机系统的横截周期条件以固定约束面定义如下：

$x_1\left(0\right)=x_1\left(T\right)=d$， ${\stackrel{\dot{}}{x}}_1^{+}\left(0\right)=-r{\stackrel{\dot{}}{x}}_1^{-}\left(T\right)$， $x_2\left(0\right)=x_2\left(T\right)$， ${\stackrel{\dot{}}{x}}_2\left(0\right)={\stackrel{\dot{}}{x}}_2\left(T\right)$ (8)

将(6)和(7)代入(8)中，可以求得初相位角和四个积分常数的表达式如下：

$M_j=\left(\frac{e_j{\omega }_{dj}\left(C_j-e_j\right)}{1-e_j{C}_j}\right)-{\omega }_{dj}$

$N_j=\frac{\left(R+1\right)e_j{\eta }_j{S}_j-R{\omega }_{dj}{e}_j\left(C_j-e_j\right)}{1-e_j{C}_j}-{\omega }_{dj}$

$h=\left(\frac{{\psi }_{22}{M}_2\left(R+1\right)\omega }{{\psi }_{11}{N}_1{\psi }_{22}{M}_2-{\psi }_{12}{N}_2{\psi }_{21}{M}_1}\right)\left(\frac{{\psi }_{11}{e}_1{S}_1}{1-e_1{C}_1}-\frac{{\psi }_{12}{e}_2{S}_2}{1-e_2{C}_2}\cdot \frac{{\psi }_{21}{M}_1}{{\psi }_{22}{M}_2}\right)$

其中 $S_j=\sin \left({\omega }_{dj}T\right)$， $C_j=\cos \left({\omega }_{dj}T\right)$， $e_j={\text{e}}^{-{\eta }_{j}T}\left(j=1,2\right)$。

当 $d=0$ 时：

${\bar{\tau }}_0=\arctan \left(\frac{hp+q}{hq-p}\right)$

当 $d\ne 0$ 时：

${\tau }_0=\arccos \left(\frac{\bar{d}\tan {\bar{\tau }}_0\pm \sqrt{{\tan }^2{\bar{\tau }}_0-{\bar{d}}^2+1}}{{\tan }^2{\bar{\tau }}_0+1}\right)$

$b_1=\frac{{\psi }_{22}{M}_2\left(R+1\right)\omega \left(p\cos {\tau }_0-q\sin {\tau }_0\right)}{{\psi }_{11}{N}_1{\psi }_{22}{M}_2-{\psi }_{12}{N}_2{\psi }_{21}{M}_1}$

$b_2=-\frac{{\psi }_{21}{M}_1\left(R+1\right)\omega \left(p\cos {\tau }_0-q\sin {\tau }_0\right)}{{\psi }_{11}{N}_1{\psi }_{22}{M}_2-{\psi }_{12}{N}_2{\psi }_{21}{M}_1}$

$a_j=\frac{e_j{b}_j{S}_j}{1-e_j{C}_j}$

其中 $p={\psi }_{11}{A}_1+{\psi }_{12}{A}_2$， $q={\psi }_{11}{B}_1+{\psi }_{12}{B}_2$， $\bar{d}=\frac{d}{hq-p}$。

3. 一类二阶可微的随机过程

Ornstein-Uhlenbeck (OU)过程的随机微分方程 [13] [17] 形式如下：

$\text{d}\chi (t)=\theta \left(\mu -\chi \left(t\right)\right)\text{d}t+\sigma \text{d}{B}_t$， $\chi \left(0\right)={\chi }_0$ (9)

其中 $\theta ,\sigma \in \mathrm{R}$ 是正参数， $\theta$ 表示系统对干扰的反映程度， $\sigma$ 表示噪声的大小， $\mu$ 为均值， $B_t$ 为维纳过程，也称布朗运动，是一种外界的随机噪声，其中，维纳过程利用高斯过程实现，每一段时间间隔内的增量服从高斯分布。

(3)式是一个一阶线性微分方程，解为Itó积分的形式为：

$\chi \left(t\right)=\mu +\left({\chi }_0-\mu \right){\text{e}}^{-\theta t}+\sigma {\displaystyle {\int }_0^t{\text{e}}^{\theta \left(s-t\right)}}\text{d}{B}_s$

通过计算得 $\chi \left(t\right)$ 的期望是：

$E\left(\chi \left(t\right)\right)=\mu +\left({\chi }_0-\mu \right){\text{e}}^{-\theta t}$

$\chi \left(t\right)$ 的方差是：

$Var\left(\chi \left(t\right)\right)=\frac{{\sigma }^2}{2\theta }\left(1-{\text{e}}^{-2\theta t}\right)$

下面利用OU过程作为基本过程，按照随机微分方程(9)的结构 [13] 首先构造1阶可微的随机微分方程：

$\frac{\text{d}{\chi }_1\left(t\right)}{\text{d}t}=\theta \left(\mu -{\chi }_1\left(t\right)\right)+\sigma \chi (\; t\; )$

其Itó形式解为：

${\chi }_1\left(t\right)=\mu +\frac{\mu }{\theta }\sigma +{\text{e}}^{-\theta t}\left(\left(\sigma {\xi }_0-\sigma \mu \right)t+C_1\right)+{\sigma }^2{\displaystyle {\int }_0^t{\text{e}}^{\theta \left(s-t\right)}\left(t-s\right)\text{d}{B}_s}$

期望和方差分别为：

$E\left({\chi }_1\left(t\right)\right)=\mu \left(1+\frac{\sigma }{\theta }\right)+{\text{e}}^{-\theta t}\left(\left({\chi }_0-\mu \right)\sigma t+C_1\right)$

$Var\left({\chi }_1\left(t\right)\right)=\frac{{\sigma }^4}{2\theta }\left(\frac{2}{{\left(2\theta \right)}^2}-{\text{e}}^{-2\theta t}\left(t^2+\frac{t}{\theta }+\frac{1}{2{\theta }^2}\right)\right)$

其中 $C_1={\chi }_1\left(0\right)-\mu \left(1+\frac{\sigma }{\theta }\right)$。

其次，构造2阶可微的随机微分方程：

$\frac{\text{d}{\chi }_2\left(t\right)}{\text{d}t}=\theta \left(\mu -{\chi }_2\left(t\right)\right)+\sigma {\chi }_1(\; t\; )$

其Itó形式解为：

${\chi }_2\left(t\right)=\mu +\mu \frac{\sigma }{\theta }\left(1+\frac{\sigma }{\theta }\right)+\left(\frac{{\left(\sigma t\right)}^2}{2}\left({\chi }_0-\mu \right)+C_1\sigma t+C_2\right){\text{e}}^{-\theta t}+{\frac{\sigma }{2}}^3{\displaystyle {\int }_0^t{\text{e}}^{\theta \left(s-t\right)}{\left(t-s\right)}^2}\text{d}{B}_s$

其中 $C_2={\chi }_2\left(0\right)-\mu \left(1+\frac{\sigma }{\theta }\left(1+\frac{\sigma }{\theta }\right)\right)$。

期望和方差分别为：

$E\left({\chi }_2\left(t\right)\right)=\mu +\mu \frac{\sigma }{\theta }\left(1+\frac{\sigma }{\theta }\right)+\left(\frac{{\left(\sigma t\right)}^2}{2}\left({\chi }_0-\mu \right)+C_1\sigma t+C_2\right){\text{e}}^{-\theta t}$

$Var\left({\chi }_2\left(t\right)\right)=\frac{{\sigma }^6}{8\theta }\left(\frac{4!}{{\left(2\theta \right)}^4}-{\text{e}}^{-2\theta t}{\displaystyle \underset{i=0}{\overset{4}{\sum }}\frac{4!}{\left(4-i\right)!}\frac{t^{4-i}}{{\left(2\theta \right)}^i}}\right)$

最后，基于 $\chi \left(t\right)$ 、 ${\chi }_1\left(t\right)$ 和 ${\chi }_2\left(t\right)$ 构造随机速度和随机加速度：

$V_2\left(t\right)=\frac{\text{d}{\chi }_2\left(t\right)}{\text{d}t}=\theta \left(\mu -{\chi }_2\left(t\right)\right)+\sigma {\chi }_1(\; t\; )$

$A_2\left(t\right)=\frac{\text{d}{V}_2\left(t\right)}{\text{d}t}=\theta \mu \left(1-\theta +\sigma \right)+{\theta }^2{\chi }_2\left(t\right)-2\theta \sigma {\chi }_1\left(t\right)+{\sigma }^2\chi (\; t\; )$

Itó形式解分别为：

$V_2\left(t\right)=\theta \left(\mu -E\left({\chi }_2\left(t\right)\right)\right)+\sigma E\left({\chi }_1\left(t\right)\right)+{\sigma }^3{\displaystyle {\int }_0^t{\mathrm{e}}^{\theta \left(s-t\right)}\left(\left(t-s\right)-\theta \frac{{\left(t-s\right)}^2}{2}\right)\text{d}{B}_s}$

$\begin{array}{c}{A}_2\left(t\right)=\theta \mu \left(1-\theta +\sigma \right)+{\theta }^{2}E\left({\gamma }_2\left(t\right)\right)-2\theta \sigma E\left({\chi }_1\left(t\right)\right)+{\sigma }^{2}E\left({\chi }_0\left(t\right)\right)\\ +{\sigma }^3{\displaystyle {\int }_0^t{\text{e}}^{\theta \left(s-t\right)}\left(1-\theta \left(t-s\right)-\theta \left(t-s\right)+\frac{{\theta }^2}{2}{\left(t-s\right)}^2\right)\text{d}{B}_s}\end{array}$

相应的期望和方差分别为：

$E\left(V_2\left(t\right)\right)=\theta \left(\mu -E\left({\chi }_2\left(t\right)\right)\right)+\sigma E\left({\chi }_1(\; t\; )\right)$

$Var\left(V_2\left(t\right)\right)={\theta }^{2}Var\left({\chi }_2\left(t\right)\right)+{\sigma }^{2}Var\left({\chi }_1\left(t\right)\right)-2\theta \sigma Cov\left({\chi }_2\left(t\right),{\chi }_1(\; t\; )\right)$

$E\left(A_2\left(t\right)\right)={\theta }^{2}E\left({\chi }_2\left(t\right)\right)-2\theta \sigma E\left({\chi }_1\left(t\right)\right)+{\sigma }^{2}E\left(\chi (\; t\; )\right)$

$\begin{array}{c}Var\left(V_2\left(t\right)\right)={\theta }^{4}Var\left({\chi }_2\left(t\right)\right)+5{\theta }^2{\sigma }^{2}Var\left({\chi }_1\left(t\right)\right)+{\sigma }^{4}Var\left(\chi \left(t\right)\right)\\ -4{\theta }^3\sigma Cov\left({\chi }_2\left(t\right),{\chi }_1\left(t\right)\right)-4\theta {\sigma }^{3}Cov({\chi }_1\left(t\right),\chi (\; t\; ))\end{array}$

式中

$Cov\left({\chi }_2\left(t\right),{\chi }_1\left(t\right)\right)=\frac{{\sigma }^5}{4\theta }\left(\frac{3}{4{\theta }^3}-{\mathrm{e}}^{-2\theta t}{\displaystyle \underset{i=0}{\overset{3}{\sum }}\frac{6}{\left(3-i\right)!}}\frac{t^{3-i}}{{\left(2\theta \right)}^i}\right)$

$Cov\left({\chi }_1\left(t\right),\chi \left(t\right)\right)=\frac{{\sigma }^3}{2\theta }\left(1-t{\mathrm{e}}^{-2\theta t}-\frac{{\mathrm{e}}^{-2\theta t}}{2\theta }\right)$

4. 随机线性化矩阵

令 $X={\left(x_1,{\stackrel{\dot{}}{x}}_1,x_2,{\stackrel{\dot{}}{x}}_2,s\right)}^{\text{T}}$， $t=s$，方程(3)转化为如下自治形式：

$\stackrel{\dot{}}{X}=F\left(X,\mu \right)$ (10)

其中 $F\left(X,\mu \right)=\left(\begin{array}{c}{v}_1\\ -2\xi v_1-x_1-k_{12}{x}_2+{\omega }^2\sin \left(\omega s+{\tau }_0\right)\\ v_2\\ \frac{-2\xi m{v}_2-k_{21}{x}_1-k_{22}{x}_2}{m}\\ 1\end{array}\right)$， $\mu$ 为间隙参数d。

确定性系统的边界函数为 $h\left(X,\mu \right)=X-d$，在确定约束面处添加随机过程 ${\chi }_2\left(t\right)$，本文随机过程的均值为0，则随机约束的边界函数为 $\tilde{h}\left(X,\mu ,t\right)=h\left(X,\mu \right)-{\chi }_2\left(t\right)$，相应地随机约束面定义为：

$\tilde{\Sigma }=\left\{X:\tilde{h}\left(X,\mu ,t\right)=0\right\}$ (11)

如图2，选择Poincaré截面为 ${\Pi }_s=\left\{\left(x_1,v_1,x_2,v_2,\theta \right)\in {\text{R}}^4\times S|\theta =0\mathrm{mod}\left(\frac{2\pi }{\omega }\right)\right\}$。首先考虑周期轨线，设 $\bar{j}\left(t\right)$ 是与碰撞面横截相交的周期轨线，从不动点 ${\bar{X}}_0\in {\Pi }_S$ 出发到达点 ${\bar{X}}_1^{*}\in \Sigma$ 需要时间 $t_1$，然后由碰撞恢复函数 $R\left(X,\mu \right)$ 将 ${\bar{X}}_1^{*}$ 点映射到 ${\bar{X}}_2^{*}=R\left({\bar{X}}_1^{*},\mu \right)$ 点，再由 ${\bar{X}}_2^{*}$ 点经过时间 $t_2$ 返回 ${\bar{X}}_0$，这里 $t_1+t_2=T$。

Figure 2. Schematic diagram of local discontinuous mapping

图2. 横截周期轨道附近局部不连续映射示意图

下面考虑扰动轨线，记为 $j\left(t\right)$，从 ${\bar{X}}_0$ 附近一点 $X_0$ 出发经过固定时间 $t_1$ 后并没有到达或者已经超越约束面 $\Sigma$，将这点记为X，如图2所示。然后根据局部不连续映射，到点 $X_3=D\left(X,\mu \right)$，再经过时间 $t_2$ 到达点 $X_T$。从点X到约束面 $\Sigma$ 所需要的时间记为 $\delta$ (此处 $\delta$ 也可能是负的)，从 $\Sigma$ 到达随机约束面 $\tilde{\Sigma }$ 所需要的时间记为 $\Delta t$。

首先，将函数 $\tilde{h}\left(\varphi \left({\bar{X}}_0,\mu ,t\right),\mu ,t\right)$ 在 $t=t_1$ 处关于t泰勒展开到一阶

$\begin{array}{l}\tilde{h}\left(\varphi \left({\bar{X}}_0,\mu ,t\right),\mu ,t\right)\\ =h\left(\varphi \left({\bar{X}}_0,\mu ,t\right),\mu \right)-{\chi }_2\left(t\right)\\ =h\left(\varphi \left({\bar{X}}_0,\mu ,t_1\right),\mu \right)-{\chi }_2\left(t_1\right)+\left(h_X\left(\varphi \left({\bar{X}}_0,\mu ,t_1\right),\mu \right)F_1^{*}-V\left(t_1\right)\right)\left(t-t_1\right)+h.o.t\\ =-{\chi }_2\left(t_1\right)+\left(h_X\left({\bar{X}}_1^{*},\mu \right)F_1^{*}-V\left(t_1\right)\right)\left(t-t_1\right)+h.o.t\end{array}$ (12)

式中 ${\bar{X}}_1^{*}=\varphi \left({\bar{X}}_0,\mu ,t_1\right)$， $F_1^{*}$ 为点 ${\bar{X}}_1^{*}$ 处的向量场。在(12)式中省去高阶项，可得

$\Delta t=\frac{{\chi }_2\left(t_1\right)}{h_X\left({\bar{X}}_1^{*},\mu \right)F_1^{*}-V\left(t_1\right)}$ (13)

下面推导随机线性化矩阵。从 $X_0$ 出发，经过一个周期内的映射为：

$\varphi \left(X_0,\mu ,T\right)=\varphi \left(R\left(\varphi \left(X_0,\mu ,t\right),\mu ,t\right),\mu ,T-t\right)$ (14)

此时时间t不再固定，而是一个随机变量，用 $t=t\left(X_0,\mu \right)$ 表示从点 $X_0$ 到随机约束面 $\tilde{\Sigma }$ 所需要的时间，有如下关系式：

$t\left(X_0,\mu \right)=t_1+\delta +\Delta t$ (15)

根据时间关系式(15)，(14)式可以写为：

$\varphi \left(X_0,\mu ,T\right)=\varphi \left(\varphi \left(\varphi \left(R\left(\varphi \left(\varphi \left(\varphi \left(X_0,\mu ,t_1\right),\mu ,\delta \right),\mu ,\Delta t\right),\mu ,t\right),\mu ,-\Delta t\right),\mu ,-\delta \right),\mu ,T-t_1\right)$ (16)

写为含局部不连续映射的形式为：

$\varphi \left(X_0,\mu ,T\right)=\varphi \left(D\left(\varphi \left(X_0,\mu ,t_1\right),\mu \right),\mu ,T-t_1\right)$ (17)

其中局部不连续映射为：

$D\left(X,\mu \right)=\varphi \left(\varphi \left(R\left(\varphi \left(\varphi \left(X,\mu ,\delta \right),\mu ,\Delta t\right),\mu ,t\right),\mu ,-\Delta t\right),\mu ,-\delta \right)$ (18)

令 $\tilde{t}\left(X,\mu \right)=\delta +\Delta t$，(18)可写为：

$D\left(X,\mu \right)=\varphi \left(R\left(\varphi \left(X,\mu ,\tilde{t}\left(X,\mu \right)\right),\mu ,t\right),\mu ,-\tilde{t}\left(X,\mu \right)\right)$ (19)

将 $D\left(X,\mu \right)$ 关于X求偏导，并且由 $\frac{\partial \varphi }{\partial X}=I$， $\frac{\partial \varphi }{\partial \mu }=0$， $\frac{\partial \varphi }{\partial t}=F$，可以得到：

$D_X\left(X,\mu \right)=R_X\left(X_1,\mu \right)+\left(R_t\left(X_1,\mu \right)+R_X\left(X_1,\mu \right)F_1-F_2\right){\tilde{t}}_X\left(X,\mu \right)$ (20)

把点 ${\bar{X}}_1^{*}$ 带入(20)式，得到：

$D_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)=R_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)+\left(R_t\left({\bar{X}}_1^{*},{\mu }^{*}\right)+R_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)F_1^{*}-F_2^{*}\right){\tilde{t}}_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)$ (21)

式中 $F_1^{*}$， $F_2$ 分别为点 ${\bar{X}}_1^{*}$ 和点 ${\bar{X}}_2^{*}$ 处的向量场。

在(21)式中只有 ${\tilde{t}}_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)$ 的表达式是未知的，下面推导这个表达式。把t看作X和 $\mu$ 的函数，则 $t=\tilde{t}\left(X,\mu \right)$ 可以看作方程 $\tilde{h}\left(\varphi \left(X,\mu ,t\right),\mu ,t\right)=0$ 的解，由隐函数定理得到：

${\tilde{t}}_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)=-\frac{{\tilde{h}}_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)}{{\tilde{h}}_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)F_1^{*}+{\tilde{h}}_t\left({\bar{X}}_1^{*},{\mu }^{*}\right)}=-\frac{h_{\hat{X}}\left({\bar{X}}_1^{*},{\mu }^{*}\right)}{h_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)F_1^{*}-V\left(t_1\right)}$ (22)

将(22)式带入(21)式得到：

$D_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)=R_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)+\frac{\left(F_2^{*}-R_t\left({\bar{X}}_1^{*},{\mu }^{*}\right)-R_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)F_1^{*}\right)h_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)}{h_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)F_1^{*}-V\left(t_1\right)}$ (23)

则在不动点 $X={\bar{X}}_0$ 处的随机线性化矩阵为：

${\varphi }_X\left({\bar{X}}_0,\mu ,T\right)={\varphi }_X\left({\bar{X}}_2^{*},\mu ,T-t_1\right)D_X\left({\bar{X}}_1^{*},{\mu }^{*}\right){\varphi }_X\left({\bar{X}}_0,\mu ,t_1\right)$ (24)

根据文献 [18]，(24)式中的 ${\varphi }_X\left({\bar{X}}_0,\mu ,t_1\right)=\Phi \left(t_1\right)$ 和 ${\varphi }_X\left({\bar{X}}_2^{*},\mu ,T-t_1\right)=\Phi \left(T-t_1\right)$。

根据上面的推导给出如下的命题1。基于此命题来研究横截周期解的稳定性及进一步的分岔行为。

命题1：当轨道横截穿越碰撞面时，系统(3)和(4)相应的含参数随机Poincaré映射的线性化矩阵为：

${\varphi }_X\left({\bar{X}}_0,\mu ,T\right)=\Phi \left(T-t_1\right)D_X\left({\bar{X}}_1^{*},{\mu }^{*}\right)\Phi \left(t_1\right)$ (25)

5. 数值模拟

选取系统的一组参数 $m=2.2$， $K_{12}=K_{21}=-0.2$， $K_{22}=8$， $\xi =0.2$， $\omega =2.2$， $R=0.8$，以间隙d作为变化参数来观察确定性系统的稳定性。根据确定性系统不动点处的线性化矩阵 [18]，通过数值计算，当 $d=d^{*}=\text{1}\text{.0738483}$，线性化矩阵的特征值分别为：

${\lambda }_1\left(d^{*}\right)=-\text{1}\text{.00000084}$， ${\lambda }_2\left(d^{*}\right)=-\text{0}\text{.20472910}$， ${\lambda }_{3,4}\left(d^{*}\right)=\text{0}\text{.37149808}\pm \text{0}\text{.42450952i}$

由此可见系统发生了倍化分岔。

确定性系统随间隙d的变化的分岔图如图3(a)，从图中可以看出，当间隙参数d在区间 $\left(1.05,1.0738\right)$ 上变化时，质块 $m_1$ 具有稳定的周期1-1运动；当d逐渐增大至分岔点 $d=d^{*}$ 时，质块 $m_1$ 发生倍周期分岔，进入稳定的周期2-2运动；随着d的进一步增大，质块 $m_1$ 继续发生倍周期分岔，进入周期4-4运动，最终倍化通向混沌运动。

(a) (b) (c) (d)

Figure 3. (a) Bifurcation diagram of deterministic system; (b) Bifurcation diagram when $\sigma =1$ ; (c) Bifurcation diagram when $\sigma =2.5$ ; (d) Bifurcation diagram when $\sigma =3.5$

图3. (a) 确定性系统的分岔图；(b) $\sigma =1$ 时系统分岔图；(c) $\sigma =2.5$ 时系统分岔图；(d) $\sigma =3.5$ 时系统分岔图

首先在低强度随机约束下作出系统以间隙d为分岔参数的分岔图，此时取随机扰动尺度 $\sigma =1$，如图3(b)，此时随机干扰对系统的分岔的影响不显著。增大随机扰动尺度，如图3(c)和图3(d)，此时可以看到分岔图发生了明显的变化，随机干扰导致系统的分岔图由确定的单线逐渐扩散为带状，同时，通过观察放大图，确定系统中孤立的分岔点演化成在原分岔点附近的分岔区间。

通过图3(c)和图3(d)并不能清晰地观察到随机约束对系统稳定性的影响，只能模糊地观察到分岔区间，下面，利用命题1中的随机线性化矩阵(25)来观察其特征值，见表1。

Table 1. Eigenvalue data corresponding to “on the unit circle” when d = 1.0738483

表1. $d=1.0738483$ 时对应“单位圆上”的特征值数据

在线性化矩阵(25)中，有随机速度和随机加速度，在模拟随机速度和随机加速度中分别以时间 $t=T/2$， $t=3T/2$， $t=5T/2$ 取值，求出间隙为 $d=1.0738483$ 时不同随机扰动尺度下的特征值，见表1，列出了系统落在“单位圆上”的特征值。观察表中的特征值发现，在此时的间隙参数下，与确定性系统相比，随着随机扰动尺度的增加，特征值基本落在单位圆外，极偶然会落在单位圆内。下面以时间取固定 $t=3T/2$ 时，求出不同随机扰动尺度下取不同间隙参数的特征值，见表2，进一步观察特征值的变化。观察发现，在 $d=1.070$ 和 $d=1.072$，随机扰动尺度为 $\sigma =1$ 时特征值落在单位圆内，随机扰动尺度逐渐增大，特征值落在单位圆外并逐渐远离单位圆，在 $d=1.074$， $d=1.076$ 和 $d=1.078$ 时，特征值均落在单位圆外，并随着随机扰动尺度增大，特征值距离单位圆越来越远，显然，系统的分岔行为提前发生了。

Table 2. Eigenvalue data corresponding to “on the unit circle” under different clearance parameters

表2. 不同间隙参数下对应“单位圆上”的特征值数据

6. 结论

本文推导了一类2阶可微的均值回复随机过程，求出了所服从的对应的分布，在此基础上研究了一类随机边界约束的二自由度悬臂梁系统的稳定性。推导得到了横截周期解含参数情形下的随机线性化矩阵。基于线性化矩阵确定了周期解发生倍化分岔的分岔点。通过对比确定性系统的分岔图以及不同扰动的随机系统的分岔图，发现随机分岔图不再是清晰的线型，随着随机扰动尺度的增大，系统的分岔行为提前发生，分岔线逐渐成为带状，分岔点也逐渐演化为分岔区间。本文的研究方法可为模拟间隙中的不确定性提供一定的思路。

基金项目

本文的研究得到国家自然科学基金(11872264)的资助，在此作者表示衷心的感谢。