航空、航天和机械等工程系统中的复杂结构大多是子部件通过螺栓法兰、过盈配合等连接形式装配而成.振动环境下，连接界面受到周期性外载荷作用，在法向可能出现预紧力振荡、界面分离和碰撞，在切向可能发生黏着、摩擦和滑动等行为.这些行为往往具有非线性、跨尺度、不确定性等特点^[1-2].研究表明，连接结构接触界面间的复杂力学行为不仅会造成结构的能量耗散，产生结构阻尼，而且会导致刚度非线性现象，从而影响整个结构的动力学特性^[3].

针对上述诸多非线性现象，学者们建立了许多接触模型来复现连接界面间的复杂力学行为.早期，一般采用库伦摩擦模型来建模连接界面的摩擦行为，但由于没有考虑界面微观黏滑行为的影响，计算的连接刚度和界面阻尼特性误差较大^[4-5].随着接触力学和摩擦学等学科的发展，许多有物理意义的连接模型被提出来，比如Valanis模型^[6]、LuGre刷毛模型^[7]和Iwan模型^[8].其中，由于Iwan模型更为系统和完善地描述了界面的微观黏滑过程，且模型的参数识别过程简单，因而被广泛关注^[9].

最早，Iwan模型采用无限多个Jenkins单元串-并联的方式用以描述土木工程结构的弹塑性迟滞行为.Jenkins单元是一个理想弹塑性模型，它是由一个线性弹簧和一个阻尼滑块串联而成.文献[10]首次将Iwan模型应用于连接界面接触力学行为的研究，并对模型中弹簧滑块单元的临界滑移力的分布函数和能量耗散的关系进行了分析.文献[11]详细推导了滑块临界滑移力均匀随机分布条件下，Iwan模型的恢复力和能量耗散的表达式.Song等^[12]在原始Iwan模型的基础上，补充了宏观滑移时的残余刚度，进一步完善了该模型.Segalman^[13]改变了原始模型中临界滑移力的分布形式，提出了基于截断幂级数分布的四参数非均匀密度函数，建立了基于该密度函数的四参数Iwan模型.最近，文献[14]在此基础上，提出了基于截断幂律分布和双脉冲函数的六参数非均匀密度函数，并推导了模型恢复力和能量耗散的解析表达式.

上述模型在建模过程中通常假设法向接触载荷为常数，不随时间变化.在振动环境下，连接界面的法向载荷一般是随时间变化的，而临界滑移力的大小与法向载荷的大小直接相关.如果在动力学分析中，不考虑法向载荷的变化及其与切向载荷的动态耦合，将不能完整地模拟界面接触摩擦行为，导致结构动态响应预测结果误差较大.

本文采用线性弹簧模拟法向载荷的变化，将其与Iwan模型结合起来，建立了界面法向和切向恢复力变化的动态耦合关系，提出了一种新的考虑法向载荷变化的Iwan模型，采用数值方法研究了简化的摩擦阻尼器系统的振动行为.针对不同激励环境，从能量耗散和频响特性的角度对含接触摩擦模型的振动系统进行动力学分析，得到了一些有意义的结论.

1 Iwan模型及其恢复力

Jenkins单元中阻尼滑块与接触面间的切向作用服从库伦摩擦定律.当弹簧受力小于阻尼滑块的临界滑移力时，Jenkins单元处于黏滞状态；相反，则Jenkins单元处于滑移状态.将Jenkins单元并联起来即可得到原始Iwan模型，如图 1所示.Iwan模型的参数包括：弹簧总刚度k_t，并假设每个Jenkins单元的弹簧刚度都相等，即为k_t/n；阻尼滑块的整体临界滑移力f_y，为了模拟连接界面的跨尺度滑动，模型假设每个阻尼滑块的临界滑移力都不同，且满足均匀分布，则每个Jenkins单元上阻尼滑块的临界滑移力为f_i^*/n.

