含润滑介质的正交各向异性结合面接触特性研究.pdf

第 36 卷第 4 期2023 年 8 月振 动 工 程 学 报Journal of Vibration EngineeringVol.36 No.4Aug.2023含润滑介质的正交各向异性结合面接触特性研究王世军1，崔圣奇1，吴敬伟2，卫娟娟1，李鹏阳1（1.西安理工大学机械与精密仪器工程学院，陕西 西安 710048；2.西安昆仑工业（集团）有限责任公司技术部，陕西 西安 710043）摘要:提出了一种含润滑介质的正交各向异性结合面法向接触刚度的分形模型。该模型基于含椭圆修正因子的Hertz接触理论，并根据考虑润滑介质的侧向泄漏的平均雷诺方程，推导出固体接触刚度与流体刚度之间的解析关系。通过分析不同因素对结合面的法向接触刚度的影响，发现当无量纲固体真实接触面积小于 0.05时，结合面的法向接触刚度受流体刚度的影响较大。随着润滑介质发生侧向泄漏，固体真实接触面积逐渐增大，固体接触刚度对结合面的法向接触刚度的影响越来越显著。给定不同预紧力，对比试验与有限元仿真获得的前三阶固有频率，其最大相对误差为4.11%，证明本文构建的模型可以准确地预测结合面的接触性能。关键词:椭球形微凸体；Hertz接触理论；法向接触刚度；平均雷诺方程；有限元法中图分类号:TH113.1；TB123 文献标志码:A 文章编号:1004-4523（2023）04-0909-14 DOI：10.16385/ki.issn.1004-4523.2023.04.004引言机械结构的有限元仿真分析中，对单个零件进行模态或应力、变形分析时，通常可以得到较为准确的结果，仿真结果与试验结果可以很接近。而整机结构是由多个零件装配而成的，对整机结构进行有限元仿真得到的结果与试验结果往往存在较大差异，主要原因在于结合面的影响。由于含润滑介质的正交各向结合面在整机结构中大量存在，即便是非滑动的固定联接面之间，通常也存在润滑介质，完全纯净的接触表面在实际的机械设备中通常并不存在，所以建立准确的结合面接触模型是提高机械整机性能分析准确性的关键1。不同的表面加工方式会产生不同的表面纹理，导致表面接触特性存在差异。根据表面纹理特征的不同，常见的机械加工表面可以分为各向同性结合面和各向异性结合面。文献 26 探讨了不含润滑介质的各向异性结合面的接触特性，文献 79 研究了含润滑介质的各向同性结合面的接触特性。在各向异性结合面的研究中，Chung 等2将微凸体形状假设为椭球体，考虑不同椭圆率对接触点接触变形的影响，构建了粗糙表面椭球形微凸体弹塑性微观接触的分形模型，但该模型并未对接触刚度进行分析。Buczkowski等3应用统计模型对各向异性粗糙表面的接触特性进行研究，认为微凸体峰顶可以用椭圆抛物线表示。王世军等4提出一种结合部法向刚度预估方法，采用接触表面的二维轮廓测量数据建立接触分析模型，预测相互接触的两个各向异性表面在加工纹理相同时的接触刚度。刘伟强等5通过映射法求解出曲面曲率，并依据 Hertz曲面接触理论，提出了椭圆抛物体形微凸体弹性接触的建模方法。谭文兵等6应用各向异性分形几何理论，将微凸体形状等效为抛物面形，对结合面的法向接触刚度进行分析。针对含润滑介质的结合面接触特性，屈重年等7基于分形理论和平均雷诺方程，推导出结合面的固体接触刚度和液体的等效接触刚度，提出了两个粗糙接触表面在含润滑介质状态下的动态刚度建模方法。高志强等8基于 GW，KKE 和 AF 模型，结合流体动力润滑方程，研究了动态载荷对固液结合面接触特性的影响。李玲等9对等效润滑介质厚度与固体刚度之间的关系进行研究，求解了流、固两种刚度的耦合解析式，建立了混合润滑状态下结合面的动态接触刚度模型。本文应用含椭圆修正因子的 Hertz 接触理论和考虑润滑介质的侧向泄漏的平均雷诺方程，研究了含润滑介质的正交各向异性结合面的接触特性，分析了相关参数对结合面的法向接触刚度的影响规律，最后将试验模型和有限元分析模型得到的固有收稿日期:2021-12-02；修订日期:2022-03-01基金项目:国家自然科学基金资助项目（52075439）。振 动 工 程 学 报第 36 卷频率对比，建立了结合面的法向接触刚度模型。1正交各向异性结合面平面磨削的表面，微观纹理表现出正交各向异性特性，两个这样的粗糙表面按照相同纹理方向组成的结合面，也表现出正交各向异性。图 1 为粗糙度Ra=0.8的平面磨削表面。对于正交各向异性结合面，其表面上带有条形纹理，且纹理具有方向性，为了使建立的模型更符合纹理特征，本文将正交各向异性结合面上微凸体的形状等效为椭球形。假设椭圆接触面的长轴与粗糙表面上顺纹理方向相同，忽略椭球形微凸体之间相互作用对接触特性的影响。2含润滑介质的结合面建模图 2（a）为含润滑介质的结合面法向接触示意图。将结合面的弹性接触等效为刚性理想平面与单个粗糙面的接触，如图 2（b）所示，其中 h表示刚性理想平面与单个粗糙表面的法向距离。依据文献10，将结合面的法向接触刚度Kn按照并联弹簧模型等效为固体接触刚度Ka与流体刚度Kl，如图 2（c）所示。结合面的法向接触刚度为：Kn=Ka+Kl（1）3结合面固体刚度建模3.1椭球形微凸体的接触模型图 3 表示两个椭球形微凸体正向接触的几何形状。原点附近的曲面表示为：z1=A1x2+B1y2+C1xy+（2）在此略去了 x 和 y 的更高次项。通过选取适当的 x，y 轴和 x1，y1轴，可消去式（2）中的 xy 项，则式（2）可改写为：z1=12R1x21+12R1y21（3）式 中 R1，R1为 接 触 曲 面z1在 原 点 处 的 主 曲 率半径。与z1接触的曲面方程为11：z2=-12R2x22+12R2y22（4）图 1 平面磨削表面Fig.1 Planar grinding surface图 2 含润滑介质的结合面法向接触示意图Fig.2 Schematic diagram of normal contact of the jointsurface with lubrication medium图 3 椭球形微凸体正向接触的几何形状Fig.3 Geometry of ellipsoidal asperities with positive contact910第 4 期王世军，等：含润滑介质的正交各向异性结合面接触特性研究式中 R2，R2为与接触曲面z1对称的接触曲面z2在原点处的主曲率半径。两椭球形微凸体之间的距离h1=z1-z2。联立式（3）和（4）得：z1-z2=12R1x21+12R1y21+12R2x22-12R2y22（5）将式（5）中含有R，R的项进行合并，如下式所示：z1-z2=Ax2+By2=12Rx2+12Ry2（6）式中 R，R为相对主曲率半径。主曲率坐标轴 x1与 x2的夹角如图 4 所示。主曲率半径与 A，B之间的关系为：A+B=12(1R1+1R1+1R2+1R2)（7）B-A=12(1R1-1R1)2+(1R2-1R2)2+2(1R1-1R1)(1R2-1R2)cos(2)12（8）B-AA+B=cos（9）引入与R，R有关的等效曲率半径12：Re=RR=12AB（10）椭圆接触面的长、短半轴a，b的参数值分别为：a=m3P4E*()A+B 13（11）b=n3P4E*()A+B 13（12）式中 m，n为主曲率半径夹角相关系数；E*为材料的复合弹性模量，且有1E*=1-21E1+1-22E2；E1和E2分别表示材料各自的弹性模量；1和2分别表示材料各自的泊松比。图 5为椭球形微凸体微观接触示意图。可以看出，微凸体承受方向为法向，大小为P的载荷。在P的作用下，微凸体产生的变形量为。3.1.1弹性阶段当 c时，为弹性变形阶段。此时微凸体弹性载荷pe、弹性接触面积se及临界变形量c可分别表示为2：pe=43E*Re1232F2-23(e)（13）se=ReF2-1(e)（14）c=(KH2E*)2ReF-232(e)（15）式中 硬度系数K=0.454+0.4，为两材料中较软材料的泊松比；H为较软材料的硬度，H=2.8y，y为较软材料的屈服强度；F2(e)为椭圆修正因子。3.1.2弹塑性过渡阶段当c 110c时，为弹塑性变形阶段。弹塑性变形又可以被分为两个阶段，当c 6c时，为弹塑性第一阶段；当6c 110c时，为弹塑性第二阶段。此时微凸体弹塑性载荷pep及弹塑性接触面积sep可分别表示为11：当1 c 6时：pep1=23KHRec 1.03(c)1.425F-232(e)（16）sep1=0.93Rec(c)1.136（17）当6 c 110时：pep2=23KHRec 1.40(c)1.263F-232(e)（18）sep2=0.94Rec(c)1.146（19）图 4 主曲率坐标轴夹角Fig.4 The angle between the coordinate axes of principal curvature图 5 椭球形微凸体微观接触示意图Fig.5 Schematic diagram of microscopic contact between ellipsoidal asperities911振 动 工 程 学 报第 36 卷3.1.3完全塑性阶段当 110c时，为塑性变形阶段。此时微凸体塑性载荷pp及塑性接触面积sp分别为11：pp=2HReF-232(e)（20）sp=2Re（21）图6为椭圆修正因子F2(e)与椭圆形状b/a，相对主曲率半径比值R/R的关系曲线。图中，F1，F2，F3为椭圆修正因子。椭圆修正因子F2(e)的表达式为：F2(e)=2()ba124e2()ba3213()ab2K2(e)-K1(e)16 K1(e)-K2(e)164e2()ba32K1(e)（22）式中 椭圆离心率e=1-b2/a2；第一类椭圆积分K1(e)=02d1-e2sin2；第 二 类 椭 圆 积 分K2(e)=021-e2sin2 d。3.2椭球形微凸体的法向接触刚度模型通过 WM 函数模拟粗糙表面轮廓，得到微凸体的变形量为1314：=23-DD-22GD-1(ln)12s2-D2（23）式中 D 为分形维数，1 D sds|=D22-D2slD2s-()D2+1（34）式中 s为微凸体接触截面积，s=2s，s为微凸体接触面积；sl为最大接触点截面积，sl=2sl，sl为最大接触点面积；为域扩展因子。D与之间的解析关系式为16：2-DD-(1+-D2)-2-DD=2-DD，1（35）结合面真实接触面积为：图 6 椭圆形状b/a、椭圆修正因子F2(e)与相对主曲率半径的比值R/RFig.6 Ratio of elliptical shape b/a，elliptical correction factor F2(e)to relative principal curvature of radius curvature R/R912第 4 期王世军，等：含润滑介质的正交各向异性结合面接触特性研究Ar=scsln(s)sds+()161D-1scscn(s)sds+()11101D-1sc()161D-1scn(s)sds+0()11101D-1scn(s)sds（36）式中 sc为微凸体临界接触截面积，sc=2sc，sc为微凸体临界接触面积。将式（34）代入式（36）得：Ar=2D2-D2-D2sl（37）结合面法向接触总载荷为17：P=Pe+Pep1+Pep2+Pp=scslpen(s)ds+()161D-1scscpep1n(s)ds+()11101D-1sc()161D-1scpep2n(s)ds+0()11101D-1scppn(s)ds（38）当1 D 2且D 1.5，结合面无量纲法向接触载荷为时：P=212-3D2G*()D-1()ln 12()2-D22D()2-DDD233-D2()3-2DAD2r()2-DD3-2D2-2D2+7D-64A3-2D2r-s3-2D2cF-162()e+24.25-1.275DG*0.85()D-1()ln 0.425()2-D22D()2-DDD23()1.425-0.925D1.03AD2r()2.8K0.150.425D-1.275s1.425-0.925Dc1-()161.425-0.925DD-1F231202()e+22.63-0.789DG*0.526()D-1()ln 0.263()2-D22D()2-DDD23()1.263-0.763D1.40AD2r()2.8K0.4740.263D-0.789s1.263-0.763Dc()161.263-0.763DD-1-()11101.263-0.763DD-1F1917502()e+2K()2-D22D()2-DD2-D2AD2rs2-D2c()11102-D2D-22-DF132()e（39）式中 Ar为结合面无量纲真实接触面积，sc为无量纲临界接触面积，G为无量纲分形粗糙度参数。当1 D 2且D=1.5时，结合面无量纲法向接触载荷为：P=23.75G*0.5()ln 0.50.0625 1.530.75 (Ar3)0.75ln(Ar30.1875sc)F-162(e)+1.03 22.34G*0.425()ln 0.4250.0625 1.50.1125()Ar30.75 ()2.8K0.15-0.6375s0.0375c 1-()160.075F231202()e+1.40 21.45G*0.263()ln 0.2630.0625 1.50.3555()Ar30.75 ()2.8K0.474-0.3945s0.1185c()160.237-()11100.237 F1917502()e+2K0.0625 1.5()Ar30.75s0.25c()11100.50.5F132()e（40）式中 为材料的塑性指数，=y/E*。固体法向接触刚度Ka为18：Ka=Ke+Kep1+Kep2=scslken(s)ds+()161D-1scsckep1n(s)ds+()11101D-1sc()161D-1sckep2n(s)ds（41）将式（28），（30），（32），（34）代入式（41）中，得结合面法向接触刚度Ka：Ka=4E*2-D2D()3-D32()1-D()2-D(s12l-sD2ls1-D2c)F-762(e)+1.03 20.825+0.15DE*0.852-D2(KH)0.1530.075D+0.275G0.15(D-1)(ln )0.075(2-D)D(1.425-0.425D)0.425(1-D)sD2ls0.425(1-D)c 1-(16)-0.425 F-256130002(e)+1.40 20.474D-0.633E*0.5262-D2(KH)0.47430.237D-0.211G0.474(D-1)(ln )0.237(2-D)D(1.263-0.263D)0.263(1-D)sD2ls0.263(1-D)c (16)-0.263-(1110)-0.263 F-3975002(e)（42）913振 动 工 程 学 报第 36 卷将式（42）量纲一化可得固体无量纲法向接触刚度：K*a=4()2-D22D()3-D32()2-D()1-D(2-DD)D2A*D2r ()2-DD1-D2-D2+3D-24A*r1-D2-s*1-D2cF-762(e)+1.0320.825+0.15D()2-D22()2.8K0.15D30.075D+0.275G0.15()D-1()ln 0.075()2-D ()1.425-0.425D0.425()1-D(2-DD)D2AD2rs0.425()1-Dc 1-(16)-0.425 F-256130002(e)+1.40 20.474D-0.633()2-D22()2.8K0.474D30.237D-0.211G0.474()D-1()ln 0.237()2-D ()1.263-0.263D0.263()1-D(2-DD)D2AD2rs0.263()1-Dc (16)-0.263-(1110)-0.263 F-3975002(e)（43）4结合面流体刚度建模4.1流体刚度建模的基本假设在流体刚度的建模中，假设润滑介质是牛顿流体，忽略润滑介质厚度方向压力的变化。由于润滑介质厚度仅百分之几毫米，压力不会发生明显变化，接触表面的润滑介质的流速与表面速度相同，沿润滑介质厚度方向的黏度为定值。4.2一般形式的雷诺方程一般形式的雷诺方程为19：x(h32fx)+y(h32fy)=6 x(Uh2)+y(Vh2)+2h2t（44）式中 U=U0+Uh为粗糙表面在 x 方向的速度之和；V=V0+Vh为粗糙表面在 y方向的速度之和；f为润滑介质挤压膜承受的压力；h2为润滑介质的厚度；为润滑介质的黏度；为润滑介质的密度；t 为时间。4.3两平行板之间的载荷模型图 7为挤压润滑介质压力分布示意图。假设两平行板在 y方向的宽度为 L且为无限长，x方向的板长为 B；且在 y方向无侧向泄漏。现将两平行板无限靠近，考虑到dh2dt为负，可将雷诺方程改写为：ddx(h3dfdx)=-12dh2dt（45）由于两板平行，所以h2不是 x的函数，则可将式（45）改写为：d2fdx2=-121h32dh2dt（46）对式（46）积分得：f=-61h32dh2dtx2+Cx+C1（47）根据边界条件：当x=0，f=fmax，dpdx=0，则 有 C=0；当x=B2，f=0，则C1=614h32dh2dtB2。将 C，C1代入式（47）中，得：f=61h32dh2dt(B24-x2)（48）则两板之间的载荷为：W=-B2B2fLdx=6L1h32dh2dt-B2B2(B24-x2)dx=L1h32dh2dtB3（49）4.4接触面积与流量因子的关系将雷诺方程改写为含有流量因子的形式20：图 7 挤压润滑介质压力分布示意图Fig.7 Schematic diagram of squeeze lubrication medium pressure distribution914第 4 期王世军，等：含润滑介质的正交各向异性结合面接触特性研究x(xfx)+y(yfy)=c12h32h2t（50）式中 x，y为压力流量因子；c为接触流量因子。圆形接触面是各向同性结合面上两球形微凸体受到载荷作用形成的。将接触面形状由椭圆形等效为圆形。假设接触面内长、短轴a，b相等且等于等效接触半径 c，即a=b=c。正交各向异性结合面在 x，y方向的压力流量因子可以等效为各向同性结合面的压力流量因子。正交各向异性的结合面的流量因子与名义接触面积Aa、固体真实接触面积Ar的关系为：x=y=ArAa（51）接触因子c表示润滑介质所占的比例，则有：c=A0Aa=Aa-ArAa（52）4.5不考虑 y方向的侧向泄漏当泄漏仅存在于 x方向时，式（45）可以表示为：d2fdx2=-cx121h32dh2dt（53）则两平行板之间的载荷为：Wf=h32cxdh2dtLB3（54）4.6考虑 y方向的侧向泄漏图 7 所示两平行板在 y 方向上的长度均为有限长。由于受到挤压，在 x 和 y 方向都会发生侧向泄漏。当考虑 y 方向的侧向泄漏时，需要在式（54）中引入由 Cameron21根据B L获得的载荷修正系数。表 1给出了B L与的对应值。无侧向泄漏时，B L=0。考虑侧向泄漏情况下的载荷可以表示为：Wf=h32cxdh2dtLB3（55）对式（55）进行修正并无量纲化：Wf=WfAaEl=AaElcxh32dh2dtLB3（56）式 中 El为 润 滑 介 质 的 弹 性 模 量，且 有1El=1(1-21E1-1-22E2)。对式（56）求导，则流体法向刚度为22：K*l=dW*fdh*2=3B2Elh4*2dh*2dt(A-1r-1)（57）式中 h*2为润滑介质的无量纲厚度，h*2=hh0。结合面无量纲法向接触总刚度为：K*n=K*a+K*l（58）5分析结果与讨论图 8为不同分形参数D下A*r与K*n的关系曲线，分 别 选 取A*r=00.5，=0.018 Pa s，G*=1 10-9，=0.1，D=1.21.9。图 9 为分形维数D=1.4时G*与K*n的 关 系 曲 线，分 别 选 取G*为1 10-9，1 10-10，1 10-11。图 10 为分形维数D=1.6时与K*n的 关 系 曲 线，分 别 选 取为0.018，0.15和0.4 Pa s。图 11 为 分 形 维 数D=1.4时与K*n的关系曲线，分别选取为0.01，0.1，1。图 12 为分形维数D=1.6时A*r与K*的关系曲线，其中K*包括结合面的无量纲法向接触总刚度K*n、固体接触刚度K*a与流体刚度K*l。图 13为润滑介质侧向泄漏对K*n的影响曲线。图 14为结合面纹理方向对K*n的影响曲线。对图814进行分析可得如下结果：（a）根据图 8 可知，当D 1.4时，随着A*r的增大K*n呈现先减小后增大的趋势。由图 8 还可以看出，当D 1.7时，K*n与A*r呈近似线性关系。（b）根据图 9 可知，随着G*的增大K*n反而减小，这是因为分形粗糙度越大，结合面越粗糙，摩擦系数越大，导致结合面的法向接触刚度降低。（c）根据图 10 可知，随着的增大K*n也在逐渐增大。（d）根据图 11可知，当增大，K*n也随之增大，材料的屈服强度增大，无量纲临界接触面积呈现减小的趋势。粗糙表面间弹性接触面积的增加，导致表 1 载荷修正系数Tab.1 Correction coefficient of loadsB L15/62/31/22/51/31/41/51/1000.4210.4980.5800.6330.7480.7900.8450.8740.9371915振 动 工 程 学 报第 36 卷结合面的法向接触刚度增大。（e）图 12 为K*n，K*a，K*l与A*r的关系曲线。当A*r 0.05时，K*l对K*n的影响较大；随着A*r逐渐增大，K*a逐渐成为影响K*n的主要因素。（f）根据图 13 可知，与考虑润滑介质侧向泄漏时的K*n相比，不考虑润滑介质侧向泄漏时的K*n较大。这是因为不考虑润滑介质侧向泄露时的法向载荷为式（56），考虑润滑介质侧向泄露时的法向载荷为式（57）。公式（57）在公式（56）的基础上引入取值范围在0.421 1之间的修正系数。但随着真实接触面积的增加，两种情况下的K*n越来越接近。（g）由图 14 可以看出，若考虑结合面加工纹理方向为各向同性，则整个结合面的法向接触载荷和法向接触刚度都将增大，但增大的趋势不明显。这是因为各向同性结合面只是正交各向异性结合面的图 8 D，A*r对K*n的影响Fig.8 Influence of D and A*r on K*n图 9 G*对K*n的影响Fig.9 Influence of G*on K*n图 10 对K*n的影响Fig.10 Influence of on K*n图 11 对K*n的影响Fig.11 Influence of on K*n916第 4 期王世军，等：含润滑介质的正交各向异性结合面接触特性研究一种特殊情况。6试验及有限元验证图 15 中的试验装置包含两个金属板与四个M20 螺栓。上、下试样用螺栓固定，两金属板接触面为具有正交各向异性特性的平面磨削表面。上金属板的材料为钢，下金属板的材料为铸铁。上下金属板长 400 mm，宽 300 mm，厚度分别为 25 mm 和 40 mm。两个金属板表面并未完全接触，在上金属板的接触表面加工凹槽，凹槽尺寸为宽 200 mm，深 5 mm，上、下金属板的材料特性参数如表 2所示。6.1试样表面形貌分析本文选用徕卡共聚焦显微镜 DCM3D 对试样真实表面进行测量。采样间隔在 x，y方向均为 1 m，z 方向均为 0 m，采样区域为1000 1000 m2。图16为试样表面的轮廓图。图 17为上、下试样的表面轮廓曲线。本文使用功率谱密度法（PSD）21得到上、下试样的分形维数 D 分别为 1.855 和 1.903，分形粗糙度G分别为7.837 10-11和2.041 10-10。表 2 金属板的机械性能Tab.2 Mechanical properties of metal planes金属板上金属板下金属板弹性模量/GPa207101泊松比0.280.245密度/（kg m3）81006980屈服强度/MPa353240硬度/MPa500700图 15 试验装置示意图（单位：mm）Fig.15 Schematic diagram of the experimental device（Unit：mm）图 16 三维正交各向异性表面形貌Fig.16 Three-dimensional orthotropic surface profile图 12 A*r与K*的关系Fig.12 Relationship between A*r and K*图 13 润滑介质侧向泄漏对K*n的影响Fig.13 Influence of the lateral leakage of lubricating medium on K*n图 14 结合面纹理方向对K*n的影响Fig.14 Influence of texture direction of joint surface on K*n917振 动 工 程 学 报第 36 卷6.2试验模态分析锤击法模态试验装置如图 18 所示。在金属板表面布置 9个激振点，并在金属板边缘均匀排布 8个加速度传感器拾振，传感器的相关参数如表 3所示。使用阿尔泰 USB2085 动态数据采样仪对数据进行采集和转化。最后对已获取信号进行时域、频域分析来获取模型的固有频率。结合面法向面压选取最大值 1.0 MPa，试验中接触面的名义接触面积为 0.08619 m2，每个螺栓上的拧紧力矩为 86.2 N m。试验振型如图 19 所示，从上到下分别为第一阶振型，第二阶振型，第三阶振型。6.3有限元建模及结果对上述试验装置进行有限元建模分析1，23。本文应用文献 2426 提出的虚拟材料的方法，将结合面刚度模型用以下参数来等效：（1）法向弹性模量：En=KEh3Aa（59）式中 h3表示虚拟接触层等效厚度；K表示结合部量纲一法向接触刚度。（2）切向弹性模量：E=D2-D2alAa()2-DE（60）（3）切变模量：图 17 试样的表面形貌曲线Fig.17 Profile curve of specimen surface图 18 锤击法试验装置Fig.18 Experimental device of hammer pulse method表 3 IC压电式加速度传感器相关参数Tab.3 Parameters of IC piezoelectric acceleration sensors仪器相关参数采样频率/Hz量程/g灵敏度/（mV g1）分辨率/g参数数值0.710000501000.0002图 19 试验振型图Fig.19 Vibration shapes from experiments918第 4 期王世军，等：含润滑介质的正交各向异性结合面接触特性研究G=KtGh3Aa（61）式中 Kt表示结合部无量纲切向接触刚度，G表示等效剪切模量。（4）法向泊松比：n=-xz（62）式 中 x表 示 法 向 应 变，分 析 得x=0，则 式（62）为 0。（5）切向泊松比：=-yx（63）式 中 y表 示 切 向 应 变，分 析 得y=0，则 式（63）为 0。（6）虚拟材料厚度为 1 mm。选择 8 节点的 SOLID185 类型单元，将模型共划分为 11140 个单元，其中虚拟薄层数量为 852 个，如图 20所示。将非线性的有限元模型线性化后，建立包含常刚度接触层的有限元整机线性模型。通过对整机线性模型的模态分析获得包含接触层影响的整机振型。图 21为有限元仿真第一阶到第三阶振型。图 22为有限元仿真计算流程。首先，建立粗糙表面接触特性理论模型。其次，基于虚拟材料法建立包含粗糙表面接触特性的有限元整机结构模型。通过非线性静力分析获得接触层在工作载荷下的接触刚度。其中，接触刚度包含固体部分的接触刚度和基于接触变形求得的流体刚度，此时的接触刚度为非线性。由于整机结构在工作时的载荷波动远小于静载荷，可以认为结构工作时的接触层刚度保持不变，这样非线性的有限元模型可以线性化。最后，静力分析获得的接触层刚度可以导入到模态分析模型中，这样得到的包含虚拟材料接触层的有限元整机模态分析模型是一个线性模型。通过对整机线性模型的模态分析可以获得包含接触层影响的整机振型和固有频率。6.4试验结果对比为了验证本文建立模型的有效性，将有限元仿真固有频率与试验获得的固有频率进行对比。表 4 中对比了拧紧力矩分别为10，20，30 Nm的有限元仿真固有频率与试验获得的固有频率。由表 4 可以看出，含润滑介质条件下的试验和有限元仿真获得的固有频率比不含润滑介质条件下的试验和有限元仿真获得的固有频率更高。这是因为在含润滑介质的正交各向异性结合面的法向接触总刚度的计算中，不仅考虑了微凸体的接触刚度，还考虑了润滑介质的刚度。图 20 含虚拟接触层的有限元模型Fig.20 Finite element model containing virtual contact layer图 21 有限元仿真振型图Fig.21 Vibration shapes from FEA图 22 有限元仿真计算流程Fig.22 Calculation process of finite element simulation919振 动 工 程 学 报第 36 卷表 5给出了含润滑介质条件下有限元仿真固有频率与试验固有频率的相对误差。由表 5 可以看出，在螺栓分别承受力矩为10，20，30 N m工况下，含润滑介质条件下有限元仿真固有频率与试验固有频率的相对误差在0.65%4.11%之间。图 23是根据含润滑介质、无润滑介质情况下固有频率之间的相对误差拟合出的曲线。由图 23 可以看出，不论是试验测量还是有限元仿真获得的固有频率都随着预紧力的增大而增大，这是因为增大螺栓的预紧力会使结合面的实际接触面积增大，导致结合面的接触刚度变大，从而使系统的固有频率增大。7结 论（1）对影响结合面的法向接触刚度的不同因素分析得出：当结合面刚发生接触时，流体刚度对结合面的法向接触刚度的影响较为明显，随着接触压力升高，固体接触刚度成为影响结合面的法向接触刚度的主要因素。图 23 不同预紧力下前三阶固有频率对比Fig.23 Comparison of the first three order natural frequencies under different preloads表 5 含润滑介质条件下有限元仿真固有频率与试验固有频率的相对误差Tab.5 Relative error of finite element simulation natural frequency and experimental natural frequency with lubricated media力矩/(Nm)102030阶数123123123含润滑介质条件下的相对误差/%3.841.060.651.433.531.494.112.410.68表 4 有限元仿真固有频率与试验获得的固有频率对比Tab.4 Comparison of finite element simulation natural frequency and experimental natural frequency力矩/(Nm)102030阶数123123123含润滑介质试验/Hz91012211688104813902134109515332201无润滑介质试验/Hz828109314539371265201598414212062含润滑介质有限元/Hz87512081677103313412102105014962186无润滑介质有限元/Hz7981055149096112362006101513862087920第 4 期王世军，等：含润滑介质的正交各向异性结合面接触特性研究（2）含润滑介质的正交各向异性结合面的法向接触刚度K*n与 D，及呈单调递增关系，与G*呈单调递减关系。（3）在含润滑介质的正交各向异性结合面的法向接触总刚度的计算中，不仅考虑了微凸体的接触刚度，还考虑了润滑介质的影响，因此，含润滑介质条件下的试验和有限元仿真获得的固有频率，比不含润滑介质条件下的试验和有限元仿真获得的固有频率有显著提高，表明整机分析中有必要考虑结合面中润滑介质的影响。（4）含润滑介质条件下，有限元仿真的前三阶振型与模态试验的前三阶振型一致，有限元仿真的固有频率与模态试验的固有频率相比，最大相对误差为 4.11%，证明本文提出的含润滑介质的正交各向异性结合面的法向接触刚度模型是合理的，可以较准确地预测结合面的接触性能。参考文献:1王世军，赵金娟.机械工程中的有限元方法 M.北京：科学出版社，2019：239-241.Wang Shijun，Zhao Jinjuan.Finite Element Methods in Mechanical EngineeringM.Beijing：Science Press，2019：239-241.2 Chung J C，Lin J F.Fractal model developed for elliptic elastic-plastic asperity microcontacts of rough surfacesJ.Journal of Tribology，2004，126（4）：646-654.3Buczkowski R，Kleiber M.Elasto-plastic statistical model of strongly anisotropic rough surfaces for finite element 3D-contact analysisJ.Computer Methods in Applied Mechanics and Engineering，2006，195（37-40）：5141-5161.4王世军，赵金娟，张慧军，等.一种结合部法向刚度的预估方法 J.机械工程学报，2011，47（21）：111-115.WANG Shijun，ZHAO Jinjuan，ZHANG Huijun，et al.A prediction method for the normal stiffness of the joint J.Journal of Mechanical Engineering，2011，47（21）：111-115.5刘伟强，张进军，洪军，等.椭圆抛物体形微凸体弹性接 触 力 学 模 型J.西 安 交 通 大 学 学 报，2015，49（10）：34-40.LIU Weiqiang，ZHANG Jinjun，HONG Jun，et al.Elastic contact mechanics model for elliptical parabolic asperity bodiesJ.Journal of Xian Jiaotong University，2015，49（10）：34-40.6谭文兵，兰国生，张学良，等.依据各向异性分形理论的固定结合面椭圆弹塑性法向接触刚度建模及仿真分析 J.固体力学学报，2021，42（1）：63-76.TAN Wenbing，LAN Guosheng，ZHANG Xueliang，et al.Modeling and simulation analysis of elliptical normal contact stiffness of fixed bond surfaces based on anisotropic fractal theoryJ.Chinese Journal