分形接触干摩擦模型及在围带叶片受迫响应分析中的应用

邱恒斌，徐自力; Qiu Hengbin; Xu Zili

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

分形接触干摩擦模型及在围带叶片受迫响应分析中的应用

- ORCID：
邱恒斌 ¹
✉
- ORCID：
徐自力 ²

1. 北京强度环境研究所，北京 100076； 2. 西安交通大学机械结构强度与振动国家重点实验室，西安 710049

最近更新：2021-08-30

DOI：10.6052/1672-6553-2020-073

摘要

为了研究含有围带接触面叶片的非线性振动响应，基于三维干摩擦微滑移模型，发展了一种分形接触干摩擦微滑移模型，用于计算考虑接触面形貌的摩擦力.在该模型中，摩擦接触面被离散成一系列接触单元，每个接触单元由一组接触点对来表征该接触单元的粘滞、滑移、分离的摩擦运动状态.采用分形几何模拟接触表面的形貌，基于分形理论和赫兹接触理论，建立接触面粗糙度、正压力、弹性模量、泊松比等参数与接触刚度和摩擦系数的关系.接触面摩擦力由接触刚度、摩擦系数和接触点对的相对位移确定.采用该发展的模型预测了真实围带叶片的受迫振动响应，研究了接触面形貌和初始正压力对围带叶片共振响应的影响.研究结果表明，本文模型能考虑接触面形貌对围带叶片非线性振动响应的影响；接触面形貌改变时，接触刚度和摩擦系数会发生变化，从而影响接触面摩擦力；在该模型中，接触刚度随着正压力的增加而增加，随着粗糙度的增加而减小；摩擦系数随着正压力增加而增加，随着粗糙度的增加，先增大后减小；接触面形貌和初始正压力对围带叶片受迫振动响应有显著影响.

关键词

分形接触; 干摩擦模型; 接触刚度; 摩擦系数; 非线性振动

2020-03-09 收到第1稿， 2020-04-19 收到修改稿.

引言

在透平机械中，叶片系统经常采用围带结构增加叶片间的干摩擦阻尼来达到减振的目的.为了研究围带干摩擦阻尼叶片的振动特性，建立一个高效准确的干摩擦模型是十分必要的.

在已有干摩擦阻尼系统的研究中，微滑移模型由于能较好的模拟实际摩擦接触行为而被广泛采用.Iwan^［

1］采用线性弹簧和库仑摩擦元组合，在接触面分布多个接触点对，提出串联和并联微滑移模型.Griffin^{［参考文献 2

百度学术}2］基于库仑摩擦定律，采用恒定摩擦系数，并采用线性弹簧表征切向接触刚度，研究了恒定正压力下的一维摩擦相对运动.Yang^{［参考文献 3

百度学术}3］提出了二维的干摩擦接触模型，采用两个无质量的弹簧表征切向和法向的接触刚度，研究了变正压力下的二维摩擦相对运动.Chen，Yang和Menq^{［参考文献 4

百度学术}4］在二维模型基础上提出了三维的摩擦接触模型.Petrov^{［参考文献 5

百度学术}5，6］采用三维摩擦模型研究了多种形式的干摩擦阻尼结构的减振特性.在已有的微滑移模型中，通常采用线性弹簧模拟接触面发生相对滑移之前的弹性变形，用库仑摩擦元模拟接触面发生滑移时的摩擦力.

接触刚度和摩擦系数是计算摩擦力的两个重要参数，在先前大部分摩擦接触模型中，摩擦系数和接触刚度大都采用经验值.实际上，摩擦系数和接触刚度均会随着接触正压力，接触面形貌，弹性模量，泊松比等参数的变化而改变.在计算接触刚度和摩擦系数过程中，接触正压力，弹性模量，泊松比等参数较容易获得，但是准确的模拟接触面形貌是较为困难的.

在以往的接触研究中，大多采用GW（Greenwood-Williamson模型）和分形理论模拟接触面形貌.Chang等^［

