引言

众所周知，蚯蚓等软体动物通过蠕动实现运动，这种运动形式看似简单，却凭借其一大优势-无足，吸引了众多学者的关注^[1-6]. 受此启发，一类新型的运动模型-振动驱动系统，近年来得到了广泛的研究^[7-25]. 这种运动模型继承了蚯蚓运动的优势-不需要外部驱动部件，如腿或轮子，因此易于实现结构的微型化，密封化. 可以预见，基于振动驱动这一思想的机器人能更好地适应多种环境，甚至有学者期望将其制成微型胶囊机器人进入人体肠道甚至血管中治疗疾病.

2005年，Chernousko^[7]研究了两个质量块系统的直线运动，通过选择合适的物理参数并合理地控制两个质量块之间的相互作用力，实现了该系统的最快平均速度.在此基础上，其首次提出了由一个内部质量块和一个外部箱型刚体构成的振动驱动系统模型^[8].在该模型中，内部质量块相对于外部箱体做周期运动，通过合理设计内部质量块的周期运动形式，外部箱体可以在内部质量块对其作用力以及支撑面对其摩擦力的共同作用下，实现周期平移运动.因而，对质量块的相对运动形式进行设计和优化以实现系统最大平均速度是该类系统研究的一个关键问题.纵观之前的研究，质量块的控制模式大致分为两类：对称控制和非对称控制.对称控制意味着内部质量块的相对运动形式具有对称性，也更容易实现，比如正弦驱动，但是往往需要环境提供能使系统发生刚体位移的异性摩擦.2012年，Fang等^[9]在研究三模块振动驱动系统时使用了正弦驱动，他们发现通过优化各模块中正弦驱动的初始相位差，能够有效地提高系统的稳态平均速度.2014年，Fang等^[10]研究了正弦驱动下，单模块振动驱动系统的运动，重点从滑移分岔的角度分析了系统的粘滑运动，加深了人们对该类系统粘滑运动的认识. 对称驱动的研究总的来说并不如非对称驱动的研究多.早在2005年Chernousko^[7]就提出了两种非对称驱动模式-两相驱动和三相驱动，对实现振动驱动系统的运动都行之有效.2006年，Chernousko^[11]研究了两相驱动模式，这种驱动模式下，内部质量块一个周期内的相对运动被分为两段，每段上的相对运动速度均为常数. 以实现系统的最大稳态平均速度为目标，对两相驱动参数进行了优化.2007年，Chernousko^[12]研究了三相驱动下振动驱动系统的运动.该种情形下，内部质量块一个周期内的相对运动被分为三段，每段上的加速度大小为常数.由于库伦干摩擦的存在，以及不同时间域内，内部质量块的相对加速度不同，加上系统速度的不断变化，因而系统会发生复杂的粘滑运动.2011年，Fang等^[13]运用平均法研究了三相驱动的单模块振动驱动系统在不存在粘滑效应时的速度.但是当考虑粘滑效应时，平均法不再适用，根据粘滑运动发生情况，系统运动可以被细分成8种情况.2011年，Fang等^[14]在研究两模块振动驱动系统时，在两个模块中也分别施加了这种三相驱动，通过调节两组振动的相位差实现了系统的较快运动速度. 可见，三相驱动在振动驱动系统的控制中占有重要的地位.此外，2006年Li等^[15]使用的四步驱动以及 2011年 Huda等^[16] 使用的四相驱动也属于非对称驱动.

除了上文中内部质量块的相对运动形式以外，另外一个对振动驱动系统的运动行为起很大影响的是系统受到的外部摩擦.2007年，Chernousko^[12]研究了非对称黏性摩擦下振动驱动系统的运动，这种环境模拟了系统在液体环境中的运动.研究针对振动驱动系统速度进行了控制参数优化，并加以实验验证. 2009年，Bolotnik等^[17]研究了平方阻尼作用下系统的运动情况，这种运动常常在系统运动速度较快、雷诺数较大的环境下发生.研究表明，在平方阻尼作用下，即使是各向同性的阻力，系统也能实现向前运动.2011年，Fang等^[13]研究了振动驱动系统在非对称库伦干摩擦下的运动情况，由于干摩擦存在，观察到了前面所述的粘滑运动.

事实上，粘滑效应至今为止仍然是振动驱动系统中的一大难题.一方面，粘滑效应使得运动情况复杂多样，难以分析，另一方面，充分利用粘滑效应也是提高系统运动性能的重要措施.粘滑运动往往都是在库伦摩擦之下发生，这一点，不仅在振动驱动系统中得到体现，而且在其他系统中也有发现并被研究^{[18, 26-28]}.库伦摩擦之所以能引起粘滑运动，是因为其不连续性.从运动状态场的角度看，这种由摩擦力引起的不连续系统，是一类典型的Filippov系统^[29-30]，而Filippov系统最突出的就是利用滑移和滑移分岔理论巧妙地分析系统中的不连续性对运动的影响. 2010年，Guardia等^[26]通过滑移分岔理论研究了一个干摩擦下的弹簧振子，研究证实系统存在两参数滑移分岔，并得到了其中一条余维一分岔曲线的解析表达式.2007年，Kowalczyk等^[27]也通过Filippov系统模型研究了干摩擦下的弹簧振子，发现了系统中的余维二分岔，并解释了可能发生的混沌现象. Filippov系统理论为研究干摩擦下的弹簧振子这类系统行为，提供了重要的分析方法^{[18, 26-28]}.受此启发，2011年，Fang等^[10]也利用Filippov系统描述了非对称干摩擦下受正弦驱动的振动驱动系统，首次从滑移分岔角度阐明了粘滑效应下系统的运动情况，优化结果显示，合理优化系统驱动参数和摩擦系数，能实现系统无后退移动.

