一种基于曲面拟合的确定结合面接触面积及刚度的方法

摘要

一种基于曲面拟合的确定结合面接触面积及刚度的方法，先采用共聚焦显微镜及三坐标测量来测量机械结合面微凸体形貌及表面波纹度和形状误差，得到的点云数据用二元高次函数进行拟合,从而得到接触表面最终解析式，然后单对微凸体接触点位置判断及接触方向判断，在进行单个接触点弹Hertz接触计算，计算在力的作用下接触变形与接触面积，最后进行接触面整体接触面积与接触刚度计算，得出总接触面积与各方向接触刚度，与传统解析法相比，本方法有更接近真实形貌的优点；与有限元法相比，本方法可扩大可计算的接触面积。

说明书

技术领域

本发明属于机械接合面力学领域，具体涉及一种基于曲面拟合的 确定结合面接触面积及刚度的方法。

背景技术

在机床、燃气轮机、汽车发动机等复杂机械系统中存在大量的结 合面，它们破坏了机械系统结构的连续性，在很大程度上影响了机械 系统的整机性能。结合面的连接特性呈现出一定的非线性，将结合面 特性参数模型引入系统整机建模过程可以有效地提高产品设计阶段 的整机性能预测能力。因此，准确构建结合面特性参数模型将为复杂 机械部件结合面优化设计提供理论基础。

结合面处的真实接触面积与应力分布是决定机械结合面电阻、热 阻、接触刚度等特性的根本原因。计算的准确程度将直接影响结合面 参数特性模型的建立。结合面两个接触表面在微观尺度下是粗糙不平 的，通常认为是由宏观形状误差、波纹度和粗糙度三种尺度的表面叠 加而成，因此当两个机加工表面相互接触时，结合面处的接触行为仅 发生在一些离散的微凸体上。而这些微凸体的接触状态决定了结合面 处的电阻热阻等特性。因此，有必要提出一种数学方法来对粗糙表面 建模，以精确地计算结合面接触刚度与接触面积，来实现结合面性能 的准确预测。

由于粗糙表面形貌复杂，微凸体形状不规则，目前国内外学者在 研究结合面微观接触行为时采用的方法是：用球或旋转抛物体等规则 几何体模拟粗糙峰，将实际的微凸体接触行为简化为这些规则几何体 的接触行为。这种对微观接触行为的研究是建立在许多假设条件的情 况下进行的，虽然模型一直在发展、完善，但现有的接触模型基本上 仍都存在以下问题：1、微凸体形状假定为球状或其他简单形状，脱 离实际；2、在研究整个表面的接触行为时，假设表面轮廓高度分布 服从某种函数分布如正态分布，采用表面粗糙度统计参数如平均曲率 半径等,并没有从真实表面出发；3、局限于在粗糙尺度层面研究结合 面处接触刚度，没有考虑几何形状即波纹度、宏观形状误差对结合面 接触刚度的影响规律。现有的有限元方法可采用真实的表面测量数据， 来进行接触问题的研究，但因为有限元计算需要划分极细的网格，严 重限制了计算面积的大小及计算效率。

发明内容

为了克服上述现有技术的缺点，本发明的目的在于提供一种基于 曲面拟合的确定结合面接触面积及刚度的方法，能够极大扩展计算面 积，并提高计算效率。

为了达到上述目的，本发明采取的技术方案为：

一种基于曲面拟合的确定结合面接触面积及刚度的方法，包括以 下步骤：

1)真实形貌的测量及函数表示

1.1)结合面形貌层次划分：将机械结合面按波长由小到大区分 为粗糙度、波纹度、形状误差三个层次，将峰与峰之间的间距小于 1mm的起伏视作粗糙度，相应单个凸起视作微凸体,并假设微凸体形 状相同；将波峰与波峰间的距离在1mm到10mm范围内的表面起伏视 作波纹度；将波峰与波峰间的距离超过10mm范围的表面起伏当作表 面形状误差处理；

1.2)微凸体函数拟合：采用共聚焦显微镜观测微凸体形貌，获 得的点云数据用来拟合微凸体，用二元高次函数来进行拟合，按刀具 加工轨迹方向，将加工轨迹方向定为u向，垂直于加工轨迹的方向定 为v向，每个微凸体拟合为:

$f_{\mathrm{ai}}=f_{\mathrm{ai}}(u,v)=a_i+\underset{t=0}{\overset{1}{\Sigma }}{b}_i{u}^t{v}^{1-t}+\underset{t=0}{\overset{2}{\Sigma }}{c}_i{u}^t{v}^{2-t}+......,((u,v)\in {\Omega }_i,i=\mathrm{1,2,3}......)$

其中：i表示粗糙表面上的第i个微凸体，Ω_i表示第i个微凸体所占 的表面区域；

1.3)波纹度及形状误差拟合：采用三坐标测量仪来测量接触面 的波纹度及形状误差，获得的点云数据用二元高次函数来进行拟合， 表面形状误差及波纹度拟合为：

$f_e=f_e(u,v)=a_e+\underset{i=0}{\overset{1}{\Sigma }}{b}_e{u}^i{v}^{1-i}+\underset{i=0}{\overset{2}{\Sigma }}{c}_e{u}^i{v}^{2-i}+......$

1.4)最终曲面的函数表达：将微凸体函数分布到表面形状误差 函数上相加，得到最终真实表面形貌的函数表达式：

$f=f(u,v)=\left(\begin{array}{cc}{f}_{a1}(u,v)+f_e(u,v)& ((u,v)\in {\Omega }_1)\\ f_{a2}(u,v)+f_e(u,v)& ((u,v)\in {\Omega }_2)\\ f_{a3}(u,v)+f_e(u,v)& ((u,v)\in {\Omega }_3)\\ ........& \end{array}\right)$

2)两接触表面接触点位置判断及单对微凸体接触的预处理

2.1)两表面间间距函数表达式：将两个接触表面函数相减，得到 接触间距公式：

δ＝f₁(u，v)-f₂(u，v)

其中，f₁(u，v)表示上表面拟合函数，f₂(u，v)表示下表面拟合函数，

2.2)接触点位置获取：接触点(u_c，v_c)满足的条件是δ(u_c，v_c)为函数极 小值点，并且δ(u_c，v_c),根据极小值求法即可求得所有可能的接触点， 再根据δ(u_c，v_c)＜0判断，即可得出所有接触点；

2.3)接触点方向求解及接触点坐标变换：设接触点的空间坐标 为(x(u_c，v_c)，y(u_c，v_c)，z(u_c，v_c))，接触方向按微分几何的方法进行求解：

$\overrightarrow{r_u^{\prime }}=\left(\begin{array}{ccc}{x}_{u_c}^{\prime }& y_{u_c}^{\prime }& z_{u_c}^{\prime }\end{array}\right),\overrightarrow{r_v^{\prime }}=\left(\begin{array}{ccc}{x}_{v_c}^{\prime }& y_{v_c}^{\prime }& z_{v_c}^{\prime }\end{array}\right),\overrightarrow{n}(u_c,v_c)=\overrightarrow{r_u^{\prime }}\times \overrightarrow{r_v^{\prime }}=\overrightarrow{n}{(x_n,y_n,z_n)}^T$

其中分接触点的两个切向量，为接触点处的法向量，之后 将接触点处附近的曲面进行坐标变换，使得切平面方向为新坐标轴的 xoy平面，接触点为新坐标的原点，将接触点处的两表面用泰勒级数 展开，忽略高阶无穷小，保留函数的二次方项，从而得到两表面在新 坐标系下的近似表达式：

$\left(\begin{array}{c}{z}_1=A_1{\bar{x}}^2+B_1{\bar{y}}^2+C_1\overline{\mathrm{xy}}\\ z_2=A_2{\bar{x}}^2+B_2{\bar{y}}^2+C_2\overline{\mathrm{xy}}\end{array}\right)$

两式相减，并通过坐标系旋转变换，消除xy交叉项，得接触间隔的表 达式为：

$z_1-z_2=A{\bar{x}}^2+B{\bar{y}}^2$

其中，A，B都是常数，经过这样的转化就可使得单对微凸体的接 触解析公式满足Hertz接触计算的要求；

