ABSTRACT

In this paper, the approximate periodic solution of a vibration system subjected to dry friction is studied. The periodic solution of this system approximately represents the Fourier series, the differential equation is transformed into nonlinear algebraic equations by Galerkin method, and the nonlinear algebraic equations are solved by Broyden method. Finally, numerical integration results are used to verify the effectiveness of the method.

Keywords:Dry Friction, Fourier Series, Galerkin Method, Broyden Method

应用于干摩擦振动系统中的Galerkin法

栗鹏

同济大学航空航天与力学学院，上海

收稿日期：2019年11月7日；录用日期：2019年11月22日；发布日期：2019年11月29日

摘 要

本文主要研究一个受干摩擦作用的振动系统的周期解的近似求解。此系统的周期解近似地表示Fourier级数，并通过Galerkin法把微分方程转换成非线性代数方程组，利用Broyden法来对此非线性代数方程组求解。最后数值积分的结果用来验证该方法的有效性。

关键词 :干摩擦，Fourier级数，Galerkin法，Broyden法

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

众所周知，干摩擦广泛地存在于各类机械系统中，比如刹车系统、离合器、轮胎运动和弓弦乐器等 [1]。干摩擦力在这些系统中可以发挥有利的作用，但也可以引起系统一些异常的振动，颤振等等。因此研究干摩擦系统的动力学行为具有理论和实际意义。

关于干摩擦动力学系统的研究方法主要有：近似解析法 [2]，数值求解和实验等等。近似解析法主要有平均法、谐波平衡法和增量谐波平衡法。Zimmermann等 [3] 利用平均法研究受干摩擦力作用的两自由度振动驱动系统。Duan等 [4] 先把干摩擦力用一个光滑的函数表示，然后利用谐波平衡法对一类扭转系统进行动力学分析。Pierre等 [5] 利用增量谐波平衡法来研究一个干摩擦阻尼振子的动力学响应。以上的方法在处理干摩擦动力学系统时都存在局限性，比如系统在运动过程中存在粘–滑状态转换时，上述方法就不适用了，而且粘–滑状态转换是干摩擦振动系统的重要现象，关于该类周期解的求解也是至关重要的。Wang等 [6] [7] 提出一个修正的增量谐波平衡法，其本质就是通过Galerkin过程把微分方程转化为非线性的代数方程组，并用Boryden法代替Newton法，并把此方法运用到非光滑的传动系统中。因此本文就是利用这种观点，我们不考虑增量过程，把解设为Fourier级数，通过Galerkin过程，微分方程变为非线性的代数方程组，最后利用Broyden法对代数方程进行求解。该类方法适用范围更广，可以用于求解弱非线性系统的周期解；也可以用于求解强非线性系统的周期解，比如系统中存在粘–滑状态转换的情况。

2. 模型的建立

本文中所考虑的干摩擦振动系统如图1所示，一个质量为m的物体，通过一个刚度系数k，粘性系数c的弹簧与一个固定端连接，而且该物体与地面之间存在干摩擦力，摩擦系数为 $\mu$，一个正弦激励u作用在该物体上，设 $u=F\sin \omega t$，F是激励振幅， $\omega$ 是激励频率。

图1. 干摩擦振动系统的示意图

根据牛顿第二定律得到：

$m\ddot{x}+c\stackrel{\dot{}}{x}+kx=F\sin \omega t+f$ (1)

其中f是干摩擦力，形式如下：

$f=-\mu mg\text{sign}\left(\stackrel{\dot{}}{x}\right)$ (2)

化简得到：

$m\ddot{x}+c\stackrel{\dot{}}{x}+kx+\mu mg\text{sign}\left(\stackrel{\dot{}}{x}\right)=F\sin \omega t$ (3)

为了减少系统的参数，引进如下参数：

${\omega }_0=\sqrt{\frac{k}{m}}$，$\xi =\frac{c{\omega }_0}{2k}$，$\tau =\omega t$，$\Omega =\frac{\omega }{{\omega }_0}$，$\alpha =\frac{F}{k}$，$\beta =\frac{\mu mg}{k}$

则方程(3)变为：

${\Omega }^2\ddot{x}+2\Omega \xi \stackrel{\dot{}}{x}+x+\beta \text{sign}\left(\stackrel{\dot{}}{x}\right)=\alpha \sin \tau$ (4)