回顾振动驱动系统十年来的研究^[7-25]，无论是内部驱动的设计，还是外部阻力的分类，抑或是研究方法的选用，都得到了长足的发展. 值得注意的是，绝大部分研究更多地是关注了摩擦力为各向异性时的情况，而对于同性摩擦的情形讨论较少.因此，本文重点研究了各向同性干摩擦下的振动驱动系统的运动.由于摩擦的各向同性，本文中将采用三相驱动这一非对称驱动以保证系统定向移动的实现.此外，干摩擦条件下，系统的不连续性会诱发粘滑运动，这种运动的复杂性也对研究工作形成了挑战.借助分岔理论可知，本文中所考虑的振动驱动系统也是Filippov系统，因此可以从滑移分岔的角度对系统的运动进行分类.这种分类有效地揭示了各向同性干摩擦下振动驱动系统多样的运动行为，对参数的设计起到一定的指导作用.最后，我们用数值方法对理论分析得到的运动分类以及运动特性进行了验证.

1 振动驱动系统 1.1 动力学方程

振动驱动系统的基本模型如图 1所示. 系统由两个部分组成，分别为质量为M的刚性箱体以及质量为m的内部质量块. 由前人的工作可知，通过m相对于M的特定振动，可以实现系统在水平直线上的定向移动，由此称其为振动驱动系统. 易知，系统有两个自由度，因而可以取两组独立的广义坐标，分别为m相对于M的位移ξ，以及M的绝对位移y，来表示系统的运动情况.

由牛顿第二运动定律，可得运动方程

$(M+m)\ddot{y}=-m\ddot{\xi }-R(\dot{y})$

(1)

式中，$R(\dot {y})$表示系统所受的库伦摩擦力，从而有