3)单对接触点处的Hertz接触力学计算

3.1)接触间距系数求解：根据Hertz任意形状曲面接触理论，对 于椭圆任意点，A，B可由两微凸体接触点处的主曲率及主曲率方向来 求得，用R,R′分别代表其中一个微凸体在接触点处的两个主曲率，w代 表两和微凸体接触点间主曲率方向夹角，系数A，B满足下式所表示的 方程，己知各微凸体主曲率大小及其主曲率方向间的夹角，即可求得 A，B的具体数值：

$\left(\begin{array}{c}A+B=\frac{1}{2}(\frac{1}{R_1}+\frac{1}{R_1^{\prime }}+\frac{1}{R_2}+\frac{1}{R_2^{\prime }})\\ B-A=\frac{1}{2}{\left(\begin{array}{c}{(\frac{1}{R_1}-\frac{1}{R_1^{\prime }})}^2+{(\frac{1}{R_2}-\frac{1}{R_2^{\prime }})}^2+\\ 2(\frac{1}{R_1}-\frac{1}{R_1^{\prime }})(\frac{1}{R_2}-\frac{1}{R_2^{\prime }})\cos 2w\end{array}\right)}^{\frac{1}{2}}\end{array}\right)$

3.2)接触区域长短半轴求解：复杂曲面点接触的接触区域近似 为椭圆，长短半轴分别为：

$\left(\begin{array}{c}a=m\sqrt[3]{\frac{3\pi }{4}\frac{P(k_1+k_2)}{(A+B)}}\\ b=n\sqrt[3]{\frac{3\pi }{4}\frac{P(k_1+k_2)}{(A+B)}}\end{array}\right)$

其中m，n是与比率(B-A)/(A+B)相关的系，数值采用接触理论的文献 Jamari J,Schipper D J.An elastic–plastic contact model of ellipsoid  bodies[J].Tribology letters,2006,21(3):262-271，a，b分别为接触区域 的长短半轴，k₁+k₂为接触点处的综合弹性模量，k由以下公式表示：

其中，E，v分别代表材料的弹性模量和泊松比，

3.3)接触点处接触法向量方向上的位移求解：垂直于接触方向 的位移定为两个方向接触位移的均值：

3.4)最大接触应力与接触面积求解：最大接触应力与接触面积 可由以下方程得出：

$\left(\begin{array}{c}{q}_0=\frac{3}{2}\frac{P}{\mathrm{\pi ab}}\\ s=\mathrm{\pi ab}\end{array}\right)$

其中P代表单对接触点所受的力，q₀为接触点处的最大接触应力， s为单对接触点处的接触面积，

4)结合面总真实接触面积与法向切向接触刚度计算

4.1)总接触面积计算：总的接触面积等于各个接触点接触面积 的总和，

4.2)各方向接触刚度计算：将各点处的所受的力累加可得最终 的合力并按各个坐标轴方向进行分解，得到垂直于整个接触 面的法向分力P_sumz及平行于接触面的两个分力P_sumx及P_sumy，各点处的位 移按矢量分解并求均值可得各方向上的平均位移某方向上的 刚度数值，按以下公式即可求出：

法向接触刚度的表达式为:两个切向方向上的刚度为: $k_z=\frac{\partial P_{\mathrm{sumx}}}{\partial {\bar{d}}_x},k_z=\frac{\partial P_{\mathrm{sumy}}}{\partial {\bar{d}}_y}.$

本发明有以下有益效果：充分考虑了结合面的形状误差的三维形 貌及微凸体真实三维形状对接触面积及接触刚度的影响，从实际粗糙 面形貌入手，摆脱了原先方法基于大量假设的缺点，与有限元法相比， 能够极大扩展计算面积，并提高计算效率。

附图说明

图1是机械结合面粗糙表面的分类示意图。

图2是单对微凸体接触及其坐标变换示意图。

图3是Hertz曲面接触接触示意图。

具体实施方式

下面结合附图对本发明进行详细描述：

一种基于曲面拟合的确定结合面接触面积及刚度的方法，包括以 下步骤：

1)真实形貌的测量及函数表示

