一种考虑界面不确定性和时变性的摩擦振动噪声预测方法与流程

本发明涉及一种摩擦振动噪声的预测方法，特别的是一种考虑界面不确定性和时变性的摩擦振动噪声预测方法。

背景技术：

摩擦振动噪声是由两个相互接触的物体在发生相对运动时引起其自激振动并向外辐射声压的一种现象。在众多中机械系统，摩擦引起的振动噪声严重影响了机械系统的精度、可靠性以及使用寿命，同时引起的噪声污染对人类、社会和生态环境造成了极大的危害。但由于其随机性、难预测性、影响因素多等特点，至今仍未有一种方法能完全消除摩擦引起的振动噪声。如今，众多研究从结构改变、材料选择和工况调整等角度来尝试改善和消除摩擦振动噪声，但摩擦学实验存在时间成本大、经济成本高等缺点，因此有必要提出一种摩擦振动噪声的预测方法，减少研究的时间和经济成本。

现有技术中，复特征值分析作为一种摩擦振动噪声的预测方法，得到了广泛的使用，该方法能够大致的预测出系统可能产生的摩擦振动的频率，但作为一种线性方法，由于摩擦过程中大量非线性因素的存在，该方法往往会导致结果的过度预测或欠预测。另一方面，摩擦接触界面的磨损随着摩擦过程必然出现并且随着时间的推移而改变，并且已有大量文献证明摩擦接触界面形貌特征对摩擦振动噪声有重要的影响。而现有的摩擦振动噪声实时预测方法通常只能在表面完整未受磨损时进行预测，并不能跟随时间的推移而进行跟踪预测。

同时，摩擦界面的磨损作为一种存在大量非线性因素的现象，现有的摩擦振动预测方法尚未将其作为影响因素列入模拟范围。同时，现有的摩擦接触界面的磨损计算方法通常以线速度、接触压力、节点位置和温度等因素作为自变量进行计算得到磨损量及方向，但通常难以对材料局部的分布、硬度等因素进行考虑，这些不确定性的因素会对接触面上各点的磨损深度和方向造成影响。如果能在振动噪声的预测中考虑磨损表面的变化，那么就能得到更为准确的预测结果，为后续的工作提供依据，奠定基础。

技术实现要素：

针对现有技术中振动噪声的预测中未考虑磨损表面的变化，从而影响预测结果的准确性的问题，本发明提供一种考虑界面不确定性和时变性的摩擦振动噪声预测方法，其目的在于：将表面磨损造成的摩擦界面随时间的变化加入模拟范围，能够更加准确的预测摩擦过程中的振动噪声现象，精确地得到其频率以及其强弱关系。在摩擦振动噪声预测方法领域，尚无与本发明类似的方法与流程。

本发明采用的技术方案如下：

一种考虑界面不确定性和时变性的摩擦振动噪声预测方法，包括如下步骤：

步骤1.针对实际的摩擦系统简化并建立有限元模型；

步骤2.设定计算时长t，并将t平均离散化为n段，初始的循环次数记为0；

步骤3.对步骤1建立的初始状态下的有限元模型进行复特征值分析，对摩擦系统中的初始表面的摩擦振动噪声结果进行预测；

步骤4.对前步骤所得有限元模型进行时域动态分析计算；

步骤5.选取磨损公式，在主程序中调用子程序对步骤4的动态分析计算结果进行信息提取，提取信息的种类根据所述磨损公式所需值来确定；

步骤6.将步骤5提取的信息带入磨损公式计算得到磨损方向及深度；

步骤7.对步骤6计算得到的磨损方向及深度进行蒙特卡洛模拟得到新的磨损方向及深度；

步骤8.根据步骤7的得到的新的磨损方向及深度更新步骤4所采用的有限元模型中的摩擦界面的节点坐标，生成新的有限元模型；

步骤9.对步骤8得到的新的有限元模型进行复特征值分析，对摩擦振动进行预测；

步骤10.使循环次数加1，在循环次数达到n的情况下结束摩擦振动噪声的实时预测；在循环次数未达到n的情况下，将步骤8得到的新的有限原模型代入步骤4中，重复进行步骤4至步骤9的操作。

优选地，步骤1中，对所述摩擦系统建立有限元模型的范围包括摩擦界面和摩擦界面周边的零件。