$R(\dot{y})=\left\{ \begin{array}{*{35}{l}} f(M+m)g\text{sgn}(\dot{y}), & \dot{y}e0 \\ R(0), & \dot{y}=0 \\ \end{array} \right.$

式中sgn (x)表示符号函数. 式(1)两边同时除以(M+m)，得到

$\ddot{y}=-\mu \ddot{\xi }-r(\dot{y})$

(2)

式中

$\mu =\frac{m}{M+m},r(\dot{y})=\left\{ \begin{array}{*{35}{l}} fg\text{sgn}(\dot{y}), & \dot{y}e0 \\ r(0), & \dot{y}=0 \\ \end{array} \right.$

需要注意的是，上面的R (0)表示速度为0时，系统所受的摩擦阻力. 因此，如果内部振动提供的驱动力不足以克服库伦摩擦，那么摩擦力就与该力等大反向，维持系统静止状态. 如果内部振动提供的驱动力大于最大静摩擦力，那么系统有运动的趋势，此时摩擦力的大小就是最大静摩擦力. 为了简便分析，本文中最大静摩擦力近似认为等于滑动摩擦力. 所以有

$r(0)=\frac{R(0)}{M+m}=\left\{ \begin{array}{*{35}{l}} -\mu \ddot{\xi }, & |\mu \ddot{\xi }|\le fg \\ -fg\text{sgn}\left( {\ddot{\xi }} \right), & |\mu \ddot{\xi }|>fg \\ \end{array} \right.$

(3)

式(2)和式(3)是振动驱动的基本力学模型.

1.2 三相驱动

假设箱体内允许质量块运动的区间长度为L，那么对于任意时刻t，有

$0\le \xi (t)\le L$

(4)

并且由振动的周期性，有

$\left. \begin{array}{*{35}{l}} \xi (0)=\xi (T)=0 \\ \dot{\xi }(0)=\dot{\xi }(T)=0 \\ \end{array} \right\}$

(5)

为了充分利用箱体的空间，同时又保证箱体和质量块之间不会发生碰撞，我们可以知道，存在某个时刻$\theta \in (0,T)$ ，使得

$\xi (\theta )=L,\dot{\xi }(\theta )=0$

(6)

在满足式(4) ∼式(6)的条件下，Chernousko^[8]提出了三相驱动模式.

三相驱动模式的基本思想是将内部质量块一个周期的相对振动在时间上分为三段，三段区间的长度分别为τ₁,τ₂ ,τ₃ . 并且，在每段时间区间内，质量块m相对于M均做匀变速直线运动，由式(5)∼式(7)易知，第一、三段中m相对运动的加速度为正常数，第二段中m相对运动的加速度为负常数.一个周期内，m相对运动的数学表达式如下

$\ddot{\xi }(t)=\left\{ \begin{matrix} {{w}_{1}}, & t\in (0,{{\tau }_{1}}) \\ -{{w}_{2}}, & t\in ({{\tau }_{1}},{{\tau }_{1}}+{{\tau }_{2}}) \\ {{w}_{3}}, & t\in ({{\tau }_{1}}+{{\tau }_{2}},T) \\ \end{matrix} \right.$

(7)

其中，$w_i > 0 (i = 1,2,3)$，$T = \tau _1 + \tau _2 + \tau _3$. 为了便于后面的分析，将式(7)改写，并将$\ddot {\xi }(t)$的表达式拓宽到整个时间域上，有

$\ddot{\xi }(t)=\left\{ \begin{align} & {{w}_{1}}-h(t-{{t}_{1}})({{w}_{2}}+{{w}_{1}})+ \\ & h(t-{{t}_{2}})({{w}_{2}}+{{w}_{3}}), t\in [0,T) \\ & \ddot{\xi }(t(\text{mod}T)), t\notin [0,T) \\ \end{align} \right.$

(8)

其中h(t)表示Heaviside函数，mod 表示取模运算. 这里，为了方便$t_1 = \tau _1$，$t_2 = \tau _1 + \tau _2 $.

三相驱动模式下内部质量块相对运动图像如图 2所示.

三相驱动中，独立的控制变量有3个，本文使用$w_i (i = 1,2,3)$作为控制的自变量. 当一组$w_i (i = 1,2,3)$选定以后，其他的变量可以根据式(4)∼式(6)导出，分别为

$\left. \begin{align} & \theta =\sqrt{\frac{2L({{w}_{1}}+{{w}_{2}})}{{{w}_{1}}{{w}_{2}}}} \\ & {{\tau }_{1}}=\sqrt{\frac{2L{{w}_{2}}}{{{w}_{1}}({{w}_{1}}+{{w}_{2}})}} \\ & {{\tau }_{2}}=\sqrt{\frac{2L}{{{w}_{2}}}}(\sqrt{\frac{{{w}_{3}}}{{{w}_{2}}+{{w}_{3}}}}+\sqrt{\frac{{{w}_{1}}}{{{w}_{1}}+{{w}_{2}}}}) \\ & {{\tau }_{3}}=\sqrt{\frac{2L{{w}_{2}}}{{{w}_{3}}({{w}_{2}}+{{w}_{3}})}} \\ \end{align} \right\}$

(9)

1.3 Filippov系统

振动驱动系统的运动情况不仅取决于内部质量块的相对振动，而且与其所受的外部阻力有关. 本文中系统所受的阻力为库伦干摩擦，干摩擦的复杂性在于其不连续性，即当系统速度方向改变时，摩擦力的变化是阶跃的. 因此，就其本质而言，本文所研究的系统，也就是式(2)和式(8) 所决定的振动驱动系统，是一种右端不连续的微分方程系统. 这种不连续的复杂性体现在力学现象上，就是系统会发生粘滑运动，而这种特殊的运动模式又让传统的分析方法显得无能为力. 为了更好地从数学上描述运动方程，以及从力学角度更好地解释各种运动情形的发生，本文采用Filippov系统的相关理论分析振动驱动系统. 这样做的好处是，Filippov滑移分岔理论能有效地解释振动驱动系统中的粘滑效应. 需要注意的是，Filippov系统中所谓的滑移运动，其实就对应了振动驱动系统粘滑运动中的黏滞状态.

采用Filippov系统分析的步骤是：建立系统一阶状态微分方程；确定滑移区域；具体讨论在滑移区域的边界是否发生滑移分岔. 当然，由于本文中采用的驱动是三相驱动，本身也存在不连续现象，所以这又增加了分析的复杂性. 基于此，文中通过分类讨论，将振动驱动系统中可能发生的各种滑移分岔现象进行了详细的罗列、归纳.

具体地，我们首先将(2)式由二阶非自治系统转化为一阶自治系统，即将t作为状态变量，引入变量替换，令

$x=\left( \begin{matrix} {{x}_{1}} \\ {{x}_{2}} \\ t \\ \end{matrix} \right)=\left( \begin{matrix} y \\ {\dot{y}} \\ t(\text{mod}T) \\ \end{matrix} \right)$

(10)

这里，${ x} \in { R}^3$，为系统的状态向量. 再引入用来定义不连续界面的函数$H({ x})=x_{2}$. 那么不连续界面

$\sum{=\{x\in {{R}^{3}}|H(x)={{x}_{2}}=0\}}$

(11)

可以将整个状态空间分为两个部分，分别为

$\left. \begin{align} & {{G}_{1}}=\{x\in {{R}^{3}}|H(x)>0 \\ & {{G}_{2}}=\{x\in {{R}^{3}}|H(x)<0 \\ \end{align} \right\}$

(12)

据此，可以将式(2)表示为

$\dot{x}=\left\{ \begin{array}{*{35}{l}} {{F}_{1}}(x,\mu ),x\in {{G}_{1}} \\ {{F}_{2}}(x,\mu ),x\in {{G}_{2}} \\ \end{array} \right.$

(13)

这里，有

${{F}_{1}}=\left( \begin{matrix} {{x}_{2}} \\ -\mu \ddot{\xi }-fg \\ 1 \\ \end{matrix} \right),{{F}_{2}}=\left( \begin{matrix} {{x}_{2}} \\ -\mu \ddot{\xi }+fg \\ 1 \\ \end{matrix} \right)$

(14)

如果在G₁和G₂两个子空间上，向量场F₁和F₂都局部地指向不连续界面Σ，这就意味着系统的流形在局部受到Σ面的约束. 其数学表述为

${{H}_{x}}{{F}_{1}}\le 0\le {{H}_{x}}F$

(15)

式中，H_x=0是 x的法向量.

分界面上满足这一条件的区域称为滑移区域(sliding region)，记为$\hat{\Sigma }$. 系统只有在滑移区域上，才会发生滑移运动，也就是出现黏滞状态.

根据Utkin等效控制理论，滑移区域上的运动由两个部分构成. 一个是F₁和F₂平均值，另一个是控制函数$\beta ({ x} )$，该函数用来将流形拉回到与滑移区域$\hat{\Sigma }$相切，这也就保证了系统的运动始终保持在滑移区域上，具体表达式如下

$\left. \begin{align} & {{F}_{s}}=\frac{{{F}_{1}}+{{F}_{2}}}{2}+\frac{{{F}_{2}}-{{F}_{1}}}{2}\beta (x) \\ & \beta (x)=-\frac{{{H}_{x}}{{F}_{1}}+{{H}_{x}}{{F}_{2}}}{{{H}_{x}}{{F}_{2}}-{{H}_{x}}{{F}_{1}}} \\ \end{align} \right\}$

(16)

将式(11)和式(14)代入式(16)，有

$\beta (x)=\frac{\mu \ddot{\xi }(t)}{fg}$

(17)

当$\beta ({ x }) = - 1$时，${ F}_s = { F}_1 $；当$\beta ({ x }) =1$时，${ F}_s = { F}_2 $. 所以滑移区域$\hat {\varSigma}$也可以表示成与式(15)等价的表达式

$\hat{\Sigma }=\{x\in \Sigma -1\le \beta (x)\le 1\}$

(18)

它的边界又可以表示为

$\partial {{\hat{\Sigma }}^{\pm }}=\{x\in \Sigma |\beta (x)=\pm 1\}$

(19)

2 运动分析

从式(7)、式(17)和式(18) 可以看出，自变量$w_i (i = 1,2,3)$的选取不仅定量上对系统的运动状态存在影响，而且定性上决定了运动的拓扑结构. 记$u_i \triangleq \mu w_i / fg$ (i= 1,2,3)，由式(17)可以看出，当u_i≤ 1时，$\vert \beta ({ x })\vert ≤1$，那么一个周期内，不连续界面在第i段(这里的第i段表示三相驱动3个时间段中的某一段，i的值可以取1,2或3)为滑移区域，系统有可能在这一时间段发生滑移运动(黏滞). 而当u_i时，$\vert \beta ({ x })\vert >1$，系统不可能在相应时间段上黏滞. 由上段分析结合式(2)和式(8)易知，当u₂ ≤1时，系统不能实现向前的稳态周期运动，所以本文只考虑u₂ > 1$时的情形.于是，由u₁和u₃与1的大小关系可知，系统可能发生4种情形的粘滑运动,见表 1.图 3中，灰色平面表示不连续界Σ，阴影部分为滑移区域$\hat {\varSigma}$，表示系统可能在区域内发生滑移运动. 滑移区域的边界$\partial \hat {\varSigma }^ + $和$\partial \hat{\varSigma }^ - $分别用实线、虚线表示. 注意观察式(7)，发现如果所有的$u_i e 1 (i = 1,2,3)$，那么三相驱动中并不存在一点严格满足式(19). 但是，当系统存在滑移区域的时候，其边界又是存在的.造成这一情况的原因是三相驱动本身的不连续性，式(7)中也可以看到在滑移区域的边界处，$\ddot {\xi}(t)$值并无定义，因此$\beta ({ x})$也无定义. 故而，我们需要在$\partial \hat {\varSigma }^ +$边界补充定义$\beta ({ x }) = 1$，在$\partial \hat {\varSigma }^ - $边界补充定义$\beta ({ x }) = -1$. 这样，当$\beta ({ x })$值由$\beta ({ x }) < - 1$变成$ - 1 ≤ \beta ({ x }) ≤1$时，边界处补充定义$\beta ({ x} ) = - 1$，表示从$\partial \hat {\varSigma }^ -$边界进入滑移区域，此时边界为虚线，反之表示系统从$\partial \hat {\varSigma }^ -$边界离开滑移区域，边界仍为虚线. 而当$\beta ({ x })$值由$\beta ({ x}) > 1$变成$ - 1 ≤ \beta({ x} ) ≤ 1$时，仿照前面分析，边界为实线，反之亦然. 以图 3(b)为例，可以由表 1判断，(t₁ ,t₂)内$\beta ({ x }) < - 1$，(t₂ ,T)内$0 < \beta ({ x }) ≤ 1$，(T,T + t₁ )内$\beta ({x }) > 1$，故而t₂ 时刻补充定义$\beta ({ x }) = - 1$，边界$\partial \hat {\varSigma }^ -$用虚线表示，T时刻补充定义$\beta ({ x }) = 1$，边界$\partial \hat {\varSigma }^{ + }$用实线表示.需要注意的是$(0,t_1 )$上$\beta ({ x} ) > 1$，(t₁ ,t₂ )内$\beta ({ x} ) < - 1$，虽然t₁ 时刻两边都不存在滑移区域，但是$\beta ({ x} )$符号不同，意味着状态向量场在t₁两边方向不同. t₁时刻左边，向量场由G₁空间指向G₂，t₁时刻右边，向量场由G₂空间指向G₁.为了体现出这种差异，t₁处用点划线表示.

2.1 情形1

此时，由表 1中的第一行数据，等价地得到μw₁ ≤ fg,μw₂ > fg,μw₃ ≤ fg.由式(17)和式(18)判定，可知不连续界面在(t₂ ,T + t₁ )的时间区域内为滑移区域，此时滑移区域示意图如图 3(a). 结合稳态运动的周期性观察可知，系统有两个可能存在滑移分岔的位置

$x_{1}^{*}=(x_{1}^{*},0,{{t}_{1}}),\quad x_{2}^{*}=(x_{1}^{*},0,{{t}_{2}})$

由于有刚体位移，可知，x₁的取值为任意值，这里都用$x_1^\ast$表示，下文中也类似. 运用Filippov判定条件来判断该点是否满足分岔点的定义和非退化条件. 由于

$\begin{array}{*{35}{l}} \underset{t\to t_{1}^{-},x\to x_{1}^{*}}{\mathop{\lim }}\,{{H}_{x}}{{F}_{s}}(x)=0 \\ \underset{t\to t_{1}^{+},x\to x_{1}^{*}}{\mathop{\lim }}\,{{H}_{x}}{{F}_{1}}(x)=\mu {{w}_{2}}-fg>0 \\ \end{array}$

所以可以补充定义$ { H}_x { F}_1 ({ x }_1^\ast ) = 0$，同样，补充定义$ { H}_x { F}_1({ x }_2^\ast ) = 0$. 这样，分岔点的定义得到满足

$\left. \begin{align} & H(x_{i}^{*})=0,{{H}_{x}}(x_{i}^{*})\ne \mathbf{1},i=1,2 \\ & {{H}_{x}}{{F}_{1}}(x_{1}^{*})=0,{{H}_{x}}{{F}_{1}}(x_{2}^{*})=0 \\ \end{align} \right\}$

(20)

对于非退化条件，借助式(8)，得到如下表达

$\left. \begin{align} & {{({{H}_{x}}{{F}_{1}})}_{x}}{{F}_{1}}(x_{1}^{*})=\mu ({{w}_{1}}+{{w}_{2}})\delta (0)\to +\infty >0 \\ & {{({{H}_{x}}{{F}_{1}})}_{x}}{{F}_{1}}(x_{2}^{*})=-\mu ({{w}_{2}}+{{w}_{3}})\delta (0)\to -\infty <0 \\ \end{align} \right\}$

(21)

注意到${ x}_{1}^{ \ast } $,${ x}_{ 2}^{ \ast } $所在的边界都是$\partial \hat {\varSigma }^ - $边界，所以根据文献^[18]，非退化条件(21)表明：在${ x}_{1}^{ \ast } $点处，系统可能发生crossing-sliding或者grazing-sliding分岔，而${x}_{2}^{ \ast } $点处，可能发生switching-sliding分岔.

分析该运动过程，首先假设x₂ (0) = 0. 考虑到(0,t₁ )上，分界面为滑移区域，初始状态x₂ (0)落在滑移区域上，因此系统发生滑移运动，即保持黏滞的状态到t₁，因而有x₂ (t₁ ) = 0. (t₁ ,t₂)上，系统加速度为正常数$\mu {{\text{w}}_{2}}-fg$，因此$x_2 (t_2 ) = (\mu w_2 - fg)\tau _2 > 0$. (t₂ ,T)上，如果系统的速度仍为正，其加速度为负常数$ - (\mu w_3 + fg)$. 据此，可以求出T时刻的速度：$x_2 (T) = (\mu w_2 - fg)\tau _2 - (\mu w_3 + fg)\tau _3$，可以证明该表达式的值为负(证明过程见附录1)，这就意味着系统将在(t₂,T)的某一段上进入滑移区域，并发生黏滞，因而实际的x₂ (T) =0. 考虑到稳态运动的周期性，可知前面的假设x₂ (0) =0正确，此处的分析符合系统稳态运动情况.

注意x₂ (t₁ )=0，x₂ (t₂ ) >0，也就是可能发生滑移分岔位置的速度符号确定，因而两点处均不会发生分岔. 此时，系统的运动状态示意图对应于图 5(a).

2.2 情形2

此时，$\mu w_1 > fg$,μw₂ > fg,μw₃ ≤ fg，不连续界面在(t₂,T)的时间域内为滑移区域，滑移区域示意图如图 3(b)所示，有两个可能的滑移分岔位置：${ x}_1^\ast = (x_1^\ast ,0,t_2 )$,${ x}_2^\ast = (x_1^\ast ,0,T)$.经验证，满足定义，并且非退化条件表明：系统在${ x}_{1}^{ \ast } $点处可能发生switching-sliding分岔，在${ x}_{ 2}^{ \ast }$点处可能发生crossing-sliding或者grazing-sliding分岔. 为了分析运动过程，假设x₂(0) = 0. (0,t₁ )上，系统的加速度为负常数 - μw₁ + fg，因此 $x_2(t_1) = (- \mu w_1 + fg)\tau _1 <0$. (t₁ ,t₂ )上，当系统速度仍为负时，加速度为正常数μw₂ + fg，可以证明经过$\Delta t_1 (\Delta t_1 < \tau _2 )$时间之后，系统速度回到0(证明过程见附录2)，但是考虑到(t₁ ,t₂ )区域不是滑移区域，因此速度会继续增大. 在($\Delta t_1,t_2 $)上，系统加速度为正常数μw₂ - fg，所以有$x_2 (t_2 ) = (\mu w_2 - fg)(\tau _2 - \Delta t_1 ) > 0$. (t₁ ,t₂ )上，当系统速度仍为正时，加速度为负常数 - (μw₃ + fg)，又可以证明，经过$\Delta t_2 (\Delta t_2 < \tau _3 )$时间之后，系统速度回到0(证明过程见附录3)，因而系统进入滑移区域，并保持黏滞状态至T，因此x₂ (T) = 0$. 又由稳态运动的周期性，可知假设x₂ (0) = 0正确，此处的分析符合稳态运动情况.

注意到x₂ (t₂ ) > 0，x₂ (T) = 0，也就是分岔位置处速度的符号确定，因而两点处均不会发生分岔. 系统的运动状态唯一，如图 5(b).

2.3 情形3

此时，μw₁ ≤ fg,μw₂ > fg,μw₃ > fg，滑移区域示意图如图 3(c)所示，有两个可能的滑移分岔位置：${ x}_1^\ast = (x_1^\ast ,0,t_1 )$,${ x}_2^\ast = (x_1^\ast ,0,T)$.经验证，满足定义，非退化条件显示：系统在${ x}_{1}^{\ast } $点处可能发生crossing-sliding或者grazing-sliding分岔，在${ x}_{2}^{\ast } $点处可能发生switching-sliding分岔. 为了分析运动过程，假设x₂ (t₁ ) = 0. (t₁ ,t₂ )上，系统的加速度为正常数μw₂ - fg，所以有$x_2 (t_2 ) = (\mu w_1 - fg)\tau _2 > 0$. (t₂ ,T)上，当系统的速度仍为正时，加速度为负常数$ - (\mu w_3 + fg)$，由此推算出的T时刻速度表达式，将其记为P₁ ，有

${{P}_{1}}=(\mu {{w}_{2}}-fg){{\tau }_{2}}-(\mu {{w}_{3}}+fg){{\tau }_{3}}$

(22)

虽然这里的表达式与情形1中x₂ (T)的表达式一样，但是需要注意的是，由于w_i (i =1,2,3)的取值范围与情形一中不同，故而得不到P₁ <0的结论，因而这里P₁的符号是不定的.根据分岔点的非退化条件，可知P₁ = 0意味着系统在${ x}_2^\ast $ 处发生switching-sliding 分岔. 当P₁ >0时，显然系统将在(T,t₁)中的某一时刻进入滑移区域，否则系统速度始终为正，所受摩擦力始终为负，由动量定理可知系统不能维持这种稳态运动.在滑移区域内，系统将保持这种黏滞状态至T + t₁ 时刻，因而x₂ (T + t₁ ) = 0.由稳态运动的周期性，可知这种情况下假设x₂ (t₁ ) = 0正确，因而分析符合实际，此时系统运动状态如图 5(c)所示.当 P₁ <0时，表明系统在(t₂ ,T)的中经过长度为$\dfrac{\mu w_2 - fg}{\mu w_3 + fg}\tau _2 $的时间后，速度减为0. 在到达T时刻之前，系统继续减速，但加速度变成 - μw₃ + fg，因此${{x}_{2}}(T)=({{\tau }_{3}}-\frac{\mu {{w}_{2}}-fg}{\mu {{w}_{3}}+fg}{{\tau }_{2}})(-\mu {{w}_{3}}+fg)<0$. 在(T,T + t₁ )的时间段上，加速度为正常数 - μw₁ +fg，将T + t₁ 时刻的速度记为P₂ ，有

$\begin{align} & {{P}_{2}}=({{\tau }_{3}}-\frac{\mu {{w}_{2}}-fg}{\mu {{w}_{3}}+fg}{{\tau }_{2}})(-\mu {{w}_{3}}+fg)+ \\ & (-\mu {{w}_{1}}+fg){{\tau }_{1}} \\ \end{align}$

(23)

根据非退化条件并结合运动，可知P₂=0意味着系统在 x₁^*处发生crossing-sliding分岔，而且由稳态运动的周期性，可知假设x₂ (t₁ ) =0正确，因此分析符合实际.

综上所述，P₁=0表示系统在x₂^*处发生switching-sliding分岔，分岔示意图如图 4(a)，其中轨迹B对应P₁ = 0，C对应P₁ <0，A对应P₁ > 0. 当P₁ <0时，P₂ = 0表示系统在x₁^*处发生crossing-sliding分岔，分岔示意图如图 4(b)，其中轨迹B对应P₂ = 0，C对应P₂ > 0，A对应P₂ <0.

因此，在这种情形(情形3)下，系统的运动状态有3种：① P₁ > 0，运动形态如图 5(c)；② P₁ <0，P₂ >0，运动形态如图 5(d)；③ P₁ <0，P₂ <0，运动形态如图 5(e).

2.4 情形4

此时，μ w₁ > fg,μ w₂ > fg,μ w₃ > fg，易知，系统不存在滑移区域，如图 3(d)所示.这种情形中，系统将不会发生黏滞. 因此，系统运动的拓扑结构唯一，且相对简单，系统做无黏滞的周期振动，运动示意图见图 5(f).

3 分岔

第2节中通过对系统运动分情形进行分析，观察到系统运动出现了滑移分岔现象，也得到了6种形式的粘滑运动.但是需要注意的是，这并不意味着这6种运动形式定性上各不相同.本节将先从滑移分岔的观点分析上述各运动属于何种滑移运动类型.然后再具体分析各种不同运动发生的条件，探究随着参数变化它们之间如何通过滑移分岔实现相互切换.

3.1 归纳

Filippov滑移分岔理论指出，右端不连续系统在不连续界面上可能会发生4种滑移分岔(crossing-sliding,grazing-sliding,switching-sliding,adding-sliding)^[29]，这4种分岔的名称也可以被借用来命名滑移运动的类型.利用图 3(a)，可以看到在图 5(a)中，系统从G₁ 空间进入滑移区域之后，并不进入G₂空间，而是保持滑移运动直到再次进入G₁空间.可以简单地理解成系统的运动状态轨迹擦过滑移区域，因而这种滑移运动称为grazing-sliding.同样可以看到，图 5(c)中的滑移运动也可以看成是一种grazing-sliding.这两种运动形态可以被归纳为同一种滑移运动类型.

观察图 5(b)，并结合图 3(b)可知，系统由G₁空间进入滑移区域后保持滑移运动至$\partial \hat {\varSigma }^{ +}$边界，然后离开滑移区域，进入G₂空间.可以简单地理解为系统的运动状态轨迹在经过一段滑移运动之后穿越了滑移区域，因而这种滑移运动称为crossing-sliding.

观察图 5(d)，并结合图 3(d)可知，系统由G₁空间直接进入G₂空间，然后又回到滑移区域. 由G₂空间进入滑移区域后，系统滑移运动至$\partial \hat {\varSigma }^ - $边界，最终回到G₁空间.可以简单地理解成，系统的状态轨迹先穿过了不连续界面，然后回到滑移区域继续运动，对于这样一个具有来回转换过程的运动，形象地称为switching-sliding.

图 5(e)和图 5(f)属于同一类运动，且都不发生滑移运动.

综上，系统存在4种定性不同的运动类型，其中3种运动中存在黏滞，1种无黏滞发生.

3.2 分岔图

为了更清晰地描述滑移分岔发生的条件，判断系统滑移运动的类型，引入符号$\tilde {\tau }_i (i = 1,2,3)$，令

${{\tilde{\tau }}_{i}}=\sqrt{\frac{fg}{2\mu L}}{{\tau }_{i}},i=1,2,3$

(24)

代入式(9)，经计算，可知$\tilde {\tau }_i $都是u_i的函数，且有如下的表达

$\left. \begin{array}{*{35}{l}} {{{\tilde{\tau }}}_{1}}=\sqrt{\frac{{{u}_{2}}}{{{u}_{1}}({{u}_{1}}+{{u}_{2}})}} \\ {{{\tilde{\tau }}}_{2}}=\sqrt{\frac{1}{{{u}_{2}}}}(\sqrt{\frac{{{u}_{3}}}{{{u}_{2}}+{{u}_{3}}}}+\sqrt{\frac{{{u}_{1}}}{{{u}_{1}}+{{u}_{2}}}}) \\ {{{\tilde{\tau }}}_{3}}=\sqrt{\frac{{{u}_{2}}}{{{u}_{3}}({{u}_{2}}+{{u}_{3}})}} \\ \end{array} \right\}$

(25)

根据式(22) $\sim$式(25)，可知

$\left. \begin{align} & {{P}_{1}}=\sqrt{2\mu Lfg}[({{u}_{2}}-1){{{\tilde{\tau }}}_{2}}-({{u}_{3}}+1){{{\tilde{\tau }}}_{3}}] \\ & {{P}_{2}}=\sqrt{2\mu Lfg}[({{{\tilde{\tau }}}_{3}}-\frac{{{u}_{2}}-1}{{{u}_{3}}+1}{{{\tilde{\tau }}}_{2}})(-{{u}_{3}}+1)+ \\ & (-{{u}_{1}}+1){{{\tilde{\tau }}}_{1}}] \\ \end{align} \right\}$