图 2为Iwan模型在周期性激励作用下得到的两种典型的迟滞回线.第1种是a-b-c-d-a所围的面积较小的滞环，此时作用于模型上外激励较小，只使部分阻尼滑块发生滑动，模型在卸载段a-b-c和加载段c-d-a都处于黏滑状态；第2种是a′-b′-c′-d′-a′所围的面积较大的滞环，这种情形下，随着外激励的增大，模型中全部的阻尼滑块都发生滑动，从而出现宏观滑移现象，其中在卸载的前段a′-b′，模型处于黏滑状态，而在卸载的后段b′-c′，模型处于宏滑状态.

骨干曲线oa指首次加载所得的恢复力位移曲线，其表达式为

${T_{oa}}\left( x \right) = \int_0^{{k_t}x} {{f^*}\varphi \left( {{f^*}} \right){\rm{d}}{f^*}} + {k_t}x\int_{{k_t}x}^\infty {\varphi \left( {{f^*}} \right){\rm{d}}{f^*}} ,$

(1)

式中φ(f^*)为阻尼滑块临界滑移力的分布函数.式(1)中第1项表示模型在滑移状态时阻尼滑块的受力，此时弹簧受力大于阻尼滑块的临界滑移力，阻尼滑块处于滑动状态；第2项表示模型在黏滞状态时弹簧的受力，此时弹簧受力小于阻尼滑块的临界滑移力，阻尼滑块处于静摩擦状态.

通过Masing准则来映射骨干函数，从而获取整个周期内的恢复力表达式.Masing准则表述为

$T\left( x \right) = - {\rm{sign}}\left( {\dot x} \right){T_{oa}}\left( {{x_0}} \right) + 2{\rm{sign}}\left( {\dot x} \right){T_{oa}}\left( {\frac{{x - {x_0}}}{2}{\rm{sign}}\left( {\dot x} \right)} \right),$

(2)

式中，x₀为速度转向时的位移，将骨干函数代入式(2)得到模型加载和卸载段的恢复力-位移关系表达式：

$\begin{array}{l} \vec T\left( x \right) = - {T_{oa}}\left( {{x_0}} \right) + 2\int_0^D {{f^*}\varphi \left( {{f^*}} \right){\rm{d}}{f^*}} + 2D\int_D^\infty {\varphi \left( {{f^*}} \right){\rm{d}}{f^*}} ,\\ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftarrow$}} \over T} \left( x \right) = {T_{oa}}\left( {{x_0}} \right) - 2\int_0^{ - D} {{f^*}\varphi \left( {{f^*}} \right){\rm{d}}{f^*}} + 2D\int_{ - D}^\infty {\varphi \left( {{f^*}} \right){\rm{d}}{f^*}} , \end{array}$

式中$D = \frac{{{k_t}\left( {x - {x_0}} \right)}}{2}$.

2 考虑法向载荷变化的Iwan模型

Iwan模型中阻尼滑块的临界滑移力与法向正压力满足库伦摩擦定律，其建立了接触面法向与切向的联系.本文的建模思想主要突出切法向的动态耦合，不重点考虑法向的建模方法，因此法向模型采用简单的线性弹簧来模拟法向正压力的变化.

2.1 恢复力-位移关系推导

假定接触面法向的预加常值载荷为N_ref，以比例参数的形式给出法向载荷的变化关系为

$\alpha \left( t \right) = \frac{{N\left( t \right)}}{{{N_{{\rm{ref}}}}}}.$

Iwan模型中阻尼滑块的临界滑移力满足均匀分布规律φ(f^*/α)=1/2f_y，且其与法向载荷满足库伦摩擦定律f_y=μN_ref，因此分布函数的变化与法向载荷的变化关系相反.

$\varphi \left( {{f^*},\alpha } \right) = \frac{{\varphi \left( {f*/\alpha } \right)}}{\alpha }.$

