库仑阻尼的一种有效处理：光顺化（SMOOTHING）

Rainyboy

中间由于毕业的原因，离开了振动论坛好久，现在有时间了，我回来啦！
分享一个干摩擦非线性计算中的技巧吧……

在分析干摩擦阻尼器时，常库仑（Coulomb）阻尼模型，即用符号函数来表示摩擦力的大小与接触面相对速度的关系：     
摩擦力 = -1×正压力×摩擦系数×sign(相对速度)， if 相对速度 != 0
摩擦力 = [-1×正压力×摩擦系数,  正压力×摩擦系数]， if 相对速度 == 0
其中，sign()是符号函数，其定义为：
       if x >0, sign(x) = 1  
       if x <0, sign(x) = -1

上述速度-力的关系是分段连续的，这给采用库仑干摩擦假设的模型的求解带来了一些困难，例如在数值积分法中，为了找到接触状态（阻滞，滑移）变化的点，需要精确而合适地设置时间步，降低了求解的效率。在这个文献中，可以找到详细的评述：
Jerry H. Griffin, 1990, A Review of Friction Damping of Turbine Blade Vibration[J].

我之前处理过一个这样的二自由度库仑摩擦模型：
图1

其中采用的就是上述分段连续的摩擦力模型，在采用非线性Newmark方法分析其时序响应时，颇费了一番心思，最后搞了一个自动拆分时间步长的方法才算是较准确地捕捉到接触状态发生变化的时刻，其时序响应（至稳态）大概是：
图2
其中最后一个稳定周期是：
图3
图中各线的意义是：
      蓝色- Id 的响应
      绿色- 相对运动
      红色- Im 的响应

最近仔细看了C. Duan 的系列文章，看到了一种对上述分段连续的摩擦力模型的光顺化处理，即：
       摩擦力 = -1×tanh(Sigema × 相对速度)
其中Sigema 只是一个认为确定的系数，可以用来调整光顺化处理之后的曲线与处理前的接近程度。

Sigma越大，则光顺后的曲线越接近经典库仑模型，如下图：
图4


于是，取Sigma = 3, 重新写了一个按照光顺的曲线计算上述二自由度模型的程序，取与图2和图3一致的参数，计算得：
图5
最后一个稳定周期为：
图6


对比图6和图3，可以发现二者的结果非常接近，其中一些细节的差异我想计算过类似问题的同志都能体会出来，我就不赘述了。然而图5，6的计算时间却远少于图2，3（大概少了80%），主要原因就是不用在程序中反复迭代和确认当前时间步的接触状态，与不需要判断接触状态了——因为对于光顺后的摩擦力模型，实际上已经不存在接触状态这一概念了，再小的摩擦力都会对应一个相对运动速度，这以情形在经典的库仑干摩擦模型中实际上应当是归于‘阻滞’（Block）状态的。这些差异对于人来说也许是巨大的，然而，对于计算程序来说，一个很小的速度和速度为零，又有什么区别呢？

最后，给出C. Duan的这篇文章：
Chengwu Duan, Rajendra Singh, 2005, Super-harmonics in a torsional system with dry friction path subject to harmonic excitation under a mean torque, Journal of Sound and Vibration.

这篇文章的亮点不是这个光顺化的处理，但是我确实是在这篇文章中第一次看到这么巧妙的处理。当然，这篇文章提出的算法如果不采用光顺处理后的速度-摩擦力关系恐怕也是行不通的，有兴趣的同志可以搜来看一看。

ChaChing

学习了!
有一想法, 将相关数据补齐, 方便有兴趣者试试学习!

qiandashuai

mark，等我有空的时候来讨论讨论

Rainyboy

回复 2 # ChaChing 的帖子

好的。

图1的无量纲动力学方程可以写为：   
   

其中，
        

在图2和3，以及4和5中采用的计算参数是一致的，为：
Mcr = 3.0
Ud = 0.1
Omega = 0.95
kesi(机械阻尼比) = 0.5%

ChaChing

回复 4 # Rainyboy 的帖子

谢谢提供
汗, 都快忘记学术上常用无量纲! 所以1F的图2也是无量纲?
另图2接近400s处好像有一突变点, 不知该如何说明?

4F的东西看了下, 真有点看不明白, 可否稍加说明?
水平有限, 从没玩这类型的模拟, 想了解清楚, 有空练习下

Rainyboy

回复 5 # ChaChing 的帖子


呵呵，是首先要列出‘普通’（带量纲）的方程，再进行无量纲化。附件就是其中推导的详细过程哈。

帖子里的图2～图5都是用无量纲模型进行的计算。

至于那一个‘突变’，不知怎么解释，在计算中也不常见，因为我比较关注响应稳定的幅值，中间过程受程序控制，可能是误差所致吧。

程瑞岩

来看看支持一下。

lex9317

今天试验了下，好像就是用tanh函数代替sign函数，我这么理解对么？确实比sign快多了。
f=zeros(2,1)
f(1)=x(2)
f(2)=-0.5*9.8*tanh(8*x(2))-200*x(1)/10
%f(2)=-0.5*9.8*sign(x(2))-200*x(1)/10
tspan=[0:0.1:8]
x0=[5;0]
[t,x]=ode23('func',tspan,x0)
plot(t,x(:,1))

Rainyboy

回复 8 # lex9317 的帖子

是的！但是你的代码似乎没贴全？

Seventy721

支持这样的讨论。谢谢楼主分享好内容。

lex9317

嗯，把定义函数的那句话忘了粘贴了，呵呵，function f=func(t,x)这个是首句

lex9317

觉得即使是这样的处理方式在多自由度系统里面，求解也是非常非常慢的吧？比如200自由度？那ansys之类的软件是怎么处理这摩擦种东西的呢？

Rainyboy

回复 12 # lex9317 的帖子

ANSYS里面的接触单元的算法比这里介绍了要复杂许多，作为非线性 问题还只能瞬态求解，所以计算起来非常的慢，还有收敛的问题。
一般要得到频响曲线的话是不会拿ANSYS去硬算的，一般会用谐波平衡方法构建频域求解的程序。谐波平衡(Harmonic Balance)只是一个求解框架，具体的处理有时/频转换(Time/Frequencty Alternating)方法，线性化(Linearise)方法，增量谐波平衡(Incremental Harmonic Balance)等方法等等，原来在帝国理工现在去了Sussex的E.P.Petrov做这个做得比较多，应该在国际上现在是一流的。你感兴趣的话可以去搜索上面的中英文关键词。

lex9317

谢谢RAINBOY啊，看来问题还有点复杂哦

lex9317

前两天我建立的那个一维模型，结果和用sign(x)符号函数做的一样，但是有个问题，就是当我初始状态时候，弹簧被拉长了很小，只有0.2，那么在初始弹力还不如摩擦力大，可是算出来依然会运动啊。我想是因为，模型是按照滑动来做的，没有考虑静摩擦，任何把静摩擦这个限制条件加进去呢？