(26)

记

$\left. \begin{array}{*{35}{l}} {{{\tilde{P}}}_{1}}=({{u}_{2}}-1){{{\tilde{\tau }}}_{2}}-({{u}_{3}}+1){{{\tilde{\tau }}}_{3}} \\ {{{\tilde{P}}}_{2}}=({{{\tilde{\tau }}}_{3}}-\frac{{{u}_{2}}-1}{{{u}_{3}}+1}{{{\tilde{\tau }}}_{2}})(-{{u}_{3}}+1)+(-{{u}_{1}}+1){{{\tilde{\tau }}}_{1}} \\ \end{array} \right\}$

(27)

由第二节的分析，不难发现，如果知道一个系统u_i (i = 1,2,3)的取值，以及P₁和P₂的正负号情况，那么系统的运动状态可以确定. 从式(26)和式(27)可以看出，P₁和P₂的正负号情况分别与$\tilde {P}_1 $和$\tilde {P}_2 $一致，且最终只取决于u_i (i = 1,2,3).故而，系统的运动状态从根本上讲，取决于3个参数 u_i (i = 1,2,3)的取值.由前文分析可知，随着这3个参数变化，系统有如下4种类型的粘滑运动.

由3.1的分析可知，grading-sliding运动类型的发生有两种可能：一种u₁ ≤1,u₂ > 1,u₃ ≤ 1，运动状态如图 5(a)所示，另一种$u_1 ≤ 1,u_2> 1,u_3 > 1,\tilde {P}_1 >0$，运动状态如图 5(c)所示. Crossing-sliding发生的条件为$u_1 > 1,u_2 > 1,u_3≤ 1$，运动状态如图 5(b)所示. Switching-sliding发生的条件为$u_1≤ 1,u_2 > 1,u_3 > 1,\tilde {P}_1 <0,\tilde {P}_2 >0$，运动状态如图 5(d)所示. 无黏滞运动的发生有两种可能：一种 $u_1 ≤1,u_2 > 1,u_3 > 1,\tilde {P}_1 <0,\tilde {P}_2 <0$，运动状态如图 5(e)所示，另一种$u_1 > 1,u_2 > 1,u_3 >1$，运动状态如图 5(f)所示. 为了便于查看，这4种粘滑运动发生的条件在表 2中予以罗列.