考虑法向模型参数α对切向Iwan模型参数的影响，进而得到模型恢复力-位移关系的骨干函数为

${T_{oa}}\left( x \right) = \int_0^{{k_t}x} {{f^*}\frac{{\varphi \left( {f*/\alpha } \right)}}{\alpha }{\rm{d}}{f^*}} + {k_t}x\int_{{k_t}x}^\infty {\frac{{\varphi \left( {f*/\alpha } \right)}}{\alpha }{\rm{d}}{f^*}} .$

(3)

将Jenkins单元中阻尼滑块的临界滑移力分布函数代入式(3)可得

${T_{oa}}\left( x \right) = \left\{ {\begin{array}{*{20}{c}} {{k_t}x - \frac{{{{\left( {{k_t}x} \right)}^2}}}{{4\alpha {f_y}}},}&{0 ＜ x \le \frac{{2\alpha {f_y}}}{{{k_t}}};}\\ {\alpha {f_y},}&{x \ge \frac{{2\alpha {f_y}}}{{{k_t}}}.} \end{array}} \right.$

以下各式中下标()₀为速度转向时的参数.当法向载荷变化时，之前使用的Masing准则中的倍数关系将不再适用，需对其进行调整.令x_r表示Jenkins单元从阻尼滑块临界滑移状态卸载至弹簧自由状态过程中的弹簧伸长量，x_t表示Jenkins单元从弹簧自由状态反向加载至阻尼滑块临界滑移状态过程中的弹簧压缩量，则存在：

$x - {x_0} = {x_r} + {x_t}.$

当法向载荷变化时，${x_t} = \frac{\alpha }{{{\alpha _0}}}{x_r}$，则${x_t} = \frac{{\left( {x - {x_0}} \right)}}{{1 + {\alpha _0}/\alpha }}.$.将式(2)中的2倍关系改为1+α₀/α即可得到调整后的Masing准则.

将骨干函数代入调整后的Masing准则可得到模型加载和卸载的恢复力-位移关系表达式：

$\begin{array}{l} \vec T\left( x \right) = - {T_{oa}}\left( {{x_0}} \right) + \left( {1 + \frac{{{\alpha _0}}}{\alpha }} \right) \cdot {T_{oa}}\left( {\frac{{\left( {x - {x_0}} \right)}}{{1 + {\alpha _0}/\alpha }}} \right),\\ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftarrow$}} \over T} \left( x \right) = {T_{oa}}\left( {{x_0}} \right) - \left( {1 + \frac{{{\alpha _0}}}{\alpha }} \right) \cdot {T_{oa}}\left( {\frac{{\left( {{x_0} - x} \right)}}{{1 + {\alpha _0}/\alpha }}} \right). \end{array}$

对恢复力-位移关系表达式进行正则化处理.引入量纲一的量参数：${G_{oa}}\left( x \right) = \frac{{{T_{oa}}\left( x \right)}}{{{f_y}}},G\left( x \right) = \frac{{T\left( x \right)}}{{{f_y}}},q = \frac{{{k_t}x}}{{{f_y}}},{q_0} = \frac{{{k_t}{x_0}}}{{{f_y}}}$.归一化后的模型骨干函数为

${G_{oa}}\left( q \right) = \left\{ {\begin{array}{*{20}{c}} {q - \frac{{{q^2}}}{{4\alpha }},}&{0 ＜ x \le 2\alpha ;}\\ {\alpha ,}&{x \ge 2\alpha .} \end{array}} \right.$

归一化后加载段和卸载段的恢复力-位移关系表达式为：

$\begin{array}{l} \vec G\left( q \right) = - {G_{oa}}\left( {{q_0}} \right) + \left( {1 + {\alpha _0}/\alpha } \right){G_{oa}}\left( {\frac{{q - {q_0}}}{{1 + {\alpha _0}/\alpha }}} \right),\\ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftarrow$}} \over G} \left( q \right) = {G_{oa}}\left( {{q_0}} \right) - \left( {1 + {\alpha _0}/\alpha } \right){G_{oa}}\left( {\frac{{{q_0} - q}}{{1 + {\alpha _0}/\alpha }}} \right). \end{array}$