3. Galerkin法

对于强非线性的动力学系统，Galerkin法可以作为对系统的周期解求解的一种近似方法。

设该系统的周期解的Fourier级数为：

$x=a_0+{\displaystyle \underset{n=1}{\overset{N}{\sum }}\left(a_n\cos \left(n\tau \right)+b_n\sin \left(n\tau \right)\right)}=V\left(\tau \right)A$ (5)

其中 $V\left(\tau \right)=\left[1,\cos \tau ,\cdots ,\cos N\tau ,\sin \tau ,\cdots ,\sin N\tau \right]$，$A={\left[a_0,a_1,\cdots ,a_N,b_1,\cdots ,b_N\right]}^{\text{T}}$ ；N是Fourier级数的阶数， $A$ 为解的Fourier级数的系数所组成的向量， ${\left[\right]}^{\text{T}}$ 表示向量的转置。

$\stackrel{\dot{}}{x}={\left(V\left(\tau \right)A\right)}^{\prime }=V\left(\tau \right)CA$ (6)

$\ddot{x}={\left(V\left(\tau \right)CA\right)}^{\prime }=V\left(\tau \right)DA$ (7)

其中 $C=\left[\begin{array}{ccc}0& 0_{1\times N}& 0_{1\times N}\\ 0_{N\times 1}& 0_{N\times N}& -S_{N\times N}\\ 0_{N\times 1}& S_{N\times N}& 0_{N\times N}\end{array}\right]$，$D=\left[\begin{array}{ccc}0& 0_{1\times N}& 0_{1\times N}\\ 0_{N\times 1}& -W_{N\times N}& 0_{N\times N}\\ 0_{N\times 1}& 0_{N\times N}& -W_{N\times N}\end{array}\right]$，$S_{N\times N}=\left[\begin{array}{ccc}1& & 0\\ & \ddots & \\ 0& & N\end{array}\right]$，$W_{N\times N}=\left[\begin{array}{ccc}1& & 0\\ & \ddots & \\ 0& & N^2\end{array}\right]$。

把(6)和(7)代入方程(4)得到：

$r={\Omega }^{2}V\left(\tau \right)DA+2\Omega \xi V\left(\tau \right)CA+V\left(\tau \right)A+\beta \text{sign}\left(V\left(\tau \right)CA\right)-\alpha \sin \tau =0$ (8)

由Galerkin过程得到系统的残差方程：

$R={\displaystyle {\int }_0^{2\pi }\left({\Omega }^{2}V\left(\tau \right)DA+2\Omega \xi V\left(\tau \right)CA+V\left(\tau \right)A+\beta \text{sign}\left(V\left(\tau \right)CA\right)-\alpha \sin \tau \right)V{\left(\tau \right)}^{\text{T}}\text{d}\tau }=0$ (9)

对于一般的Galerkin法来说，求解残差方程一般利用解析的方法，进而得到Fourier级数的系数。但是对于非线性方程组(9)来说，由于符号sign函数的存在，方程(9)中的 $a_0,a_1,\cdots ,a_N,b_1,\cdots ,b_N$ 并不能解析求出，而且方程(9)对于系数 $a_0,a_1,\cdots ,a_N,b_1,\cdots ,b_N$ 的导数也不能显示表示，尤其系统中存在stick运动的情形，数值求解中的Newton法也就不再适用于此方程的求解。因此可以利用Broyden法来对方程(9)进行数值求解，迭代过程中系统的Jacobian矩阵可以数值给出。其迭代过程如下：

第一步：给出近似的初值 $A_0$ 和Jacobian矩阵 ${{\rm B}}_0$，并由此初值开始，对系统进行迭代，直到求出满足精度的解；

第二步：假设方程已经进行了第i次迭代，得到方程的解为 $A_i$，把 $A_i$ 代入(5)中，并把(5)在一个周

期内进行离散化，即把 $\left[0,2\pi \right]$ 等分为H个小区间，则 ${\tau }_j=j\Delta \tau$，$j=1,2,\cdots ,H$，其中 $\Delta \tau =\frac{2\pi }{H}$，$H\ge 2N+1$，

则对应的向量为：

$\bar{x}={\left[x\left({\tau }_1\right),x\left({\tau }_2\right),\cdots ,x\left({\tau }_H\right)\right]}^{\text{T}}=\Gamma A_i$ (10)

