基于Iwan 模型的接合面切向响应建模1)

占旺龙 李 卫 黄 平

*(深圳技术大学中德智能制造学院,深圳 518118)

†(华南理工大学机械与汽车工程学院,广州 510640)

引言

机械工程结构中广泛存在着以螺纹连接、法兰盘和铆接为代表的搭接接头形式,如图1 所示.搭接连接的典型特征是接合面在法向预紧力作用下发生接触,由于表面粗糙度的存在使得接触界面上承受着非均匀的压力分布.当连接件受到切向力作用时,接触面间的摩擦力会阻碍界面间的宏观相对运动.假设接触界面上摩擦系数处处相等,则显然接触界面上各处抵抗切向相对运动的能力各不相同,法向接触压力大的区域抵抗切向相对运动的能力强,而法向接触压力小的区域则抵抗切向相对运动的能力弱.由于界面接触压力大小与粗糙峰高度正相关,因此在切向力作用下微凸体高度小的先发生滑移,而高度大的粗糙峰后发生滑移,当切向力增大到一定程度时,滑移区将覆盖整个接触界面从而导致连接失效.在往复振荡切向载荷作用下,接触界面上的黏着区与滑移区会交替变化,力与位移关系在笛卡尔坐标系中形成迟滞曲线,迟滞曲线所围成的面积即为振荡加载过程中的能量耗散[1-2].

图1 螺栓搭接接头Fig.1 Bolted lap joint

在传统的机械动力学研究中,一般都采用等效的弹簧阻尼元件来模拟连接界面间的接触力学行为而忽略接合面的非线性特征[3],模型中的参数需要通过实验进行标定.然而,由于连接界面间的非线性特征,在某种条件下标定得到的模型参数并不一定适用于其他实验条件[4].有限元方法也被广泛应用于简单装配结构的建模分析工作中,但是若想要准确地了解接合面的滑移行为则需要更加细密的网格,这将大大增大计算时间.与此同时,对于由随机粗糙面组成的接合面,要建立精确的有限元模型还存在很大的困难,因此开展搭接接头非线性特性研究并建立基于实际可测物理参数的理论模型具有重要的意义[5-7].

目前用于描述接合面非线性力学行为的代表模型有:Iwan 模型[8]、LuGre 模型[9]和Valanis 模型[10]等.但是这些模型大部分都是唯象的,模型中的参数靠实验标定而往往缺少明确的物理含义,因此也就无法反应出接合面粗糙程度对宏观力学行为的影响.Iwan 模型[11]最早用来描述土木结构的弹塑性迟滞行为,而后Segalman[8]通过定义连接界面临界滑移力密度分布函数,将该模型推广至连接界面黏滑行为建模.在小幅振荡周期性载荷作用下,单个加载周期内的能量耗散值D与切应力幅值F0之间的关系可以表示为

式中,ν 和χ 通过实验数据拟合得到,χ 范围在-1 到0 之间[12].利用该关系,Segalman 定义的临界滑移力分布函数为

式中,H(·)为Heaviside 单位阶跃函数;φmax为发生宏观滑移时的位移;S为界面开始发生宏观滑移时力-位移曲线斜率不连续程度,R和S的具体计算方法见文献[8].

从式(2)可以看出,当φ 增大时,ρ(φ)会逐渐减小;当φ=φmax时,ρ(φ)会发生突变;而当φ ＞φmax时,ρ(φ)=0,其对应的图像形式如图2(a)所示,其为含截断幂律分布和单脉冲的四参数非均匀密度函数[8].Song 等[13]和张相盟等[14]基于均匀密度函数建立起摩擦接头Iwan 模型(如图2(b)所示).其他一些比较常见的临界滑移力分布函数如图2(c)和图2(d)所示,具体解释可参阅文献[15].而后Argatov 等[16]、王东等[17]利用Kragelsky-Demkin 粗糙面接触理论求解了螺栓连接接合面间的非线性力学行为.Brake[15]还将Iwan 模型扩展到连接界面发生宏观滑移后的阶段.另外李一堃等[18-19]提出一种含截断幂律分布和双脉冲的非均匀密度函数的六参数Iwan 模型.