优选地，步骤2中，时域上的离散化方法为平均处理，离散化后的循环次数n，每次计算步长时间t0和总时长t满足以下关系：

其中每次计算步长时间t0应远小于总时长t。

优选地，采用的复特征值分析方法满足以下公式：

(λ²[m]+λ[c]+kns){φ}＝0，

式中，[m]为系统的质量矩阵，[c]为系统的阻尼矩阵，[kns]为系统的非对称刚度矩阵，λ为特征值，{φ}为其对应的特征向量；

得到的复特征值计算结果表示为：

λ＝α+jω，

式中，ω为特征值的虚部，代表振动频率，α为特征值的实部，代表该振动频率的稳定性；实部大于0，说明对应的振动频率不稳定，实部越大说明系统在对应的振动频率下产生不稳定振动噪声的倾向越大。

优选地，步骤6中，有限元模型中摩擦界面上任一节点上磨损方向的确定方法为该节点与相邻节点所组成的向量的法向量之和，磨损方向w用公式表示为：

其中，磨损节点表示为m0，磨损节点附近有j个相邻节点，分别表示为m1，m2,…,mj，磨损节点与相邻节点所组成的向量分别为所述向量所对应的法向量为n1，n2,…,nj，以与接触压力法向夹角小于180为正方向。

优选地，磨损深度计算公式的通用形式表述为：

(h,θ)＝f(s,p,t,…)，

其中，h为磨损深度，θ为磨损方向，s，p，t及省略表示的参数为磨损公式中的变量。磨损深度的计算方法与影响变量可由用户自行定义，可选择但不局限于现有的磨损公式，用户可根据实际情况编写或改变磨损公式及其内在变量。

优选地，步骤5中，从结果文件中提取信息的种类根据磨损公式所需值来确定，包括节点坐标、温度场和接触压力。

优选地，步骤7中，将步骤6中得到磨损深度和方向的值进行随机化处理，生成附近的随机值作为新的磨损深度和方向；随机磨损深度和方向的获得方法为用步骤6中得到的磨损深度和方向带进行monte-carlo处理，计算公式表示为：

其中，

var[]为样本方差,e[]为均值/原始值，h为步骤6中得到的磨损深度，θ为步骤6中得到的磨损方向，hnew为新的磨损深度，θnew为新的磨损方向。

优选地，步骤8中新的有限元模型的摩擦界面的节点坐标由以下计算公式得到：

(x1,y1,z1)＝(x0,y0,z0)+(xh,yh,zh)

其中(x1,y1,z1)为新的有限元模型的摩擦界面的节点坐标，(x0,y0,z0)为步骤4所采用的有限元模型中的摩擦界面的节点坐标，(xh,yh,zh)为各个方向上的磨损深度，为该节点上磨损深度的向量表述形式，并有：

优选地，步骤8中有限元模型的更新方法通过将原命令流文件中的接触面节点坐标进行替换实现。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1)本发明在进行摩擦振动噪声预测的过程中考虑了摩擦过程中表面形貌的变化，使得预测结果能更接近真实情况；

2)本发明在原有的静态预测方法上考虑了摩擦振动噪声的时变特性，通过随时间改变摩擦接触面，能对不同磨损情况下的摩擦振动噪声进行预测，改善了原复特征值方法只能预测初始状态下摩擦振动噪声的缺陷；

3)原复特征值预测方法为线性分析，经常会导致过预测和欠预测的情况发生，在引入摩擦界面不确定性磨损这一非线性项后，有利于改善过预测和欠预测现象。

4)摩擦实验具有时间金钱成本大的特点，利用本发明的预测方法可以大量减少实验数量和研究成本，其预测结果也可以作为结构修饰改善摩擦振动噪声的评判标准，为后续的优化工作等提供了极大的便利；

5)本发明具有可控变量多的特点，用户可根据自身需求选择合适的设置，如计算时长、磨损公式、时域动态分析的求解器和随机化处理函数等；

6)本发明在简单的初始设置后为全自动循环计算，使用方便，对用户友好，同时不失准确性。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1为本发明实施例1采用的摩擦试验台及其有限元模型；

图2为本发明中考虑摩擦界面不确定性和时变特性的摩擦振动噪声实时预测方法的流程图；