7］基于GW接触模型，提出金属表面静摩擦系数的计算模型，得出摩擦系数与正压力的关系.尤晋闽等^{［参考文献 8

百度学术}8］研究了分形接触表面的摩擦模型，得出了静摩擦系数与接触面材料参数、接触面分形参数和正压力的关系.Jones ^{［参考文献 9

百度学术}9］利用GW模型研究了法向载荷、切向载荷、表面形貌等对切向刚度的影响.Jiang等^{［参考文献 10

百度学术}10］采用分形几何模拟接触面形貌，研究了机械结合面的接触刚度，理论计算结果与实验结果相接近.GW模型获得的表面特征参数受限于分辨率及采样长度等，不能准确反映出粗糙表面的全部信息^{［参考文献 11

百度学术}11］.粗糙表面轮廓具有自相似性，且与尺度无关，因此可以由分形几何来表征，利用表面分形特性建立接触模型，可使表面接触的分析结果具有确定性和唯一性^{［参考文献 11

百度学术}11，12］.

本文基于三维干摩擦微滑移模型，发展了一种分形接触干摩擦微滑移模型，用于计算考虑接触面形貌的摩擦力.采用分形几何模拟接触表面的形貌，基于分形理论和赫兹接触理论，建立接触面粗糙度、正压力、弹性模量、泊松比等参数与接触刚度和摩擦系数的关系.由接触刚度、摩擦系数和接触点对的相对位移计算摩擦力.采用该发展的模型研究真实围带叶片的受迫振动响应.

1 分形接触

1.1　分形接触表面

分形具有自相似性和尺寸独立性，使其能够描述具有多尺度特性的粗糙表面.具有分形特征的三维粗糙表面轮廓曲线可由W-M（Weierstrass-Mandelbrot）函数来模拟，该函数具有处处连续，处处不可导，具有无限细节特征^［

13］.该数学表达式为：

