CN106529035A - 一种基于三维分形理论的栓接结构动态特性分析方法 - Google Patents

Abstract

一种基于三维分形理论的栓接结构动态特性分析方法，现代工业飞速发展，对机床加工精度要求也越来越高，而机床部件间结合部接触质量对机床加工精度有很大影响，研究表明，机床总刚度的60％～80％和总阻尼的90％，来自于结合部，而机床的整体刚度与阻尼又决定了机床的加工精度，因此准确预估结合部接触刚度与阻尼能够为研究机床整体动态性能提供理论基础。本发明针对栓接结合面提供了一种基于三维分形表面同时考虑了弹塑性变形阶段的接触刚度阻尼模型，通过MATLAB数值仿真软件与ANSYS有限元分析软件相结合预估栓接结合面动态性能的方法。

Description

一种基于三维分形理论的栓接结构动态特性分析方法

技术领域

本发明涉及到一种基于三维分形理论的栓接结构动态特性分析方法，属于结构动态性能研究领域，具体涉及通过建立三维分形模型估计结合面刚度阻尼特性参数，进而达到预估结构整体动态性能的分析方法。

背景技术

栓接结合面是机械零件中常见的一种结合面，而结合面的接触质量在很大程度上影响着机械结构的整体性能，研究表明，结合面刚度和阻尼均影响着机床结构整体的刚度和阻尼，其接触质量成为影响机床加工精度的重要因素之一。因此如何通过准确估计结合部表面参数达到预估结构整体动态性能的研究是具有深远意义的，尤其在现代的机床结构设计中，精确辨识结合面参数(包括接触刚度和接触阻尼)是对机床整体进行准确建模、性能预测和优化设计的关键环节。

20世纪60年代，人们认识到结合面接触质量对机床结构整体的重要影响后，国内外学者开始进行结合面动力特性研究，半个世纪以来，已经有多种识别估计结合面参数的方法，早期人们直接测试结合面处的力与位移的关系来研究结合面刚度，这种识别方法最直接却因每次实验无法考虑多种影响参数而受到限制。此外，由J.H.Wang和C.M.Liou提出的频响函数法是一种通过利用频响函数来提取结合面参数的方法，该方法同事考虑了频响函数测试时误差的影响。后来，人们通过有限元理论对微观结合面参数进行建模的方法对结构整体动态性能进行分析，通过此方法实现对结构性能的预测和结构的优化。

早期基于有限元理论的微接触方法，国外学者研究了多种利用结合部表面测量数据通过数学建模来描述结合面参数的方法。1966年，Greenwood与Williamson共同提出了基于统计分析的接触模型，即G-W模型。这种模型首次考虑表面形貌参数，将表面的高度分布视为随机变量，对接触理论的研究具有重要影响。此后Whitehouse和Archard又共同建立了W-A模型，虽然较G-W模型有所改善，但依旧是基于统计分析的基础上，这就使得在表面接触分析中，其各统计学参数都受到仪器分辨率的影响。

自1990年A.Majumdar和C.L.Tien利用W-M函数对粗糙表面进行分形表征，A.Majumdar和B.Bhushan在此基础上建立了M-B分形接触模型，这种模型与G-W模型有着相同的假设，将粗糙表面之间的接触简化为粗糙表面与刚性理想平面的接触(典型的微凸体接触模型如图1所示)，其最大特点是，粗糙表面微观形貌特征的表征参数—分形维数D和分形特征长度尺寸参数G具有尺度独立性，而不受仪器分辨率的影响，并与取样长度无关。各国学者就一直在这个研究内容上进行改进、完善和扩展。该方法考虑了结合面压力分布不均的情况，能够较准确地估计结合面刚度阻尼特性参数，对动态性能预测和结构优化具有深远的意义。

发明内容

本发明旨在提供一种基于三维分形理论的栓接结构动态特性分析方法。该方法的主要特点是通过在分形理论中同时考虑弹塑性变形与域拓展因子的影响建立了精确的结合面刚度阻尼模型；基于不均匀分布压力和分形刚度阻尼模型，在结合面间建立刚度阻尼单元矩阵，从而实现对栓接结构动态特性的精确预估。

本发明是采用以下技术手段实现的：

S1、基于M-B分形模型，同时考虑弹塑性变形及域拓展因子的影响，建立精确的三维分形刚度阻尼模型。

S2、选取栓接结构模型进行ANSYS建立实体模型并进行静力分析获取栓接结合面的压力分布云图，并提取结合面各节点处的压力值。