图2 临界滑移力密度分布函数图[20]Fig.2 Density distribution function of critical slip force[20]

本文针对前人研究中临界滑移力分布函数需要进行参数标定的缺陷,从接触面统计接触模型出发,推导出新的临界滑移力分布函数,在此基础上对典型工程摩擦副表面进行数值计算并对结果进行分析,最后还将推导出的函数用于与已发表的实验结果进行对比,结果表明理论推导出的分布函数与实验结果具有很好的一致性,该研究成果可为搭接接头的切向响应研究提供参考.

1 Iwan 模型

组成Iwan 模型的两种基本单元如图3 所示,每个基本单元由一个线性弹簧元件和一个库仑摩阻片按照串联或者并联方式组成.并联单元允许在位移不发生改变的情况下改变力的大小;而串联模型允许位移发生改变而力不发生变化.

一般可以将随机粗糙面的接触问题简化成刚性平面与球体高度随机分布的弹性体之间的接触,如图4(a)所示.由于球体高度随机分布,那么每个球体所受到的法向力则会是球体高度的函数.根据库仑摩擦定律,即摩擦副间的摩擦力与法向力成正比,这样由于每个微凸体所承受的法向载荷不同而使得其所能承受的最大切应力也各不相同.与Iwan 模型对比发现[16],可以将摩擦面间的切向接触问题用串-并联单元来描述,如图4(b)所示.该模型由若干个Jenkins 单元并联组成,而每个Jenkins 单元又是由弹性常数为k的弹簧和阈值为qi的库仑摩阻片串联而成.通过类比可以将图4(a)中每个微凸体的接触行为用图4(b)中的一个Jenkins 单元进行描述.

图3 Iwan 模型基本单元Fig.3 Basic element of Iwan model

图4 (a)粗糙面接触模型;(b)简化并-串联Iwan 模型Fig.4 (a)Rough surface contact model;(b)Simplified parallel-series Iwan model

Iwan 模型的工作原理如下:当单个Jenkins 单元所受到的切向力Fi小于摩阻片滑动阈值qi时,弹簧将会发生弹性形变从而承受切向力作用;而当切向力Fi大于摩阻片滑动阈值qi时,摩阻片将会发生滑动,此时该单元所能承受的最大切向力为qi.假设在力F作用下,上表面相对于下表面切向相对位移为x,根据上述描述的Jenkins 切向滑动行为,可以用数学关系表示为

由于微凸体高度随机分布,根据类比可以假定摩阻片滑动力阈值也随机分布.假设其摩阻片临界滑移力密度分布函数为ρ(q),这样当＞0 时,整个体系切向力与相对位移之间的关系可以表示为

式中,ρ(q)dq为滑动力阈值在(q,q+dq)之间的Jenkins 单元数.式(4)中等号右边第一项为处于滑动状态下Jenkins 单元所产生的总切向力,此时所对应的微凸体发生滑动;第二项为处于弹性变形状态下Jenkins 单元所产生的总切向力,此时微凸体处于黏着状态.若式(4)等号右端第二项为零,此时界面发生宏观滑移;而若该项不为零,此时接触界面部分微凸体处于滑移状态.因此摩擦副界面间的运动可以分为两个阶段:

(1)微滑移阶段:当外界施加的切向载荷较小时,接触面承受法向载荷较小的局部会发生滑动,而其他部分处于黏着状态;

(2)宏观滑移阶段:当外界施加的切向载荷大到一定值的时候,接触面间会发生明显的宏观滑动,此时整个接触面都发生相对滑动.

利用以下变量关系式对式(4)进行变换

式中,σ 为微凸体高度分布方均根.经过上述变量代换后,式(4)中变量的量纲将发生变化,由于σ 量纲为长度,因此和φ 为无量纲量,同理的量纲为力.代入式(4)后可得

对上式求一阶导数,可以得到摩擦副间的切向刚度kt为

对式(6)求二阶导数可得