2.2 模型运动状态的判据

周期性激励作用下，模型的运动状态是不断变化的，而各运动状态间相互转变的判据对于正确地描述模型的特性起关键作用.表 1给出了模型各种运动状态间相互转变时的判据.

在恢复力计算过程中，速度$\dot x\left( t \right) = 0$的判断至关重要.在这一点处，模型的加载方向和运动状态都会发生突变，比如在速度由正向转变为负向过程中，模型的运动状态将可能从t_i-1时刻正向的宏滑状态变为t_i时刻负向的微滑状态.

此外，模型从黏滑状态转为宏滑状态的过程中，宏滑状态的判断也是关键步骤.当模型处于黏滑状态时，其恢复力小于阻尼滑块的临界滑移力，而当模型处于宏滑状态时，其恢复力等于阻尼滑块的临界滑移力.因此，在恢复力计算过程中，需要判断恢复力的大小，并将其限制为小于等于临界滑移力，如果不进行此操作，那么将得不到正确的恢复力-位移关系.

当法向载荷小于零时，连接界面发生分离，恢复力归零.需要注意的是，分离状态时，计算程序中状态转换点的判据发生变化，由速度为零的点变为分离状态的末端点.这与连接界面相接触时有所不同.

2.3 准静态过程

本文给出准静态情况下考虑法向载荷变化的Iwan模型预测的滞回曲线，并将不同法向振荡载荷幅值下对应的结果进行对比.计算参数为：切向弹簧总刚度k_t=100 kN/m，法向弹簧刚度k_n=100 kN/m，接触面滑动摩擦系数μ=0.1，给定切向位移x=0.5sin(2πt)，mm, 法向载荷由预加常值载荷N₀和振荡载荷N₁sin(2πt+δ)两部分叠加而成，其中N₀=100 N.图 3、4分别给出了切法向激励在同相(δ=0°)和异相(δ=90°)情况下，不同法向振荡载荷幅值下对应的模型迟滞曲线变化的示意图，图例中参数η表示法向振荡载荷幅值与预加常值载荷的比值.

当η=0，此时模型退化为原始Iwan模型.原始的Iwan模型不考虑法向载荷的变化，模型参数在加载和卸载过程保持不变，因此其迟滞曲线是对称的；而考虑法向载荷变化的改进Iwan模型则不同，法向载荷的时变性引起了模型的临界滑移力的变化，从而导致宏观滑动的转折点和宏观滑动行为出现较大差异，其迟滞曲线变为非对称的.因此，是否考虑法向载荷的变化对切向摩擦行为的预测有重要影响.

当η＜1，不论法向载荷如何变化，连接界面间总是相互接触的.模型迟滞回线的形状随着法向振荡载荷的变化而变化，其变化规律与切法向的激励形式有关.法向载荷的变化导致模型的宏观滑移力随时间发生变化，从而造成模型各种运动状态的区域发生变化，比如同相激励下，随着η的逐渐增大，迟滞曲线中卸载段表示模型微滑运动状态的区域逐渐增大，表示模型宏滑运动状态的区域逐渐减小；加载段模型微滑运动状态的区域逐渐减小，宏滑运动状态的区域逐渐增大.而异相激励下的变化规律则与之相反.

当η≥1，此时法向载荷会出现小于零的情况，意味着连接界面间会发生间歇性的接触和分离.图 3、4中S所指的区域即为η=1.2时模型发生分离的情形，η越大，滞回曲线中分离状态的区域也越大；且分离状态出现的时刻与切法向激励的相位有关.

3 系统振动分析