1.1)结合面形貌层次划分：如图1所示，将机械结合面按波长 由小到大区分为粗糙度、波纹度、形状误差三个层次，将峰与峰之间 的间距小于1mm的起伏视作粗糙度，相应单个凸起视作微凸体,并假 设微凸体形状相同；将波峰与波峰间的距离在1mm到10mm范围内的 表面起伏视作波纹度；将波峰与波峰间的距离超过10mm范围的表面 起伏当作表面形状误差处理；

1.2)微凸体函数拟合：采用共聚焦显微镜观测微凸体形貌，获 得的点云数据用来拟合微凸体，用二元高次函数来进行拟合，按刀具 加工轨迹方向，将加工轨迹方向定为u向，垂直于加工轨迹的方向定 为v向，每个微凸体拟合为:

$f_{\mathrm{ai}}=f_{\mathrm{ai}}(u,v)=a_i+\underset{t=0}{\overset{1}{\Sigma }}{b}_i{u}^t{v}^{1-t}+\underset{t=0}{\overset{2}{\Sigma }}{c}_i{u}^t{v}^{2-t}+......,((u,v)\in {\Omega }_i,i=\mathrm{1,2,3}......)$

其中i表示粗糙表面上的第i个微凸体，Ω_i表示第i个微凸体所占的 表面区域；

1.3)波纹度及形状误差拟合：采用三坐标测量仪来测量接触面的 波纹度及形状误差，获得的点云数据用二元高次函数来进行拟合，表 面形状误差及波纹度拟合为：

$f_e=f_e(u,v)=a_e+\underset{i=0}{\overset{1}{\Sigma }}{b}_e{u}^i{v}^{1-i}+\underset{i=0}{\overset{2}{\Sigma }}{c}_e{u}^i{v}^{2-i}+......$

1.4)最终曲面的函数表达：将微凸体函数分布到表面形状误差 函数上相加，得到最终真实表面形貌的函数表达式，

$f=f(u,v)=\left(\begin{array}{cc}{f}_{a1}(u,v)+f_e(u,v)& ((u,v)\in {\Omega }_1)\\ f_{a2}(u,v)+f_e(u,v)& ((u,v)\in {\Omega }_2)\\ f_{a3}(u,v)+f_e(u,v)& ((u,v)\in {\Omega }_3)\\ ........& \end{array}\right)$

2)两接触表面接触点力学判断及单对微凸体接触的预处理

2.1)两表面间间距函数表达式：将两个接触表面函数相减，得到 接触间距公式：

δ＝f₁(u，v)-f₂(u，v)

其中，f₁(u，v)表示上表面拟合函数，f₂(u，v)表示下表面拟合函数，

2.2)接触点位置获取：接触点(u_c，v_c)满足的条件是δ(u_c，v_c)为函数极 小值点，并且δ(u_c，v_c)＜0,根据极小值求法即可求得所有可能的接触点， 再根据δ(u_c，v_c)＜0判断，即可得出所有接触点；

2.3)接触点方向求解及接触点坐标变换：设接触点的空间坐标 为(x(u_c，v_c)，y(u_c，v_c)，z(u_c，v_c))，接触方向按微分几何的方法进行求解：

$\overrightarrow{r_u^{\prime }}=\left(\begin{array}{ccc}{x}_{u_c}^{\prime }& y_{u_c}^{\prime }& z_{u_c}^{\prime }\end{array}\right),\overrightarrow{r_v^{\prime }}=\left(\begin{array}{ccc}{x}_{v_c}^{\prime }& y_{v_c}^{\prime }& z_{v_c}^{\prime }\end{array}\right),\overrightarrow{n}(u_c,v_c)=\overrightarrow{r_u^{\prime }}\times \overrightarrow{r_v^{\prime }}=\overrightarrow{n}{(x_n,y_n,z_n)}^T$

其中分接触点的两个切向量，为接触点处的法向量，之后 将接触点处附近的曲面进行坐标变换，如图2所示，使得切平面方向 为新坐标轴的xoy平面，接触点为新坐标的原点，将接触点处的两表 面用泰勒级数展开，忽略高阶无穷小，保留函数的二次方项，从而得 到两表面在新坐标系下的近似表达式：