S3、与S1建立的数学分形模型相结合对栓接结构模型进行自由模态分析以获取结构整体的各阶固有频率和振型图，并对结构变形进行分析预测。

S4、在以上模态分析的基础之上进行谐响应分析，并提取结构上某一点的频响图。

本发明的特点是通过在分形理论中同时考虑弹塑性变形与域拓展因子的影响建立了精确的结合面刚度阻尼模型；基于不均匀分布压力和分形刚度阻尼模型，在结合面间建立刚度阻尼单元矩阵，从而实现对栓接结构动态特性的精确预估。

附图说明

图1典型的微凸体接触模型图；

图2栓接结合部动态性能分析流程图；

图3典型栓接结构件示例；

图4栓接结合面约束及压力分布图；

图5频率与幅值图；

具体实施方式

以下结合附图2-5对本发明作进一步详细说明。

步骤(1)建立三维结合面刚度阻尼模型

基于M-B分形理论，结合赫兹理论，同时考虑弹塑性变形和域拓展因子ψ，通过对处于不同变形区域的单个微凸体法向载荷进行积分可得到总弹性法向载荷、弹塑性法向载荷及总塑性法向载荷分别如下：

其中H_G1,H_G2均为与材料属性及结合面分形参数相关的系数，H为较软材料的硬度，H＝2.8Y(Y为屈服强度值)；k为与泊松比相关的参数，k＝0.454+0.41ν；a′_1c,a′_2c分别为弹性、弹塑性及塑性变形间临界横截面积。

则结合面总法向载荷可表示为F＝F_e+F_ep+F_p。

法向刚度建模中，在弹性变形与弹塑性变形区域法向接触刚度分别为：

切向刚度建模中，在弹性变形与弹塑性变形区域切向接触刚度分别为：

式中G'为结合部等效剪切模量，1/G'＝(2-ν₁)/G₁+(2-ν₂)/G₂；H₁,H₂通过公式得到。

综上，结合面法向与切向总接触刚度分别为：K_n＝K_ne+K_nep，K_t＝K_te+K_tep。

法向阻尼建模中，结合面在弹性和塑性变形阶段应变能分别为：

则法向阻尼因子为法向阻尼为

切向阻尼建模中，结合面在弹性和弹塑性变形阶段周期总能量W₁,W与耗散能量W_d1,W_d2分别为：

式中H₃通过公式得到，

则切向阻尼因子为切向阻尼为

步骤(2)有限元静力分析及节点压力值提取

在建立模型结合面部分时，采用了TARGET170单元和CONTACT174单元作为目标单元和接触单元，进行结合面部分的静力分析，为了避免分析复杂并缩短计算时间，栓接接触部分网格划分密度较密集，其他部分较粗糙的方法，以提高分析效率。

步骤(3)进行自由模态分析得到结构固有频率及振型

将静力分析提取的结合部接触面的节点压强与单元面积的乘积作为结合面总法向载荷代入第一步建立的数学分形模型中通过MATLAB进行计算，以获得每一个结合面接触部分各节点单元的刚度及阻尼值。将所得的各节点的等效刚度和阻尼通过ANSYS MATRIX27单元与各节点一一对应建立刚度和阻尼单元，并进行自由模态分析，以得到结构的固有频率和振型。

步骤(4)进行谐响应分析得到结构频响图

在进行谐响应分析时，对结构进行激励，并求取在该激励下对应点的频响函数。

Claims (2)

1.一种基于三维分形理论的栓接结构动态特性分析方法，其特征在于：

S1、基于M-B分形模型，同时考虑弹塑性变形及域拓展因子的影响，建立精确的三维分形刚度阻尼模型；

S2、选取栓接结构模型进行ANSYS建立实体模型并进行静力分析获取栓接结合面的压力分布云图，并提取结合面各节点处的压力值；

S3、与S1建立的数学分形模型相结合对栓接结构模型进行自由模态分析以获取结构整体的各阶固有频率和振型图，并对结构变形进行分析预测；

S4、在以上模态分析的基础之上进行谐响应分析，并提取结构上某一点的频响图。

2.根据权利要求1所述的一种基于三维分形理论的栓接结构动态特性分析方法，其特征在于：

步骤(1)建立三维结合面刚度阻尼模型；

基于M-B分形理论，结合赫兹理论，同时考虑弹塑性变形和域拓展因子ψ，通过对处于不同变形区域的单个微凸体法向载荷进行积分可得到总弹性法向载荷、弹塑性法向载荷及总塑性法向载荷分别如下：