摩擦阻尼器常被用来降低复杂连接结构的振动响应，广泛应用于航空发动机领域.本文建立如图 5所示的简化的航空发动机叶片根部摩擦阻尼器模型，利用本文所提出的摩擦接触模型模拟接触面间的力学行为，通过数值仿真的方法研究系统的振动行为.

系统参数设为：质量块m=10 kg，刚度k=3.9×10⁵ N/m，切向刚度k_t=3×10⁵ N/m，法向刚度k_n=3×10⁵ N/m，接触面滑动摩擦系数μ=0.5.系统振动微分方程为

$m\ddot x\left( t \right) + kx\left( t \right) = F\left( t \right) - T\left( {x,\dot x} \right),$

方程两边同时乘以k_t/(mf_y)，得到归一化的方程为

$\ddot q\left( t \right) + \omega _0^2q = \frac{{{k_t}}}{{m{f_y}}}F\left( t \right) - \frac{{{k_t}}}{m}G\left( {q,\dot q} \right),$

式中${\omega _0} = \sqrt {k/m} $.

假定法向载荷和切向载荷的形式为：

$\begin{array}{*{20}{c}} {N\left( t \right) = {N_0} + {N_1}{\rm{sin}}\left( {{\omega _d}t + \delta } \right),}\\ {F\left( t \right) = {F_1}{\rm{sin}}\left( {{\omega _d}t} \right).} \end{array}$

式中：F₁为切向载荷的幅值，ω_d为切法向载荷角频率.

引入量纲一的量参数：$\eta = \frac{{{N_1}}}{{{N_0}}},\beta = \frac{{\mu {N_0}}}{{{F_1}}},\gamma = \frac{{{\omega _d}}}{{{\omega _0}}}$.针对不同的激励条件(体现在参数η、β、γ、δ的变化)，通过大量的数值仿真，研究连接界面各种运动状态下，振动系统的动力学特性.

3个重要的参数：激励频率、切向激励幅值和法向力都会不同程度地影响摩擦阻尼器耗散振动能量的性能^[15].因此本文在一定范围内选取上述3种参数的一些特征离散值来研究振动系统的动力学特性.取N₀=100 N，相位δ=0°，η取4种典型的情形，分别为N₁相对很小的情形η=0.01；适中且连接界面未分离的情形η=0.50；N₁适中且连接界面分离的情形η=1.50；和N₁相对很大的情形η=10.00，而F₁的取值范围则需以激起两个近似线性频率为准.采用数值方法求解摩擦阻尼器系统的动力学方程，对振动系统进行频响分析，结果如图 6所示.图 6中纵坐标表示正则化后的稳态响应.

首先，研究切向激励对系统动态特性的影响.同一法向振荡载荷作用下，振动系统的频响曲线都表现出典型的非线性特征：随着切向激励幅值的增大，振动系统的刚度出现明显的软化现象，其固有频率从高频向低频渐变，响应峰值先降低后增大.

当β=0.1时，切向激励相对较大，此时模型运动以宏观滑动运动形式为主，系统为近似线性系统，其频率比约为1，系统的动响应很大.

当β=1 000时，切向激励相对较小，模型中大多数滑块处于黏滞状态，系统为近似线性系统，其频率比因连接界面分离而不同，当连接界面未发生分离时，此近似线性频率比为1.33，如图 6(a)、(b)所示；而当连接界面分离时，该值将会减小，且法向振荡载荷越小，频率比则越小, 如图 6(c)、(d)所示.此时，由于大多数滑块都未滑动，导致系统的阻尼性能较差，动响应幅值很高.这两种情况下，接触界面间要么滑移状态占主导，要么黏滞状态占主导，其滞回曲线为一扁平状的环线，近似直线段，因此导致系统耗散的振动能量很小，动响应幅值则很高.

当法向相对正压力适中时，振动系统的迟滞回线比较规则、“饱满”，不存在某一运动状态占绝对主导的情形.相对“饱满”的迟滞曲线所围成的面积也就相对越大，说明其耗散的振动能量越多.