图3为本发明实施例1中实验初始的法向振动加速度时域信号；

图4为本发明实施例1中实验末尾的法向振动加速度时域信号；

图5为本发明实施例1中实验初始的振动加速度频谱；

图6为本发明实施例1中实验末尾的振动加速度频谱；

图7为本发明在实施例1中，实验与预测的振动噪声的频域结果比较；

图8为本发明在实施例1中，预测的振动噪声强度变化结果；

图9为本发明实施例2模拟结束时摩擦接触面的磨损情况。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本发明提出的是一种考虑界面不确定性和时变性的摩擦振动噪声预测方法，下面参照附图1至图9与文字，以具体实例进行详细的说明讲解。值得注意的是，此处所描述的具体实施仅用以解释本发明，并不用于限定本发明。

本方法的实施方法主要分为两部分，基于有限元方法的有限元软件和执行循环主程序的软件。本实例中选用的有限元软件为abaqus，主程序在matlab中执行，主程序中的读取结果文件功能由python实现，在主程序中调用。

本发明的实施过程为：在时域上进行离散化，将总分析时长划分为均匀的n个时间间隔；在空间上进行离散化，建立离散的有限元模型；对原模型采用复特征值分析方法，得到初始的振动噪声预测结果；采用时域动态分析法，对结果文件采用从外部进行特征参数提取，如线速度、接触应力、温度等；结合选定的磨损公式与随机函数对磨损深度和磨损方向进行不确定性模拟；根据结果确定磨损节点的位移变化量；并依据磨损结果更新有限元模型；判断是否达到预设总分析时长，若未达到则重新采用复特征值分析预测摩擦振动噪声并循环进行，若已达到总分析时长则完成分析。

下面通过具体的实施例对本发明的技术方案做进一步的说明。

实施例1

图1所示为本实例的摩擦试验台示意图，本实例中的对摩件为销-盘结构，摩擦盘以90rad/s的速度按所示方向旋转，摩擦块在夹具与法向力的作用下始终与盘面贴合，摩擦块材料为汽车制动器上所使用的复合材料，摩擦盘所用材料为蠕墨铸铁。相较于盘式样，块式样的弹性模量更小，且在摩擦过程中盘面的接触位置始终在改变而摩擦块的接触面始终保持接触状态，因此块上的磨损较盘面要明显得多，因此本实施例将只考虑块上的磨损，忽略摩擦面上磨损对振动噪声的影响。同时由于系统结构对摩擦振动噪声有重大的影响，因此在有限元建模时应该考虑整体的实验装置而不仅是两对磨件。实验中的动摩擦系数为0.25。本实施例将以此为工况对本发明进行阐述。

步骤(1)、设定参数，计算总时长t，每次计算时长t0，此处磨损公式选用的是自定义的摩擦公式：h＝k·p·r·v^0.5，其中h为磨损深度，k为磨损系数，p为节点上法向压力，r为摩擦半径，v为相对滑动速度。磨损方向由相邻节点位置决定，随机化的方法采用monte-carlo模拟，分布方差设为均值的0.001；

步骤(2)、在有限元软件abaqus中针对实验台建立有限元模型，并分别生成用于复特征值分析的inp文件和用于动态分析的inp文件。其中动态分析时长为t0；

步骤(3)、在matlab中调用abaqus对模型进行复特征值分析，完成得到dat文件；

步骤(4)、在matlab中调用python程序读取dat文件，通过搜索关键字的方式获得复特征值结果，提取其中实部大于0的量并保存这些值；

步骤(5)、在matlab中调用abaqus对模型进行动态分析，得到odb文件；

步骤(6)、在matlab中调python程序读取odb结果文件，根据磨损公式需要提取其接触面节点坐标、节点上的接触压力p；

步骤(7)、每个节点上的接触压力为p，通过节点坐标可计算得到单位时间滑动距离s，s的计算公式可表示为：其中xi,yi为接触面上的节点i的坐标，(xd,yd)为摩擦盘的旋转中心坐标，ω为摩擦盘的旋转角速度；

步骤(8)、根据上一步得到的p、s可以得到每一个节点上的磨损深度hi＝k·si·pi。由于是面-面接触形式，摩损方向皆为垂直于接触表面，以平行于摩擦块切入面的方向为x轴正向，垂直于摩擦块切入面的方向为y轴正向，垂直于盘面的方向作为z轴正向，用向量的形式可将磨损方向表示为θ＝(0,0,1)；