$F_e=\left\{\begin{array}{cc}\frac{2^{(11-2D)/2}}{3{\mathrm{\π}}^{(4-D)/2}}\·\frac{D-1}{5-2D}{\left(ln\mathrm{\γ}\right)}^{1/2}{G}^{(D-2)}{\mathrm{E\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}(a_l^{\′(5-2D)/2}-{a_{1c}^{\′}}^{(5-2D)/2})& (D\≠2.5)\\ 2{\mathrm{\π}}^{-3/4}{\left(ln\mathrm{\γ}\right)}^{1/2}{G}^{1/2}{\mathrm{E\ψ}}^{1/4}{a}_l^{\′3/4}ln\left(\frac{a_l^{\′}}{a_{1c}^{\′}}\right)& (D=2.5)\end{array}\right.$

$F_{ep}=\left\{\begin{array}{cc}\frac{H_{G2}(D-1){\mathrm{\ψ}}^{(3-D)/2}{a}_l^{(D-1)/2}}{2(2.26-0.88D)}\·({a_{1c}^{\′}}^{(2.26-0.88D)}-{a_{2c}^{\′}}^{(2.26-0.88D)})& D\≠113/44\\ \frac{69}{88}{H}_{G2}{\mathrm{\ψ}}^{19/88}{a}_l^{\′69/88}ln\left(\frac{a_{1c}^{\′}}{a_{2c}^{\′}}\right)& D=113/44\end{array}\right.$

$F_p=\frac{H(D-1)a_l^{\′(D-1)/2}}{3-D}{\mathrm{\ψ}}^{(3-D)/2}{a}_{2c}^{\′(3-D)/2}$

其中H_G1,H_G2均为与材料属性及结合面分形参数相关的系数，

H为较软材料的硬度，H＝2.8Y(Y为屈服强度值)；k为与泊松比相关的参数，k＝0.454+0.41ν；a′_1c,a′_2c分别为弹性、弹塑性及塑性变形间临界横截面积；

则结合面总法向载荷可表示为F＝F_e+F_ep+F_p；

法向刚度建模中，在弹性变形与弹塑性变形区域法向接触刚度分别为：

$K_{ne}=\frac{2\sqrt{2}E(4-D)(D-1)}{3\sqrt{\mathrm{\π}}(3-D)(2-D)}{\mathrm{\ψ}}^{(3-D)/2}{a_l}^{\′(D-1)/2}({a_l}^{\′(2-D)/2}-{a_{1c}}^{\′(2-D)/2})$

$k_{nep}=\frac{{\mathrm{df}}_{ep}/{\mathrm{da}}^{\′}}{d\mathrm{\δ}/{\mathrm{da}}^{\′}}=\frac{H_{G2}{\mathrm{\π}}^{(3-D)/2}(1.76-0.38D)}{2^{3-D}{G}^{(D-2)}{\left(ln\mathrm{\γ}\right)}^{1/2}(3-D)}{a}^{\′(0.26+0.12D)}$

切向刚度建模中，在弹性变形与弹塑性变形区域切向接触刚度分别为：

$K_{te}=\frac{4G^{\′}(D-1)}{\sqrt{2\mathrm{\π}}}{\mathrm{\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}{\∫}_{a_{1c}^{\′}}^{{a_l}^{\′}}{a}^{\′(-D/2)}{(1-H_1{a}^{\′(D-2)/2})}^{1/3}{\mathrm{da}}^{\′}$

$K_{tep}=\frac{4G^{\′}(D-1)}{\sqrt{\mathrm{\π}}{H_{G2}}^{1/2}}{\mathrm{\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}{\∫}_{a_{1c}^{\′}}^{{a_l}^{\′}}{a}^{\′(-D/2)}{(1-H_2{a}^{\′(0.44-0.22D)})}^{1/3}{\mathrm{da}}^{\′}$

式中G'为结合部等效剪切模量，1/G'＝(2-ν₁)/G₁+(2-ν₂)/G₂；H₁,H₂通过公式得到；

综上，结合面法向与切向总接触刚度分别为：K_n＝K_ne+K_nep，K_t＝K_te+K_tep；

法向阻尼建模中，结合面在弹性和塑性变形阶段应变能分别为：

$W_e=\left\{\begin{array}{cc}\frac{2^{(19-4D)/2}{\mathrm{EG}}^{2(D-2)}\left(ln\mathrm{\γ}\right)(3-D)(D-1)}{3{\mathrm{\π}}^{(7-2D)/2}(7-2D)(8-3D)}{\mathrm{\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}(a_l^{\′(8-3D)/2}-{a_{1c}^{\′}}^{(8-3D)/2})& D\≠8/3\\ \frac{2^{19/6}{\mathrm{EG}}^{4/3}\left(ln\mathrm{\γ}\right)}{3{\mathrm{\π}}^{5/6}}{\mathrm{\ψ}}^{1/6}{a}_l^{\′5/6}ln\left(\frac{a_l^{\′}}{a_{1c}^{\′}}\right)& D=8/3\end{array}\right.$

$W_p=\frac{2^{(3-D)}{\mathrm{HG}}^{(D-2)}{\left(ln\mathrm{\γ}\right)}^{1/2}(D-1)}{(5-D){\mathrm{\π}}^{(3-D)/2}}{\mathrm{\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}{a_{2c}^{\′}}^{(3-D)}$

则法向阻尼因子为法向阻尼为

切向阻尼建模中，结合面在弹性和弹塑性变形阶段周期总能量W₁,W与耗散能量W_d1,W_d2分别为：

$W_1=\frac{3{\mathrm{\π\μ}}^2{H}_3(D-1){\mathrm{\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}}{16G^{\′}}{\∫}_{a_{1c}^{\′}}^{a_l^{\′}}\[\frac{3}{5}+{(1-H_1{a}^{\′(D-2)/2})}^{3/5}-{(1-H_1{a}^{\′(D-2)/2})}^{2/3}\]a^{\′(7-3D)/2}{\mathrm{da}}^{\′}$

$\begin{array}{c}{W}_2=\frac{3{\mathrm{\π}}^{3/2}{\mathrm{\μ}}^2{H_{G1}}^{0.5}{H_{G2}}^2(D-1){\mathrm{\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}}{16G^{\′}}\\ \×{\∫}_{a_{2c}^{\′}}^{a_{1c}^{\′}}\[\frac{3}{5}+{(1-H_2{a}^{\′(0.44-0.22D)})}^{3/5}-{(1-H_2{a}^{\′(0.44-0.22D)})}^{2/3}\]a^{\′(1.92-0.96D)}{\mathrm{da}}^{\′}\end{array}$

$\begin{array}{c}{W}_{d1}=\frac{9{\mathrm{\π}}^{1/2}{\mathrm{\μ}}^2{H}_3(D-1){\mathrm{\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}}{\mathrm{20}{G}^{\′}}\\ \×{\∫}_{a_{1c}^{\′}}^{a_l^{\′}}\{1-{(1-H_1{a}^{\′(D-2)/2})}^{5/3}-\frac{5H_1}{6}{a}^{\′(D-2)/2}\[1-{(1-H_1{a}^{\′(D-2)/2})}^{2/3}\]\}{a}^{\′(7-3D)/2}{\mathrm{da}}^{\′}\end{array}$

$\begin{array}{c}{W}_{d2}=\frac{9{\mathrm{\π}}^{1/2}{\mathrm{\μ}}^2{H_{G1}}^{0.5}{H_{G2}}^2(D-1){\mathrm{\ψ}}^{(3-D)/2}{a}_l^{\′(D-1)/2}}{20G^{\′}}\\ \×{\∫}_{a_{2c}^{\′}}^{a_{1c}^{\′}}\{1-{(1-H_2{a}^{\′(0.44-0.22D)})}^{5/3}-\frac{5}{6}{H}_2{a}^{\′(0.44-0.22D)}{\[1-(1-H_2{a}^{\′(0.44-0.22D)})\]}^{2/3}\}{a}^{\′(1.92-0.96D)}{\mathrm{da}}^{\′}\end{array}$

式中H₃通过公式得到，

则切向阻尼因子为切向阻尼为

步骤(2)有限元静力分析及节点压力值提取

在建立模型结合面部分时，采用了TARGET170单元和CONTACT174单元作为目标单元和接触单元，进行结合面部分的静力分析，为了避免分析复杂并缩短计算时间，栓接接触部分网格划分密度较密集，其他部分较粗糙的方法，以提高分析效率；

步骤(3)进行自由模态分析得到结构固有频率及振型

将静力分析提取的结合部接触面的节点压强与单元面积的乘积作为结合面总法向载荷代入第一步建立的数学分形模型中通过MATLAB进行计算，以获得每一个结合面接触部分各节点单元的刚度及阻尼值；将所得的各节点的等效刚度和阻尼通过ANSYS MATRIX 27单元与各节点一一对应建立刚度和阻尼单元，并进行自由模态分析，以得到结构的固有频率和振型；

步骤(4)进行谐响应分析得到结构频响图

在进行谐响应分析时，对结构进行激励，并求取在该激励下对应点的频响函数。