对于一般的摩擦副而言,其上下表面间能承受的最大切向力为摩擦力,即Fmax=µW.从上式可以看出变量φ 的期望为滑动摩擦力,并且当φ →∞时,有若将式(9)两端同时除以摩擦力可得

从上述分析可知,若要将Iwan 模型成功地应用于粗糙表面之间的切向接触问题研究,首先要解决的问题应该是推导出正确的分布函数该分布函数要满足恒大于或等于零及式(7)～式(9)的条件.

2 周期性加载力-位移曲线

将求得的临界滑移力密度分布函数代入式(6)便可得到切向力与切向位移之间的关系.一般来说,分布函数为分段函数,因此在计算加载曲线时要考虑φ 成立的区间范围,例如,当时,切向力计算式为

对于其他区间范围也同样可以用上述方法求出.由于式(11)形式较为复杂,特别是当表面粗糙峰处于非高斯分布时,为得到切向力与切向位移之间的关系,可以用Simpson 数值积分方法进行求解.

当摩擦副承受位移幅值为a的周期性位移激励时,在位移-力图上可以得到如图5 所示的迟滞曲线.其中脊线函数可以根据式(6)得出,而正、反向加载曲线可以通过Masing 模型[21]求出.当体系由图5 所示的F(·)的(a,F(a))位置反向加载到(-a,F(-a))位置时,反向加载力-位移曲线方程为

图5 力-位移迟滞曲线Fig.5 Force-displacement hysteresis curve

将式(12)和式(13)代入式(14),即可得到不同位移幅值下单个加载周期内的能量耗散值.

3 临界滑移力分布函数推导

如图4(a)所示,粗糙面之间的接触问题可以等价成刚性平面与覆盖球体但高度随机分布的弹性面之间的接触.一般用轮廓仪测量得到的是σ,但在现有的统计学模型中用到的是σa,两者存在如下关系[22]

其中,β=σRη,R为微凸体曲率半径,η 为微凸体分布密度,对一般工程表面来说,β 范围为0.02～0.06.当微凸体高度服从正态分布时,分布函数可以写成

对z进行无量纲化,令z*=z/σ,则分布函数可以写成

塑性指数ψ 为[23]

式中,E为综合弹性模量,,E1,2为摩擦副材料的弹性模量,ν1,2为摩擦副材料的泊松比,H为较软材料硬度值,K最大应力接触系数.

根据粗糙表面ZMC 弹塑性接触模型[23],当微凸体法向干涉量ω ＜ω1(ω=z-d,d为摩擦面间的法向接触间距)时,微凸体发生完全弹性变形,此时法向载荷We与微凸体法向形变量ω 的关系为

当ω ＞ω2时,微凸体发生完全塑性变形,此时法向载荷Wp与微凸体法向形变量ω 的关系为

当ω ∈[ω1,ω2]时,微凸体间既有弹性变形,又存在塑性变形,利用样板函数插值方法可到法向载荷Wep与微凸体法向形变量ω 的关系为

式中,ka=1-2K/3,其中K为最大接触应力系数,K=0.454+0.41ν(ν 为材料泊松比).ω1为微凸体开始发生弹塑性变形时的法向临界干涉量,ω1=(πKH/E)2R;ω2为开始发生完全塑性变形时的法向干涉量阈值,ω2=54ω1[24].

利用统计学方法对式(19)～式(21)进行积分可以得到粗糙表面接触时的法向载荷与变形量的关系为

式中,φ(z)为微凸体高度分布函数,An为名义接触面积,具体推导过程可以参见文献[23-24].由于法向载荷W与法向接触距离d之间为单调关系,即接触间距越小,接触面间的法向载荷越大.因此在给定法向载荷的条件下,可以通过二分法计算得到给定条件下的法向接触间距d.

根据库仑摩擦定律,摩擦界面间的摩擦力为

式中所有带星号(*)的变量都是用高度分布均方根σ 进行无量纲化后的变量,µ为滑动摩擦系数.对比式(23)与式(9),由于积分计算与积分变量无关,因此将ω*替换为φ 便可得到分布函数,为