\begin{array}{l} z (x) = L {(\frac{G}{L})}^{(D - 1)} {(\frac{l n γ}{M})}^{1 / 2} \sum_{m = 1}^{M} \sum_{n = 0}^{n_{m a x}} γ^{(D - 2) n} \times {c o s ϕ_{m, n} - \\ c o s [\frac{2 π γ^{n} (x^{2} + y^{2})^{1 / 2}}{L} \times c o s (a r c t a n (\frac{y}{x}) - \frac{π m}{M}) \\ + ϕ_{m, n}} \end{array}

（1）

式中： $z (x)$ 为粗糙表面轮廓的高度； $x, y$ 为粗糙表面轮廓的位移坐标； $G$ 为反映表面轮廓 $z$ 的特征尺度系数， $G$ 越大越粗糙； $D$ 为表面轮廓的分形维数，表征表面轮廓的不规则性和复杂程度，且 $1 < D < 2$ ； $n$ 为空间频率序数为，最低频率序数为0，最高频率序数为 $n_{m a x} = i n t [l o g (L / L_{g}) / l o g γ]$ ， $L_{g}$ 为截止长度，通常近似为材料的原子间距离；M为曲面褶皱的重叠数； $γ^{n}$ 表示轮廓曲线的空间频率，一般 $γ$ =1.5比较符合高频谱密度和相位随机的情况； $c o s ϕ_{m, n}$ 是随机相位，取值范围为［0， $2 π$ ］.图1为D=1.3时模拟的分形接触表面.

图1 三维分形接触表面模拟

Fig.1 Simulated 3D fractal surfaces

1.2　接触刚度

围带接触面的接触刚度源于接触表面的一系列微凸体，研究表明^［

13］可将两个粗糙表面的接触等效为一个粗糙表面和一个光滑的刚性平面接触的如图2.

图2 粗糙表面微凸体与刚性平面的接触

Fig.2 Contact between an asperity of composite rough surface and a rigid plane

则该系统的等效弹性模量为：

E^{*} = [(1 - ν_{1}^{2}) / E_{1} + (1 - ν_{2}^{2}) / E_{2}]^{- 1}

（2）

式中 $v_{1}$ 、 $v_{2}$ 和 $E_{1}$ 、 $E_{2}$ 分别为两个接触微凸体的泊松比和弹性模量.微凸体变形量由式（1）W-M分形函数^［

13］决定：

δ = 2 G^{(D - 1)} {(l n γ)}^{1 / 2} {(2 r')}^{(2 - D)}

（3）

Majumdar和Bhushan^［

12］研究表明接触点的面积分布密度函数与海洋面上的岛屿面积分布密度函数相似，得到切割面积

a'

的统计学分布函数为：

n (a') = \frac{D}{2} (a_{l}^{'})^{D / 2} {(a')}^{- (D + 2) / 2}

（4）

式中， $a' = π {(r')}^{2}$ 是微接触的切割面积， $a_{l}'$ 是最大微接触切割面积.微凸体的曲率半径为 $R^{2} = {(R - δ)}^{2} + {(r')}^{2}$ ，忽略二阶小量可写为 ${(r')}^{2} = 2 R δ$ ，代入到式（3）可得：

R = \frac{{(a')}^{D / 2}}{2^{(4 - D)} π^{D / 2} G^{(D - 1)} {(l n γ)}^{1 / 2}}

（5）

分形维数D、分形参数G与接触面粗糙度的关系可近似为^［

14］：

D = 1.515 / R a^{0.088}

，

G = 10^{(- 8.259 / R a^{0.088})}

.Kogut等^{［参考文献 15

百度学术}15］将微凸体的变形可分为四个阶段：完全弹性、弹塑性和完全塑性变形.即当

δ < δ_{c}

时微凸体处于完全弹性状态；当

δ_{c} \leq δ \leq 6 δ_{c}

时微凸体处于第一弹塑性状态；当

6 δ_{c} < δ \leq 110 δ_{c}

时微凸体处于第二弹塑性状态；当

δ > 110 δ_{c}

时微凸体处于完全塑性状态.微凸体处于各个接触状态的接触载荷如表1所示^{［参考文献 15

百度学术}15］.

表1 接触变形及对应的接触载荷

Table 1 Contact deformation and corresponding contact load

Deformation	Contact load
$δ < δ_{c}$	$F_{e} = 4 / 3 E^{*} R^{1 / 2} δ^{3 / 2}$
$δ_{c} \leq δ \leq 6 δ_{c}$	$F_{e p 1} = 2 / 3 K H π R δ_{c} \times 1.03 {(δ / δ_{c})}^{1.425}$
$6 δ_{c} < δ \leq 110 δ_{c}$	$F_{n e p 2} = 2 / 3 K H π R δ_{c} \times 1.4 {(δ / δ_{c})}^{1.263}$
$δ > δ_{c}$	$F_{p} = 2 π R δ H$

微凸体接触面由完全弹性到弹塑性变形的临界变形量为：

δ_{c} = {(\frac{π K H}{2 E^{*}})}^{2} R

（6）

式中， $H$ 是材料的硬度， $K$ 是材料的硬度系数， $K = 0.454 + 0.41 v$ ， $v$ 是软材料的泊松比.将式（3）和式（5）代入式（6），微凸体由完全弹性到弹塑性变形的临界切割面积可表示为：

a_{c}' = {[2^{7 - 2 D} π^{D - 1} G^{2 D - 2} l n γ {(\frac{E^{*}}{K H})}^{2}]}^{1 / (D - 1)}

（7）

微凸体的法向接触刚度：

k_{n} = \frac{d F_{e}}{d δ} + \frac{d F_{e p 1}}{d δ} + \frac{d F_{e p 2}}{d δ}

（8）

根据文献［

16］得到微凸体的切向接触刚度：

k_{τ} = \frac{4 G r}{2 - v}

（9）

整个接触面的接触刚度可认为是处于完全弹性、第一和第二弹塑性区的单个微凸体接触刚度的总合.因此整个接触面的法向接触刚度可表达为：

\begin{array}{l} K_{n} = \int_{a'_{c}}^{a'_{l}} k_{n e} n (a') d a' + \int_{6^{1 / (1 - D)} a'_{c}}^{a'_{c}} k_{n e p 1} n (a') d a \\ + \int_{110^{1 / (1 - D)} a'_{c}}^{6^{1 / (1 - D)} a'_{c}} k_{n e p 2} n (a') d a' \end{array}

（10）

同样的，整个接触面的切向接触刚度可表达为：

K_{τ} = \int_{110^{1 / (1 - D)} a'_{c}}^{a'_{l}} k_{τ} n (a') d a'

（11）

切向和法向接触刚度随正压力和粗糙度的变化如图3和图4所示.由图可知，接触刚度随着正压力的增加而增加，随着粗糙度的增加而减小.接触表面越粗糙意味着接触表面越尖锐，接触表面微凸体在受到接触载荷时越容易发生塑性变形，因此接触刚度会减小.

图3 接触刚度随正压力的变化

Fig.3 Contact stiffness vs normal pressure

图4 接触刚度随粗糙度的变化

Fig.4 Contact stiffness vs roughness

1.3　摩擦系数

当微凸体同时受到法向载荷F和切向载荷T共同作用时，微凸体表面的应力场目前还没有完整的解析解，现有文献大多基于一定的假设条件对表达式进行简化.本文采用文献^［

8］中提出的假设，认为微凸体的屈服发生在微凸体表面接触点的边缘，接触点边缘的应力为：

T = \frac{8 σ_{s} a}{π (6 - 3 v)} + \frac{8 (2 v - 1)}{π (6 - 3 v)} F

（12）

式中 $σ_{s}$ 为较软材料的屈服强度.在法向载荷作用下，接触界面上处于塑性变形的微凸体由于受到局部接触载荷作用而发生塑性流动，将不能继续承受切向载荷^［

17］.因此在计算静摩擦力只包括前两个状态的微凸体，当接触面处于临界滑动状态时，接触面将达到完全屈服，此时整个接触面所有微凸体所承受的切向载荷即为最大静摩擦力

T_{t}

：

T_{t} = \int_{a'_{c}}^{a'_{l}} T n (a') d a' + \int_{6^{1 / (1 - D)} a'_{c}}^{a'_{c}} T n (a') d a'

（13）

接触面的法向载荷为：

\begin{array}{l} F_{n} = \int_{a'_{c}}^{a'_{l}} F_{e} n (a') d a' + \int_{6^{1 / (1 - D)} a'_{c}}^{a'_{c}} F_{e p 1} n (a') d a' + \\ \int_{110^{1 / (1 - D)} a'_{c}}^{6^{1 / (1 - D)} a'_{c}} F_{e p 2} n (a') d a + \int_{a'_{s}}^{110^{1 / (1 - D)} a'_{c}} F_{p} n (a') d a' \end{array}

（14）

因此摩擦系数可以由式（15）得到：

μ = \frac{T_{t}}{F_{n}}

（15）

接触面摩擦系数随正压力和粗糙度的变化如图5和图6所示.由图可知，摩擦系数随着正压力的增加而增加，随着粗糙度的增加而先增大后减小.

图5 摩擦系数随正压力的变化

Fig.5 Friction coefficient vs normal pressure

图6 摩擦系数随粗糙度的变化

Fig.6 Contact stiffness vs roughness

2 干摩擦模型

叶片振动时，围带接触面会发生三维相对运动，产生摩擦约束力.为了模拟实际摩擦接触过程中的微滑移现象，将围带接触区域被离散成一系列接触单元，每个接触单元由一组接触点对来表征该接触单元在摩擦接触运动过程中粘滞、滑移和分离的摩擦运动状态.

接触面上单个接触点对的模型如图7所示.接触点对的接触刚度由式（10）和式（11）确定，摩擦系数由式（15）确定.A，B两点分别为互相接触的两个平面上对应的接触点，用两个无质量的线性弹簧表征接触刚度.弹簧 $K_{t}$ 表征接触面的切向接触刚度，弹簧 $K_{n}$ 表征接触面间的法向接触刚度.当接触面做三维相对运动时，假设点A与面1始终保持粘滞，而点B会沿着面2产生粘滞或滑移的运动状态.在面2的接触面上建立局部坐标系，在局部坐标系中，点A和点B的位移分别定义为 $(x_{k}^{A}, y_{k}^{A}, z_{k}^{A})$ 和 $(x_{k}^{B}, y_{k}^{B}, z_{k}^{B})$ ，点A和点B之间的相对位移定义为

法向相对位移：

z_{k} = z_{k}^{A} - z_{k}^{B}

（16）

切向相对位移：

u_{k} = [\begin{matrix} x_{k}^{A} - x_{k}^{B} \\ y_{k}^{A} - y_{k}^{B} \end{matrix}]

（17）

图7 单个摩擦接触点对模型

Fig.7 The model of a friction pair

接触面间的法向接触正压力由接触面初始正压力 $n_{0}$ 、法向接触刚度 $k_{n}$ 和法向相对位移 $z_{k}$ 确定，可由下式来计算：

f_{n}^{k} = \{\begin{array}{r} \begin{matrix} n_{0} + k_{n} z_{k} & w h e n z_{k} \geq - n_{0} / k_{n} \end{matrix} \\ 0 w h e n z_{k} < - n_{0} / k_{n} \end{array}

（18）

若接触面初始正压力 $n_{0}$ 为负，表示接触面存在初始间隙； $f_{n}^{k}$ 等于零，则表示接触面在此时刻发生分离，接触点摩擦力为零.接触面摩擦约束力就可由点A和点B间的相对位移计算获得.切向摩擦力 $f_{t}^{k}$ 定义为：

f_{t}^{k} = {f_{x}^{k}, f_{y}^{k}}

（19）

定义点A与点B之间的切向拉力 $T_{k}$ 为：

T_{k} = \sqrt[]{u_{k}^{T} (K_{t})^{T} K_{t} u_{k}}

（20）

如果 $T_{k}$ 小于临界摩擦力 $μ f_{n}^{k}$ ，则点B保持不动，若 $T_{k}$ 大于 $μ f_{n}^{k}$ ，那么点B将向点A滑动直至拉力 $T_{k}$ 等于临界摩擦力 $μ f_{n}^{k}$ .则切向摩擦力可以表示为：

f_{t}^{k} = \{\begin{array}{l} K_{t} u_{k} w h e n T_{k} < μ f_{n}^{k} \\ μ f_{n}^{k} u_{k} / |u_{k}| w h e n T_{k} \geq μ f_{n}^{k} \end{array}

（21）

3 围带叶片响应

某真实围带叶片的有限元模型如图8所示，采用该发展的分形接触三维干摩擦模型和高阶谐波平衡法^［

18］预测该围带叶片的受迫振动响应.对叶片轮盘施加固定约束，在叶片轴向上施加

p (t) = 50 s i n ω t

的激振力，接触面刚度及摩擦系数采用分形接触模型计算，选取前5阶谐波和前5阶模态进行响应计算.为了研究接触面粗糙度对围带叶片减振效果的影响，本文选取了表面粗糙度为Ra=1.6μm和Ra=2.2μm的两种情况，计算两种表面粗糙度在不同初始正压力下叶片系统的振动响应，两种粗糙度下叶片系统的幅频曲线如图9和图10所示.

图8 围带叶片的有限元模型

Fig.8 The finite element model of a shroud blade

图9 $R_{a} = 2.2 μ m$ 时的幅频曲线

Fig.9 Amplitude-frequency curves at Ra=2.2μm

图10 $R_{a} = 1.6 μ m$ 时的幅频曲线

Fig.10 Amplitude-frequency curves at Ra=1.6μm

由图9和图10可知，当接触面粗糙度不同时，叶片的幅频曲线、共振幅值和共振频率均发生变化，说明叶片的振动响应对接触面粗糙度较为敏感.就本文选取参数计算所得结果，在相同的初始正压力下，表面粗糙度1.6μm计算的共振幅值较表面粗糙度2.2μm的小.同时发现该真实围带叶片的响应没有出现摩擦锁死的动力学现象，这是因为实际机组为了在振动过程中达到减振的目的，会选取合适的初始正压力，防止出现摩擦锁死现象，这与机组实际运行的情况相符.

4 结论

本文基于三维干摩擦微滑移模型，发展了一种分形接触干摩擦微滑移模型，用于计算考虑接触面形貌的摩擦力.在该模型中，摩擦接触面被离散成一系列接触单元，每个接触单元由一组接触点对来表征该接触单元的粘滞、滑移、分离的摩擦运动状态.采用分形几何模拟接触面形貌，基于分形理论和赫兹接触理论，建立接触面粗糙度、正压力、弹性模量、泊松比等参数与接触刚度和摩擦系数的关系.由接触刚度、摩擦系数和接触点对的相对位移计算摩擦力.

结果表明：该分形接触三维干摩擦模型能较好地模拟接触面形貌，有效地考虑接触面形貌、正压力等参数对摩擦力的影响.其中接触刚度随着正压力的增加而增加，随着粗糙度的增加而减小；摩擦系数随着正压力增加而增加，随着粗糙度的增加，先增大后减小；存在最优的初始正压力使叶片共振幅值最小，共振频率随着初始正压力的增加而增加；接触面形貌对围带叶片受迫振动响应有显著影响.

参考文献

Iwan W D. On a class of models for the yielding behavior of continuous and composite systems. Journal of Applied Mechanics， 1967，89： 612~617 [百度学术]

Griffin J H， Friction damping of resonant stresses in gasturbine engine airfoils. Journal of Engineering for Gas Turbines and Power，1980， 102： 329~333 [百度学术]

Yang B D， Chu M L， Meng C H. Stick-slip-Separation analysis and non-linear stiffness and damping characterization of friction contacts having variable normal load.Journal of Sound and Vibration，1998， 210（2）：461~481 [百度学术]

Chen J J， Yang B D， Menq C H. Periodic forced response of structures having three-dimensional frictional constraints. Journal of Sound and Vibration， 2000， 229 （4）：775~792 [百度学术]

Petrov E P. A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks. Journal of Turbomachinery，126， 2004：175~183 [百度学术]

Petrov E P. Explicit finite element models of friction dampers in forced response analysis of bladed disks. Journal of Engineering for Gas Turbines and Power， 2008， 130（1）：1~11 [百度学术]

Chang W R， Etsion I， Bogy D B. Static friction coefficient model for metallic rough surfaces. ASME Journal of Tribology， 1988， 110（1）：57~ 63 [百度学术]

You J M，Chen T N. A static friction model for the contact of fractal surfaces， Proceedings of the Insititution of Mechanical Engineering. Part J： Journal of Engineering Tribology. 2010， 224： 513~518 [百度学术]

Jones R E. A Greenwood-Williamsion model of small-scale friction. Journal of Applied Mechanics， 2007， 74（1）：31~40 [百度学术]

Jiang S， Zheng Y， Zhu H. A contact stiffness model of machined plane joint based on fractal theory. Journal of Tribology， 2010，132：1~7 [百度学术]

Majumdar A， Tien C L. Fractal characterization and simulation of rough surfaces. Wear， 1990，136：313~327 [百度学术]

Majumdar A， Bhushan B. Fractal model of elastic–plastic contact between rough surfaces. Journal of Tribology，1991， 113（1）：1~11 [百度学术]

Yan W， Komvopoulos K. Contact analysis of elastic-plastic fractal surfaces. Journal of Applied Physics，1998，84（7）：3617~3624 [百度学术]

方兵．精密数控机床及其典型结合面理论建模与实验研究［博士学位论文］．长春：吉林大学，2012（Fang B. Theoretical modeling and experimental study of precision CNC machine tools and its typical interface ［Ph.D Thesis］.Changchun： Jilin University， 2012 （in Chinese）） [百度学术]

Kogut L，Etsion I. Elastic-plastic contact analysis of a sphere and a rigid flat. ASME Journal of Applied Mechanics， 2002，69：657~662 [百度学术]

Johnson K L. Contact mechanics.Cambridge：Cambridge University Press，1985 [百度学术]

盛选禹，雒建斌，温诗铸. 基于分形接触的静摩擦系数预测. 中国机械工程，1998，9（7）：16~18 [百度学术]

Sheng X Y， Luo J B， Wen S Z.Static friction coefficient model based on fractal contact. China Mechanical Engineering， 1998，9（7）：16~18（in Chinese） [百度学术]

邱恒斌，徐自力，刘雅琳，等.一种求解含围带阻尼成圈叶片振动响应的高效方法. 西安交通大学学报，2016，50（11）：1~6 [百度学术]

Qiu H B， Xu Z L， Liu Y L， et al. An efficient method for solving vibration response of continuous covered blades with damped shroud. Journal of Xi’an Jiaotong Univeristy， 2016，50（11）：1~6（in Chinese） [百度学术]

分形接触干摩擦模型及在围带叶片受迫响应分析中的应用

摘要

关键词

引言