根据表 2，可以绘制出当u₂给定后，振动驱动系统在u₁-u₃参数平面的分岔图. 如图 6所示，本文分别选取了两种具有代表性的情况,图 6 (a) u₂=3，图 6 (b) u₂=5. 图中，点划线表示$\tilde {P}_1 =0$，虚线表示$\tilde {P}_2 = 0$，直线u₁=1和u₃=1都用实线表示. 图中， R_i (i = 1,2,3,4)表示发生不同粘滑运动的区域，下标i的值与表 2中的序号对应.对比图 6(a)和图 6(b)，发现图 6(a)中不存在$\tilde {P}_1 = 0$这条线，而图 6 (b)中存在；相应地，图 6 (a)中的R₁是1片区域，而图 6 (b)中的R₁ 是2片区域. 由此可知，当u₂的值较小(通过数值计算，观察到大约u₂ <3.5)时，$\tilde {P}_1 <0$恒成立，看不到所谓的switching-sliding分岔，而当u₂的值较大(大约u₂>3.5)时，switching-sliding分岔可能会发生.

当u₂较大时，运动分类情况较为典型也较为丰富. 以图 5(b)为例，从左下角的R₁出发逆时针绕一圈，研究各分岔线或分界线对运动行为的影响. 从R₁出发，当u₁的值增加时，系统跨过u₁=1分界线进入R₂ (箭头I)，运动类型由grazing-sliding (图 5(a))变成crossing-sliding (图 5(b)). 这是因为随着u₁增加，系统在第一段时间上的黏滞区域逐步地变为非黏滞区域，无法再保证系统不发生后退运动. 从R₂出发，当u₃的值增加时，系统跨过u₃ = 1分界线进入R₄ (箭头Ⅱ)，运动类型由crossing-sliding(图 5(b))变成无黏滞运动(图 5(f)). 这是因为，随着u₃的增加，仅存的黏滞区域也逐步变成非黏滞区域，因而系统将不再黏滞. 注意到$\tilde {P}_1 = 0$与$\tilde {P}_2 =0$两条曲线在u₁相交，这是因为当u₁=1时，系统在第一段上加速度为0(见式(2))，所以$\tilde {P}_1 =0$成立，也就意味着$\tilde {P}_2 = 0$成立. 由此可知，从R₄出发，系统有两种选择，一种是跨过u₁ = 1进入R₁ (箭头Ⅲ)，运动类型再次变成grazing-sliding (图 5(c)). 这是因为随着u₁减小，系统在第一段时间的非黏滞区域又变成黏滞区域，因而保证系统不会后退，而且，这个区域中$\tilde {P}_1 >0$，系统在第三段时间内会由于惯性保持速度为正，不会向后运动. 另外一种，系统从R₄出发，经过$\tilde {P}_2 =0$分岔线进入R₃区域(箭头IV)，运动类型变成switching-sliding (图 5(d))，这是因为跨过虚线，$\tilde {P}_2$由负变正，意味着系统在第一段时间上速度由负值开始逐渐增大，且一定会在t₁时刻之前，达到速度为0.考虑到第一段时间区间上为滑移区域，因而这会引起switching-sliding运动. 从R₁出发，跨过$\tilde {P}_1 = 0$分岔线，系统也会进入R₃(箭头V)，这是因为点划线左边$\tilde {P}_1 <0$，意味着x₂ (0) <0，而第一段时间区间为滑移区域，因而会发生switching-sliding运动. 最后，从R₃出发，随着u₃减小，系统会回到R₁区域(箭头VI).这是因为，系统在第三段时间上的非粘滑区域变成了粘滑区域，再次保证了系统不会向后运动.