4 实验结果验证

Eriten 等[25]对不同预紧力作用下螺栓搭接接头的切向行为进行了实验研究,组成摩擦副的材料为420 型不锈钢,接触表面经铣削加工而成.材料参数为:弹性模量E=200 GPa,泊松比ν=0.24,硬度H=5.825 GPa,名义接触面积An=156 mm2,表面粗糙度参数σ/R=0.088 8,β=0.023,σ=2.677 μm,ψ=4.725,摩擦系数µ=0.5,具体实验装置及过程可以参考文献[25].需要注意的是,在文献[25]中摩擦系数为一与法向载荷有关的参数,但在本文中设定摩擦系数为发生宏观滑动时摩擦力与法向载荷的比值,认为其为一恒定量,本文取摩擦副发生宏观滑动时摩擦系数为计算值,约为0.5.

图6 给出了不同法向载荷下对摩擦副进行切向循环位移加载得到的迟滞曲线图.其中实线为用本文中模型计算得到的迟滞曲线,点状线的为文献[25]中实验得到的曲线.从图中可以看出理论计算值与实验结果很吻合,从而证明了本文模型的可行性.

图7 给出了不同法向载荷下单个加载周期内的能量耗散实验值与理论计算结果图,其中实验值是统计10 个加载周期内得到的统计箱式图,黑色菱形为计算结果.从图中可以看出计算结果与实验结构很吻合,进一步证明了本文提出模型的有效性.

5 数值计算结果与讨论

下面以钢-钢接触为例,计算中所用到的材料物理参数为[23]:弹性模量E1=E2=207 GPa,泊松比ν1=ν2=0.29,硬度H1=H2=1.96 GPa,名义接触面积An=100 mm2,摩擦系数µ=0.5.由式(24)可知滑移力密度分布函数与形貌参数σ和R都有关,根据McCool[22]的研究,曲率半径R、粗糙峰分布密度η 及方均根σ 与测量的表面形貌高度z(x)之间为耦合关系.因此对于一般的工程应用表面,形貌参数可以用两个参数描述,即:β和σ/R.而这两个参数直接决定塑性指数ψ 的大小,所以下面仅讨论不同塑性指数ψ 对本构关系的影响.表1 给出典型工程表面的形貌参数及对应的塑性指数值.

图6 不同法向载荷及位移幅值下的迟滞曲线实验与数值结果对比Fig.6 Experimental and numerical results comparision of hysteresis curves under different normal loads and displacement amplitudes

图7 单个周期内的能量耗散值实验与数值结果对比(箱图为实验值,黑色菱形为计算值)Fig.7 Experimental and numerical results comparision of energy dissipation per cycle(the box shows the experimental value,and the black diamond is the numerical value)

表1 典型工程表面参数及塑性指数Table 1 Typical engineering surface parameters and plasticity index