$\left(\begin{array}{c}{z}_1=A_1{\bar{x}}^2+B_1{\bar{y}}^2+C_1\overline{\mathrm{xy}}\\ z_2=A_2{\bar{x}}^2+B_2{\bar{y}}^2+C_2\overline{\mathrm{xy}}\end{array}\right)$

两式相减，并通过坐标系旋转变换，消除xy交叉项，得接触间隔的表 达式为：

$z_1-z_2=A{\bar{x}}^2+B{\bar{y}}^2$

其中，A，B都是常数，经过这样的转化就可使得单对微凸体的接 触解析公式满足Hertz接触计算的要求；

3)单对接触点处的Hertz接触力学计算

3.1)接触间距系数求解：根据Hertz任意形状曲面接触理论，对 于椭圆任意点，A，B可由两微凸体接触点处的主曲率及主曲率方向来 求得，如图3所示，用R,R′分别代表其中一个微凸体在接触点处的两 个主曲率，w代表两和微凸体接触点间主曲率方向夹角，系数A，B满 足下式所表示的方程，己知各微凸体主曲率大小及其主曲率方向间的 夹角，即可求得A，B的具体数值：

$\left(\begin{array}{c}A+B=\frac{1}{2}(\frac{1}{R_1}+\frac{1}{R_1^{\prime }}+\frac{1}{R_2}+\frac{1}{R_2^{\prime }})\\ B-A=\frac{1}{2}{\left(\begin{array}{c}{(\frac{1}{R_1}-\frac{1}{R_1^{\prime }})}^2+{(\frac{1}{R_2}-\frac{1}{R_2^{\prime }})}^2+\\ 2(\frac{1}{R_1}-\frac{1}{R_1^{\prime }})(\frac{1}{R_2}-\frac{1}{R_2^{\prime }})\cos 2w\end{array}\right)}^{\frac{1}{2}}\end{array}\right)$

3.2)接触区域长短半轴求解：复杂曲面点接触的接触区域近似 为椭圆，长短半轴分别为：

$\left(\begin{array}{c}a=m\sqrt[3]{\frac{3\pi }{4}\frac{P(k_1+k_2)}{(A+B)}}\\ b=n\sqrt[3]{\frac{3\pi }{4}\frac{P(k_1+k_2)}{(A+B)}}\end{array}\right)$

其中m，n是与比率(B-A)/(A+B)相关的系，数值采用接触理论的文献 Jamari J,Schipper D J.An elastic–plastic contact model of ellipsoid bodies[J].Tribology letters,2006,21(3):262-271，a，b分别为接触区 域的长短半轴。k₁+k₂为接触点处的综合弹性模量，k由以下公式表 示：

其中，E，v分别代表材料的弹性模量和泊松比，

3.3)接触点处接触法向量方向上的位移求解：垂直于接触方向 的位移定为两个方向接触位移的均值：

3.4)最大接触应力与接触面积求解：最大接触应力与接触面积 可由以下方程得出：

$\left(\begin{array}{c}{q}_0=\frac{3}{2}\frac{P}{\mathrm{\pi ab}}\\ s=\mathrm{\pi ab}\end{array}\right)$

其中P代表单对接触点所受的力，q₀为接触点处的最大接触应力， s为单对接触点处的接触面积，

4)结合面总真实接触面积与法向切向接触刚度计算

4.1)总接触面积计算：总的接触面积等于各个接触点接触面积 的总和，

4.2)各方向接触刚度计算：将各点处的所受的力累加可得最终 的合力并按各个坐标轴方向进行分解，得到垂直于整个接触 面的法向分力P_sumz及平行于接触面的两个分力P_sumx及P_sumy，各点处的位 移按矢量分解并求均值可得各方向上的平均位移某方向上的 刚度数值，按以下公式即可求出：

法向接触刚度的表达式为:两个切向方向上的刚度为: $k_z=\frac{\partial P_{\mathrm{sumx}}}{\partial {\bar{d}}_x},k_z=\frac{\partial P_{\mathrm{sumy}}}{\partial {\bar{d}}_y}.$