4 数值验证

由第3节可知，决定系统分岔特性的变量是驱动参数u_i (i = 1,2,3)，为了验证分岔结果是否正确，以及运动特性是否与理论分析一致，有必要对其进行数值计算的验证.本节中物理参数固定，唯一改变的是驱动参数. 这里选定L =1 m,μ = 0.333 3,g =10 m/s²,f=0.2，通过改变控制驱动参数w_i (i = 1,2,3)，可以得到所有类型的运动特性图形. 本节中取u₂ =5，这样系统的粘滑运动类型可以由图 6(b)直观地看出. 图 7通过数值计算，得到了系统在不同区域的速度时间图.通过比较表 3和图 7，发现理论分析出的运动类型和数值计算的运动形态吻合. 通过对比图 5和图 7，发现运动形态的理论分析和数值计算结果吻合. 由此，验证了本文中理论分析的正确性.

5 结论

本文研究了受三相驱动的振动驱动系统在水平面上的直线运动. 不同于以往的研究，本文重点关注了系统在各向同性干摩擦下的运动. 这种环境与物理实际比较吻合，研究工作具有潜在的实用价值.

含库伦干摩擦的振动驱动系统是一种Filippov系统，因而本文采用滑移分岔的理论分析了其粘滑运动. 在只考虑系统向前运动的情况下，可以将运动分成4种基本情形. 对4种情形分类讨论，并加以归纳，最终得到4种不同的粘滑运动. 借助理论分析得到的分岔条件，我们给出了在驱动参数w_i(i = 1,2,3)变化时，系统运动的分岔图. 通过数值模拟，运动分类以及各种运动形态的正确性都得到了验证.

通过本文的滑移分岔分析，可以看到，三相振动驱动系统在各向同性干摩擦环境下具有丰富的动力学行为，值得进一步的实验验证和探究. 此外，从结构上改善和提升振动驱动系统，使其更好地适应各向同性干摩擦环境，也是未来值得思考的问题.