步骤(9)、以上一步的磨损深度与角度作为分布均值，进行monte-carlo模拟得到新的磨损深度以及方向：hnew＝pmh(h)，θnew＝pmθ(θ)。

步骤(10)、计算新的接触面节点坐标，计算公式为：

其中，(x0,y0,z0)为原节点坐标，(x1,y1,z1)为新节点坐标；

步骤(11)、根据所得新接触面节点坐标，通过字符串替代的方式改写复特征值分析和动态分析的inp文件，实现模型的更新；

步骤(12)、重复进行步骤(3)-(10)直至达到预设时间。

由图7可见，在未考虑摩擦界面的变化下预测频率与实际的实验所测的频率之差为27.2hz,当采用本方法进行预测后，对摩擦振动的频率预测误差由27.2hz降低至23.5hz，在频率上提高了预测的精确度，

图3，4，6表明，随着摩擦界面的变化，摩擦振动的强度逐渐上升。复特征值实部是一种能反映系统不稳定摩擦振动噪声的强度的参数，当采用本方法考虑摩擦界面的变化后，如图8所示，系统的复特征值实部随着摩擦界面的变化而上升，与实验所示一致，因此该方法能同时在频率和强度两个方向上增强了对摩擦振动噪声的预测能力。

实施例2

本实施例中所阐述的为对摩副为球-面样品时的预测方法，对摩球以一定速度在对摩面上往复移动，对摩球在夹具与法向力的作用下始终与摩擦面贴合。同样的，由于在摩擦过程中摩擦面的接触位置始终在改变而对摩球面始终保持接触状态，因此球样品上的磨损较盘面要大得多，因此忽略摩擦面上磨损对振动噪声的影响。具体实施步骤如下：

步骤1、设定参数，计算总时长t，每次计算时长t0，此处磨损公式选用的是自定义的摩擦公式：h＝k·p·r·v^0.5，其中h为磨损深度，k为磨损系数，p为节点上法向压力，r为摩擦半径，v为相对滑动速度。磨损方向由相邻节点位置决定，球-面接触下磨损方向如图9所示。

步骤2、在有限元软件abaqus中针对实验台建立有限元模型，其中对球-面的接触点做一定的简化。简化方式为在接触点位置以一极小平面代替曲面，方便后续分析避免模拟过程中不收敛的情况。

步骤3、并分别生成用于复特征值分析的inp文件和用于动态分析的inp文件。其中动态分析时长为t0。

步骤4、在matlab中调用abaqus对模型进行复特征值分析，完成得到dat文件。

步骤5、在matlab中调用python程序读取dat文件，通过搜索关键字的方式获得复特征值结果，提取其中实部大于0的量并保存这些值。

步骤6、在matlab中调用abaqus对模型进行动态分析，得到odb文件。

步骤7、在matlab中调用数据提取子程序读取odb结果文件，根据磨损公式需要提取其接触面节点坐标、节点上的接触压力p。

步骤8、每个节点上的接触压力为p，通过节点坐标可计算得到单位时间滑动距离s，s的计算公式可表示为：其中xi,yi为接触面上的节点i的坐标，(xd,yd)为摩擦盘的旋转中心坐标，ω为摩擦盘的旋转角速度。

步骤9、根据上一步得到的p、s可以得到每一个节点上的磨损深度hi＝k·si·pi。并根据图所示确定其磨损方向；

步骤10、以上一步的磨损深度与角度作为分布均值，以上一步的磨损深度与角度作为分布均值，进行monte-carlo模拟得到新的磨损深度以及方向：hnew＝pmh(h)，θnew＝pmθ(θ)。

步骤11、计算新的接触面节点坐标，计算公式为：

其中，(x0,y0,z0)为原节点坐标，(x1,y1,z1)为新节点坐标。

步骤12、根据所得新接触面节点坐标，通过字符串替代的方式改写复特征值分析和动态分析的inp文件，实现模型的更新。

步骤13、重复进行步骤3-10直至达到预设时间。

以上所述实施例仅表达了本申请的具体实施方式，其描述较为具体和详细，但并不能因此而理解为对本申请保护范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请技术方案构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。