再研究法向载荷对系统动态特性的影响.图 6(a)、(b)分别是连接界面相接触时的结果，两者较为接近，主要的差异出现在切向激励幅值适中的情形下，同一激励对应的固有频率和响应都会有微小的差异.

当η=1.50时，连接界面发生分离，系统的非线性程度明显增强，频响曲线越加扭曲.当法向振荡载荷适中时，频响曲线出现了“跳跃”，动响应幅值随激励频率的增加出现突然的增大，此时频响曲线可能出现共振峰向低频段偏转的情形，这种现象在实验或者解析解中能够观测到，但是数值计算过程不能复现该现象^[16]，这就意味着图中的共振峰对应的频率可能略大于其实际峰值频率；此外，当切向激励较小时，频响曲线出现了两个共振峰，且峰值频率比连接界面相接触时小，造成这种现象的原因是界面分离导致的模型运动状态的突变，如图 7所示.将切向小载荷下的滞回曲线近似为直线段，其中S表示界面分离状态，箭头表示滞回曲线绘制的路径.图 7(a)中模型正向加载至A后，模型的运动状态发生突变，恢复力变为零，连接界面发生分离；而图 7(b)则是模型分离但未发生运动状态突变时的示意图.由于模型运动状态的突变，使得模型的非线性程度增强，从而导致振动系统的频响特性复杂化.当η=10.00时，模型的法向振荡进一步增强.随着η的增加，频响共振峰值频率逐渐减小.

在实际的工程结构设计过程中，设计工作者们更加关心的是处于两个近似线性频率之间的频段，在这一频段内，接触界面间的滑动摩擦特性会耗散更多的振动能量，从而为复杂连接结构提供阻尼，降低其动响应幅值.

连接界面的法向分离增强了连接结构的非线性程度，造成系统的阻尼性能降低，使得其稳态动响应幅值增大，瞬态过程延长，且瞬态动响应幅值也增大；在某些情形下，可能会发生非常复杂的非线性现象.因此，实际工程中应尽量避免出现法向分离状态或法向振荡载荷过大的情况.

4 相位的影响

给定切向激励幅值和激励频率，改变法向激励的相位来分别求解连接界面在接触和分离状态时振动系统的滞回曲线.图 8(a)给出了连接界面接触时的滞回曲线，结果显示，随着法向激励相位的增加，模型滞环中表示卸载段黏滑状态和加载段宏滑状态的曲线都逐渐变长，而表示卸载段宏滑状态和加载段黏滑状态的曲线则逐渐变短.图 8(b)给出了连接界面分离时的滞回曲线，图中S表示模型分离状态的区域.分离状态的区域随着法向激励相位的变化而发生变化，有的发生在加载段，有的发生在卸载段，有的则跨越了两个阶段，而且不同激励相位对应的滞回曲线中的分离状态区域的大小也不相同.

当δ=90°时，滞环所围的面积达到最大，其中表示模型运动状态的曲线长度较为均匀；而当δ=0°和δ=180°时，滞环的面积则最小.法向激励相位通过改变法向力的时域分布来改变模型中阻尼滑块临界滑移力的分布，从而造成滞环中表示运动状态的曲线长度和滞环面积发生变化，随即导致振动系统的阻尼性能变化.

5 结论

1) 本文的模型相比于原始Iwan模型能够更加完善、准确地模拟接触界面间的切法向耦合的力学行为.

2) 考虑法向载荷变化与否对连接界面接触力学行为的特性有重要影响.考虑法向载荷变化时，切向恢复力与位移的迟滞回线表现出非对称的特征.

3) 通过对简化摩擦阻尼器系统的动力学分析，发现法向载荷幅值变化比对系统的能量耗散特性有重要影响.连接界面的法向分离增强了系统的非线性程度，造成系统的阻尼性能降低，使得其稳态动响应幅值增大，甚至可能会发生非常复杂的非线性现象.