高斯色噪声激励下的高速转子轴承系统的P-分岔分析

P-Bifurcation Analysis of High-Speed Rotor-Bearing System under Excited byGaussian Color Noise

摘要: 对高速转子轴承系统引入高斯色噪声，根据统一色噪声近似原理将色噪声白化，得到白噪声激励的等效非线性模型，运用拟不可积Hamilton理论随机平均法得到Itô微分方程，并求出FPK方程和对应的概率密度函数。最后根据概率密度函数分析随机P-分岔，并进行数值模拟验证。

Abstract: The Gaussian color noise is introduced into the High-Speed Rotor-Bearing system, and the equivalent nonlinear model of white noise excitation is obtained according to the uniform colored noise approximation theory. The Itô differential equation is obtained by the stochastic average method of quasi-non-integrable Hamilton theory, and the FPK equation and the corresponding probability density function are obtained. Finally, the stochastic P-bifurcation is analyzed according to the probability density function and verified by numerical simulation.

1. 引言

高速转子–轴承系统作为现代工业机械中的核心部件，广泛应用于航空航天、能源发电和精密制造等领域。其运行的安全性和可靠性直接关系到整个设备的性能和寿命。然而，在实际运行环境中，此类系统不可避免地会受到各种随机因素的影响，例如：材料质地的不均匀、温度和气压的波动、流体湍流以及轴承润滑膜的随机扰动等。这些随机扰动并非简单的背景噪声，而是能够深刻改变系统动力学行为的激励源，可能诱发复杂的非线性现象和随机失稳，从而对设备造成潜在的破坏。因此，深入研究高速转子–轴承系统在随机激励下的动力学特性，尤其是其稳定性与随机分岔行为，具有重要的理论意义和工程价值。

近年来，不少学者对转子系统进行了广泛深入的研究。文献[1]考虑油膜力，分析了双盘转子系统的动力学，基于分岔图、Poincare图、最大Lyapunov指数研究裂纹深度定子刚度对转子系统的非线性影响；文献[2]考虑了受外力激励的转子系统，结合非线性动力学提出了用于阻尼优化的观点，并通过数值模拟验证了结论的有效性，文献[3]对具有一定几何结构和力学的高速转子进行分析，结合理论和数值仿真研究了转子系统的动力学特性。文献[4]对转子轴承系统引入高斯白噪声，根据随机动力学理论建立随机激励耗散系统，对系统进行稳定性和分岔分析，并进行仿真验证。

现实环境中，任何机械转子系统都不可避免受随机因素影响，如风速、温度等不可测因素，文献[5]根据非线性动力学理论分析了正弦扰动因素作用下的永磁同步电动机的非线性动力学特征，文献[6]讨论了噪声干激励下的PMSM模型的随机分岔，基于包映射进行数值模拟了演化规律。基于上述分析，考虑随机激励因素具有现实意义。实际上，随机动力系统应用于机械、生物、金融等诸多领域，文献[7]对调速器系统引入高斯色噪声，由随机平均法和极坐标变换分析了系统稳定性和分岔，并进行Monte Carlo模拟验证，文献[8]考虑色噪声激励的水轮机调节系统，得到系统随机分岔条件，并模拟不同噪声下的分岔情况；文献[9]对一类SIR传染病模型引入白噪声，分析了系统稳定性的条件和随机P-分岔，并由数值模拟验证，文献[10]分析了白噪声激励的网络传染病模型得到随机灭绝和持久的条件；文献[11]基于经济周期动力模型引入色噪声，应用Hamilton理论体系由Lyapunov指数分析了系统稳定性条件。在物理应用中，噪声激励耗散的Hamilton系统可表示非线性随机动力学[12]，通常Legendre方程可变换得到Hamilton方程．文献[13]介绍了几种非高斯过程的Hamilton方法并给出如Van der Pol振子系统、Rayleigh振子系统的算例，得到概率密度函数，文献[14]讨论了谐和项与宽带噪声联合激励下的拟可积Hamilton系统最优时滞控制问题，得到Hamilton函数。受到上述文献的启发，本文对高速转子轴承系统引入色噪声，建立随机激励下耗散的Hamilton系统并讨论不同噪声强度下的系统分岔情况。

已有的研究工作主要通过P-分岔分析来揭示系统在噪声激励下的行为转变，这为理解系统的长期统计行为提供了宝贵的见解。然而，单纯的P-分岔分析未能全面揭示系统平凡解的根本稳定性。例如，在某些情况下，即使拓扑结构发生变化，系统也可能在更严格的意义上保持稳定。因此，将P-分岔与D-分岔分析相结合，可以提供对随机分岔现象更为全面和深刻的理解，实现对系统稳定性判定的相互印证。这是本报告的核心出发点和主要研究缺口。