假设微凸体高度满足正态分布,根据上述推导公式得出如图8 所示的在不同塑性指数和不同法向载荷下的摩阻片滑动阈值分布函数和切向刚度图.其中图8(a)和图8(b)分别为塑性指数为0.7 和2.5 时的临界滑移力分布函数图.由式(5)中第二式可知,对于给定的粗糙表面,φ 是摩阻片滑动阈值q的度量值,φ 越大即摩阻片滑动阈值q越大.另外由式(10)可知,的数学期望为1,即图8(a)和图8(b)图所示分布曲线的期望值为1.对比图8(a)和图8(b)图中不同法向载荷下的分布函数可以发现,随着法向载荷的减小,出现概率最大的摩阻片滑动阈值也随之减小,这主要是因为法向载荷越小,表面间接触间距d变小,此时单个微凸体承担的法向载荷减小从而使其发生相对滑动所需的切向力也随之减小,即摩阻片滑动阈值减小.对于给定的法向载荷,虽然φ 的取值范围为[0,+∞],但当φ ＞1 时,迅速趋于零,这是因为为了保证期望为恒定值,必须在达到一定程度后迅速减小到零,法向载荷越小,衰减速度越快.从微凸体高度分布函数也可以进行解释,因为虽然正态分布中自变量的取值范围为[-∞,+∞],但在±3σ范围内的微凸体高度所占比例达到99.73%,在该范围以外的可以忽略不计,另外实际工程表面亦不存在无穷高度的微凸体,这些因素都保证了具有水平渐近线.对比相同载荷下不同塑性指数时的分布曲线可以看出,塑性指数越大,φ 的众数越小,众数出现的概率越大.这主要是因为随着塑性指数的增大,高度较高的微凸体将发生塑性变形,该状态下法向干涉量增大但法向载荷保持不变,这样就会导致更多其他的微凸体发生接触但实际承担的法向载荷较小,从而导致随着塑性指数的增大而众数减小,同时出现的概率增大.图8(c)和图8(d)为不同法向载荷下的切向刚度随切向位移的变化关系图,由式(7)可知图中单条曲线为在给定载荷下临界滑移力分布函数的尾部分布积分值.从图中可以看出随着φ 的增大,切向刚度迅速减小到零,此时摩擦副间将发生宏观滑动.φ=0 所对应的量即为初始切向刚度,从图中可以看出随着法向载荷的增大,初始切向刚度也随之增大,这也与文献[26]中的结论相一致.同时塑性指数越大,初始切向接触刚度也越大.

图8 不同法向载荷及塑性指数下的分布函数及切向刚度Fig.8 Distribution function and tangential stiffness under different normal load and plasticity index

图9 为在塑性指数ψ=0.7 时无量纲切向力与切向位移关系图,从图中可以看出加载曲线为凸函数,随着切向位移的增大,切向力也同样随之增大,切向力增大速率逐渐减小并趋于零,表明切向力逐渐趋于滑动摩擦力,此时摩擦副发生宏观滑动.法向载荷越大,发生宏观滑动所需的切向位移越大.

图9 无量纲化切向力-位移曲线(ψ=0.7)Fig.9 Dimensionless tangential force-displacement curve(ψ=0.7)

图10 为当塑性指数ψ=0.7 及法向载荷W=100 N 条件下对摩擦副进行周期性加载时不同位移幅值下用Masing 准则求解的力-位移迟滞曲线,由图中可以看出不同的位移幅值下,脊线曲线都相同.随着位移幅值的增大,迟滞曲线所围成的面积增大,表明单个周期内能量耗散值随位移幅值的增大而增大.

图10 不同切向位移幅值下的迟滞曲线(ψ=0.7,W=100 N)Fig.10 Hysteresis curves under different tangential displacement amplitudes(ψ=0.7,W=100 N)

6 总结

本文对以螺栓连接为代表的搭接接头为例,对在法向预紧力作用下的粗糙面接触问题进行研究.基于Iwan 模型研究界面接触的非线性力学性质,重点推导了Iwan 模型中摩阻片临界滑移力密度分布函数.相比于前人研究,本文推导出的分布函数不再是唯象描述,而是从材料性能参数和表面粗糙度参数推导得出.在求出临界滑移力密度分布函数的基础上,进而求解出接触面间的切向接触响应和在周期性位移加载条件下的迟滞曲线及能量耗散.主要结论如下:

(1)临界滑移力函数开始迅速上升,到达最高点后迅速收敛到零.收敛速度和最大值出现位置与法向载荷、塑性指数及粗糙度参数有关.

(2)界面切向接触刚度是滑移力分布函数的尾部分布.刚开始时切向刚度最大,进一步增大切向位移时,切向刚度逐渐减小到零,此时界面发生宏观滑动.初始切向刚度与法向载荷、粗糙度参数及塑性指数有关.对于确定的接触表面,法向力越大,初始切向刚度越大;初始切向刚度同样也随着塑性指数的增大而增大.

(3)使用本文推导出的函数计算得到的单位加载周期内的迟滞曲线与实验结果很吻合,计算出的单个周期内的能量耗散也与实验结果吻合,证明模型的正确性.