其中 $\Gamma =\left[\begin{array}{ccccccc}1& \cos {\tau }_1& \cdots & \cos N{\tau }_1& \sin {\tau }_1& \cdots & \sin N{\tau }_1\\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 1& \cos {\tau }_H& \cdots & \cos N{\tau }_H& \sin {\tau }_H& \cdots & \sin N{\tau }_H\end{array}\right]$，把(10)代入(8)得到：

$\bar{r}={\left[r\left({\tau }_1\right),r\left({\tau }_2\right),\cdots ,r\left({\tau }_H\right)\right]}^{\text{T}}=\Gamma r_F$，$r_F={\left[r_0^i,\text{2}{r}_1^i,\cdots ,\text{2}{r}_N^i,\text{2}{q}_1^i,\cdots ,\text{2}{q}_N^i\right]}^{\text{T}}$。

第三步：对于r利用离散的Fourier变换得到： $r_F={\Gamma }^{-1}\bar{r}$，其中

${\Gamma }^{-1}=\frac{1}{H}\left(\begin{array}{ccc}1& \cdots & 1\\ 2\cos {\tau }_1& \cdots & 2\cos \left({\tau }_H\right)\\ ⋮& \ddots & ⋮\\ 2\cos \left(N{\tau }_1\right)& \cdots & 2\cos \left(N{\tau }_H\right)\\ 2\sin {\tau }_1& \cdots & 2\sin \left({\tau }_H\right)\\ ⋮& \ddots & ⋮\\ 2\sin \left(N{\tau }_1\right)& \cdots & 2\sin \left(N{\tau }_H\right)\end{array}\right)$，

则：

$R_i=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }r{T}^{\text{T}}\text{d}t}={\left[r_0^i,r_1^i,\cdots ,r_N^i,q_1^i,\cdots ,q_N^i\right]}^{\text{T}}$ (11)

第四步： $\Delta A_i=-B_i^{-1}{R}_i$，则：

$A_{i+1}=A_i+\Delta A_i$ (12)

若 $\Vert \Delta A_i\Vert <\delta$ 成立，则迭代结束，这里的 $\Vert \Delta A_i\Vert$ 是当前增量的模， $\delta$ 为容许误差。

若不满足，则下一次迭代的近似Jacobian矩阵为：

$B_{i+1}=B_i+\frac{\left(R_{i+1}-R_i-B_i\Delta A_i\right)\Delta A_i^{\text{T}}}{\Delta A_i^{\text{T}}\Delta A_i}$ (13)

一直重复上面的过程，直至收敛。

4. 数值验证

本文中考虑的系统具有强非线性和非光滑的特性，为了保证系统解的精确性，这里Fourier级数的阶数取为40。并且利用数值积分的方法来验证该方法的有效性。

当 $\Omega =1$，$\xi =0.1$，$\beta =1$，$\alpha =1.6$ 时，系统的稳态速度与时间的关系如图2所示，Galerkin表示本文的方法，*表示数值积分的方法。可以看到，系统在运动过程中不存在粘–滑状态转换，仅有滑动状态，两种方法的计算结果可以有比较好的吻合。

当 $\Omega =1$，$\xi =0.1$，$\beta =1$，$\alpha =1.4$ 时，系统的稳态速度与时间的关系如图3所示。可以看到：系统在这组参数下，其运动过程中存在粘–滑状态的转换，两种方法计算的结果也有很好的吻合。

图2. 系统的稳态速度与时间的关系

图3. 系统的稳态速度与时间的关系

5. 结论

本文主要研究了Galerkin法在一类干摩擦振动系统中的应用。其主要思想是：把系统的周期解设为Fourier级数，Galerkin过程把系统的微分方程转化为非线性的代数方程组，然后利用Broyden法求解非线性代数方程组，并利用数值积分的结果验证该方法的有效性。但是Boryden法有其自身的局限性，就是该方法强烈地依赖于给定的初始值和初始的Jacobain矩阵。因此寻找更有效的算法来求解此类代数方程组是将来更进一步的研究内容。

文章引用

栗 鹏. 应用于干摩擦振动系统中的Galerkin法
Galerkin Method Applied to the Vibration System with Dry Friction[J]. 声学与振动, 2019, 07(04): 155-160. https://doi.org/10.12677/OJAV.2019.74017

参考文献