2. 随机模型的建立

考虑转子系统材料不均及易受外界因素如温度、气压等影响，把影响因素视作高斯色噪声，引入随机项建立如下无阻尼的运动方程：

$M\ddot{q}\left(t\right)+C\dot{q}\left(t\right)+Kq=F_c\eta \left(t\right)$ (1)

方程(1)中：$\eta \left(t\right)$ 为高斯色噪声，且$\langle \eta \left(t\right)\rangle =0$ ，$\langle \eta \left(t_1\right)\eta \left(t_2\right)\rangle =\frac{D}{\tau }{\text{e}}^{-\frac{\left|t_1-t_2\right|}{\tau }}$ ；$2D$ 为$\eta \left(t\right)$ 的噪声强度；$q$ 表示广义坐标，$M$ 表示广义质量，$C$ 表示广义阻尼，$K$ 表示广义刚度，$F_c$ 表示广义强迫力。

令$m=\frac{C}{M}$ ，$n=\frac{K}{M}$ ，$k=\frac{f}{M}$ ，由色噪声近似原理[15]，系统模型(1)可近似表示为以下方程：

$\{\begin{array}{l}\ddot{q}\left(t\right)+m\dot{q}\left(t\right)+nq=k\eta \left(t\right)\\ \dot{\eta }\left(t\right)=-\eta \left(t\right)/\tau +\xi \left(t\right)/\tau \end{array}$ (2)

其中：$\xi \left(t\right)$ 为高斯色噪声，且$\langle \xi \left(t\right)\rangle =0$ ，$\langle \xi \left(t_1\right)\xi \left(t_2\right)\rangle =2D\delta \left(t_1-t_2\right)$ 。

为了方便计算，把方程(2)表示为一阶微分形式，令$p=\dot{q}$ ，得

$\{\begin{array}{l}\dot{q}=p\\ \dot{p}\left(t\right)=-mp-nq+k\eta \left(t\right)\\ \dot{\eta }\left(t\right)=-\eta \left(t\right)/\tau +\xi \left(t\right)/\tau \end{array}$ (3)

方程(3)中2式求$t$ 的微商，再利用3式消去变量$\eta$ ，可得

$\ddot{p}+\left[{\tau }^{-1}-c^{\prime }\left(p\right)\right]\dot{p}-{\tau }^{-1}c\left(p\right)=k{\tau }^{-1}\xi \left(t\right)$ (4)

其中：$c\left(p\right)=-mp-nq$ 。

由于$p$ 表示广义运动坐标，参数m、n可设定，则$c\left(p\right)$ 为有界项，即当$\tau \to 0$ 时$\xi \left(p,\tau \right)\to \infty$ ，所以模型(4)是一个过阻尼过程，可省去惯性项$\ddot{p}$ 。

令$\alpha \left(p,\tau \right)={\tau }^{-1}-c^{\prime }\left(p\right)$ ，则方程(4)可化为：

$\dot{p}=\frac{{\tau }^{-1}c\left(p\right)}{\alpha \left(p,\tau \right)}+\frac{k{\tau }^{-1}}{\alpha \left(p,\tau \right)}\xi \left(t\right)$ (5)

将方程(5)代入原模型(1)可得噪声白化的高速转子轴承系统模型：

$\{\begin{array}{l}\dot{q}=p\\ \dot{p}=-\frac{{\tau }^{-1}m}{\alpha \left(p,\tau \right)}p-\frac{{\tau }^{-1}n}{\alpha \left(p,\tau \right)}q+\frac{k{\tau }^{-1}}{\alpha \left(p,\tau \right)}\xi \left(t\right)\end{array}$ (6)

其中$\alpha \left(p,\tau \right)$ 为关于$p$ 的非线性函数，在邻域$U\left(0\right)$ 内，把方程(6)中的$\frac{{\tau }^{-1}}{\alpha \left(p,\tau \right)}$ 进行泰勒展开，等效近似为关于$p$ 的线性形式：

$\frac{{\tau }^{-1}}{\alpha \left(p,\tau \right)}=\frac{{\tau }^{-1}}{{\tau }^{-1}+m}-\frac{{\tau }^{-1}n}{{\left({\tau }^{-1}+m\right)}^2}p+\cdots \approx k_1-k_{2}p$ (7)

记$k_1=\frac{{\tau }^{-1}}{{\tau }^{-1}+m}$ 、$k_2=\frac{{\tau }^{-1}n}{{\left({\tau }^{-1}+m\right)}^2}$ ，把方程(7)代入方程(6)中，得到等效非线性近似模型：

$\{\begin{array}{l}\dot{q}=p\\ \dot{p}=-k_{1}nq+\left[k_{2}nq+k_{2}mp-k_{1}m\right]p+k\left(k_1-k_{2}p\right)\xi \left(t\right)\end{array}$ (8)

引入哈密顿函数$\dot{q}=\frac{\partial \overline{H}}{\partial p}$ $\dot{p}=-\frac{\partial \overline{H}}{\partial q}$ ，上述模型可转化为单自由度的哈密顿系统模型：

$\{\begin{array}{l}\dot{q}=\frac{\partial \overline{H}}{\partial p}\\ \dot{p}=-k_{1}n\frac{\partial \overline{H}}{\partial q}+\left[k_{2}nq+k_{2}mp-k_{1}m\right]\frac{\partial \overline{H}}{\partial p}+k\left(k_1-k_{2}p\right)\xi \left(t\right)\end{array}$ (9)

其中：$q$ 为广义位移；$p$ 为广义速度；$\overline{H}\left(q,p\right)$ 为系统总能量，将模型(9)化为Stratonovich随机微分方程：

$\{\begin{array}{l}\text{d}q=\frac{\partial \overline{H}}{\partial p}\text{d}t\\ \text{d}p=-\left[k_{1}n\frac{\partial \overline{H}}{\partial q}-\left(k_{2}nq+k_{2}mp-k_{1}m\right)\frac{\partial \overline{H}}{\partial p}\right]\text{d}t\\ +\frac{1}{2}\sqrt{2D}k\left(k_1-k_{2}p\right)\frac{\partial \sqrt{2D}k\left(k_1-k_{2}p\right)}{\partial p}+\sqrt{2D}k\left(k_1-k_{2}p\right)\circ \text{d}B\left(t\right)\end{array}$ (10)

模型(10)中的Wong-Zakai修正项为：

$\frac{1}{2}\sqrt{2D}k\left(k_1-k_{2}p\right)\frac{\partial \sqrt{2D}k\left(k_1-k_{2}p\right)}{\partial p}=-D{k}^2{k}_1{k}_2+D{k}^2{k}_2^{2}p$

式中$-D{k}^2{k}_1{k}_2$ 为保守成分的修正项；$D{k}^2{k}_2^{2}p$ 为消散成分的修正项，则系统(10)进一步表示为：

$\{\begin{array}{l}\text{d}q=\frac{\partial H}{\partial p}\text{d}t\\ \text{d}p=-\left[\frac{\partial H}{\partial q}+m\left(p,q\right)\frac{\partial H}{\partial p}\right]\text{d}t+\sqrt{2D}k\left(k_1-k_{2}p\right)\text{d}B\left(t\right)\end{array}$ (11)

其中$\frac{\partial H}{\partial q}=k_{1}n\frac{\partial \overline{H}}{\partial q}+D{k}^2{k}_1{k}_2$ ；$\frac{\partial H}{\partial p}=\frac{\partial \overline{H}}{\partial p}$ ；$m\left(p,q\right)=-\left[k_{2}nq+k_{2}mp-k_{1}m+D{k}^2{k}_2^2\right]$ 。

此时系统的哈密顿函数为：

$H=\frac{1}{2}{p}^2+\frac{1}{2}{k}_{1}n{q}^2+D{k}^2{k}_1{k}_{2}q$ (12)

根据拟不可积Hamilton理论，方程弱收敛于一维的Itô微分方程，且系统相关的Itô微分方程可表示为：

$\text{d}H=m\left(H\right)\text{d}t+\sigma \left(H\right)\text{d}B\left(t\right)$

记$u=\frac{1}{2}{k}_{1}n$ ，$v=D{k}^2{k}_1{k}_2$ ，则公式(12)进行恒等变换得

$2H+\frac{v^2}{2u}=p^2+{\left({\left(2u\right)}^{1/2}q+\frac{v}{{\left(2u\right)}^{1/2}}\right)}^2$

令$p={\left(2H+\frac{v^2}{2u}\right)}^{1/2}\cos \theta$ ，$q=\left[{\left(2H+\frac{v^2}{2u}\right)}^{1/2}\sin \theta -\frac{v}{{\left(2u\right)}^{1/2}}\right]/{\left(2u\right)}^{1/2}$ ，根据拟哈密顿理论随机平均法得：

$\begin{array}{c}m\left(H\right)=\frac{1}{T\left(H\right)}{\displaystyle {\int }_{\Omega }\left[-m\left(p,q\right)+D{k}^2{\left(k_1-k_{2}p\right)}^2\right]/p\text{d}q}\\ =\frac{k_{2}nv}{8\sqrt{2}{u}^2}\left(4H^{2}u+v^{2}H\right)-\frac{k_{1}m}{4u}\left(4H^{2}u+v^{2}H\right)\\ +\frac{Dk{k}_2^2}{2u}\left(4H^{2}u+v^{2}H\right)+Dk{k}_1^2\\ =m_1{H}^2+m_{2}H+m_3\end{array}$

$\begin{array}{c}{\sigma }^2\left(H\right)=\frac{1}{T\left(H\right)}{\displaystyle {\int }_{\Omega }\left(2D{k}^2{\left(k_1-k_{2}p\right)}^2\right)/p\text{d}q}\\ =\frac{Dk{k}_2^2}{2u}\left(4H^{2}u+v^{2}H\right)+2Dk{k}_1^2\\ =n_1{H}^2+n_{2}H+n_3\end{array}$

其中：$T\left(H\right)={\displaystyle {\int }_{\Omega }\frac{1}{\partial H/\partial p}\text{d}q}=\frac{1}{\sqrt{2u}}\pi$ ，$\Omega =\left\{q|H\left(q,0\right)\le H\right\}$ ；

$m_1=\frac{k_{2}nv}{2\sqrt{2}u}+2Dk{k}_2^2-k_{1}m$ ；$m_2=\frac{k_{2}n{v}^3}{8\sqrt{2}{u}^2}+\frac{Dk{k}_2^2{v}^2}{2u}-\frac{k_{1}m{v}^2}{4u}$ ；$m_3=Dk{k}_1^2$

$n_1=\text{2}Dk{k}_2^2$ ；$n_2=\frac{Dk{k}_2^2{v}^2}{2u}$ ；$n_3=2Dk{k}_1^2$ 。

3. P-分岔分析

已知系统的平均Itô微分方程可表示为：

$\text{d}H=m\left(H\right)\text{d}t+\sigma \left(H\right)\text{d}B\left(t\right)$ (13)

概率密度法是研究系统P-分岔的有效方法，思路为首先根据Itô微分方程得到FPK方程进而求出平稳概率密度函数，由(13)得FPK方程为：

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

方程(14)的解为：

$p\left(H\right)=C\exp \left(-\lambda \left(H\right)\right)$

其中：

$\begin{array}{c}\lambda \left(H\right)={\displaystyle {\int }_0^H\left\{\left[\frac{\text{d}{\sigma }^2\left(H\right)}{\text{d}H}-2m\left(H\right)\right]/{\sigma }^2\left(H\right)\right\}\text{d}H}\\ ={\displaystyle {\int }_0^H\left[2n_{1}H+n_2-2m_1{H}^2-2m_{2}H-2m_3\right]/\left(n_1{H}^2+n_{2}H+n_3\right)\text{d}H}\\ ={\displaystyle {\int }_0^H-\frac{2m_1}{n_1}+\frac{\left(\frac{2m_1-2m_2}{n_1}+\frac{2m_1{m}_2}{n_1^2}\right)H+\frac{n_2-2m_3}{n_1}+\frac{2m_1{m}_2}{n_1{n}_3}}{H^2+\frac{n_2}{n_1}H+\frac{n_3}{n_1}}\text{d}H}\\ =-\frac{2m_1}{n_1}H+{\displaystyle {\int }_0^H\frac{MH+N}{H^2+UH+V}\text{d}H}\\ =-\frac{2m_1}{n_1}H+{\displaystyle {\int }_0^H\frac{u_1}{H+u_2}+\frac{v_1}{H+v_2}\text{d}H}\\ =-\frac{2m_1}{n_1}H+u_1\ln \left|H+u_2\right|+v_1\ln \left|H+v_2\right|-u_1\ln \left|u_2\right|-v_1\ln \left|v_2\right|\end{array}$

记：

$\begin{array}{l}M=\frac{2m_1-2m_2}{n_1}+\frac{2m_1{m}_2}{n_1^2},N=\frac{n_2-2m_3}{n_1}+\frac{2m_1{m}_2}{n_1{n}_3},U=\frac{n_2}{n_1},V=\frac{n_3}{n_1}\\ u_1=\frac{u_{2}M-N}{\sqrt{U^2-4V}},v_1=\frac{-v_{2}M+N}{\sqrt{U^2-4V}},u_2=\frac{U-\sqrt{U^2-4V}}{2},v_2=\frac{V-\sqrt{U^2-4V}}{2}\end{array}$

进而，

$\begin{array}{c}p\left(H\right)=C\exp \left(-\lambda \left(H\right)\right)\\ =C{\left|u_2\right|}^{u_1}\cdot {\left|v_2\right|}^{v_1}\cdot {\left|H+u_2\right|}^{-u_1}\cdot {\left|H+v_2\right|}^{-v_1}\exp \left(\frac{2m_1}{n_1}H\right)\end{array}$

其中$C$ 为归一化常数，此时$p,q$ 的联合概率密度为：

$\begin{array}{c}p\left(p,q\right)=C{\left|u_2\right|}^{u_1}\cdot {\left|v_2\right|}^{v_1}\cdot {\left|H+u_2\right|}^{-u_1}\cdot {\left|H+v_2\right|}^{-v_1}\exp \left(\frac{2m_1}{n_1}H\right)\\ =C{\left|u_2\right|}^{u_1}\cdot {\left|v_2\right|}^{v_1}\cdot {\left|\frac{1}{2}{p}^2+u{q}^2+vq+u_2\right|}^{-u_1}\cdot {\left|\frac{1}{2}{p}^2+u{q}^2+vq+v_2\right|}^{-v_1}\\ \cdot \exp \left[\frac{2m_1}{n_1}\left(\frac{1}{2}{p}^2+u{q}^2+vq\right)\right]\end{array}$

取一组参数$m=4$ $n=2$ $k=8$ ，当相关时间为$\tau =0.1$ 时，噪声强度$D$ 分别取0.5，0.7，0.9，1.0时的联合概率密度模拟图像，如下。

图1(a)为$D=0.5$ 时联合概率密度函数，图像呈单峰状，且最高点为有限值，系统较大可能处于稳定状态；如图1(b)所示，当噪声强度增加到$D=0.7$ 时，单峰最高点下降，周围出现凸起变成火山口形状，概率密度函数出现新的极值，即系统发生了P-分岔；如图1(c)，随噪声强度的继续增加，单峰周围凸起的幅度也不断增大，火山口形状更加显现，此时系统可能变得不稳定；图1(d)为$D=1.0$ 时联合概率密度函数，如图所示，除了单峰周围形成较为明显的火山口外，单峰中心也开始凹陷形成一个新的火山口，概率密度函数极值再次增加，即系统又一次发生了P-分岔。联合概率密度函数的变化反映系统性态的改变，说明一定的噪声强度对系统稳定性具有破坏作用。

(a) $D=0.5$ (b) $D=0.7$

(c) $D=0.9$ (d) $D=1.0$

Figure 1. Joint probability density function under different noise intensities

图1. 不同噪声强度下的联合概率密度函数

Figure 2. Stationary probability density functions under different noise intensities

图2. 不同噪声强度下的平稳概率密度函数

如图所示，图2表示不同噪声强度下的平稳概率密度函数，当噪声强度$D=0.6$ 时概率密度函数有一个单峰极大值，当噪声强度$D=0.7$ 时，概率密度函数极大值减少，单峰坡度变缓并出现新的极值点，极值点的增加说明系统发生P-分岔，也说明系统性态发生改变，此时系统可能变得不稳定，当噪声强度增加到0.8，0.9时，这种变化趋势更加显现，系统失稳更加明显，这也与上述结论相互对应。

4. 结论

转子系统材料不均和外界温度、气压等因素可能会使系统不稳定发生随机分岔，本文把影响因素视作高斯色噪声，讨论了色噪声激励下高速转子轴承系统的随机P-分岔分析。首先根据统一色噪声原理将噪声白化，运用拟不可积哈密顿理论随机平均法得到Itô微分方程。然后根据概率密度函数分析了P-分岔，并进行数值模拟，结果表明一定噪声强度会对系统性态产生影响，且随噪声强度的增大影响更为显现，说明较大的噪声强度使系统变得不稳定。

基金项目

广东理工学院科技项目(2024QNZK004)，广东理工学院高等教育教学改革项目(JXGG2024151)，广东理工学院课堂教学改革项目(KTJXGG202579)。

NOTES

^*通讯作者。