接触刚度单元划分密度对计算精度的影响

接触刚度单元划分密度对计算精度的影响
作者：李初晔王海涛马岩
摘要：有限元分析中通常采用接触单元和耦合节点位移两种方法处理组合体计算问题，对于明显存在相对滑动的组合结构，采用节点耦合的计算方法已不适用，必须采用接触单元来处理此类问题。

接触计算时零件之间的接触穿透会丧失部分计算精度，穿透值的大小主要与接触刚度有关，对于弯曲问题，单元划分密度也会对穿透值的大小产生影响。

本文基于ANSYS 的接触计算理论，讨论了接触刚度、单元划分密度对有限元计算精度的影响，期望找到提高接触计算精度的方法。

前言
零件组合结构分析的难点是如何正确处理零件连接部位之间的约束过渡，通过此约束条件的建立完成工作力在整机零件之间的传递。

有限元分析软件通常采用接触单元和耦合节点位移两种方法处理此类问题，节点耦合是点对点的位移协调关系，即一组耦合节点在变形过程中沿所定义的自由度方向变形保持一致。

显然只要两点距离足够接近，并且各零件在变形过程中没有相对滑动，这种定义是符合实际的。

然而对于明显存在相对滑动的组合结构，采用节点耦合的计算方法已不适用，必须采用接触单元来处理此类问题。

接触计算时零件之间的接触穿透会丧失部分计算精度，穿透值的大小主要与接触刚度有关，对于弯曲问题，单元划分密度会影响最大穿透值的测量基准位置，间接对穿透值的大小产生影响。

本文的研究内容基于ANSYS的接触计算理论，讨论了接触刚度、单元划分密度对有限元计算精度的影响，期望找到提高接触计算精度的方法。

1 问题描述
图1的两板之间采用接触单元连接，力F使其压缩，板厚度a，弹性模量Ｅ，泊淞比c，板间间隙e，研究接触刚度及单元划分密度对计算精度的影响。

有限元分析采用的模型参数：L1=800mm，L2=1200mm，a=300mm，b=800mm，e=8mm。

材料特性参数：弹性模量210000MPa，泊淞比0.3。

A，B两板之间存在8mm间隙，采用接触单元连接，板A在力F作用下逐渐靠近板B。

在开始阶段，由于接触单元与目标单元的距离大于接触发生的公差范围，接触单元处于“开”状态，A，B之间不发生力的传递，这与标准弹簧有本质区别；当距离处于接触公差范围带内，接触单元处于“闭合”状态，A，B之间开始传递力，此后，当接触面和目标面出现空间位置交错时，发生接触穿透现象。

2 ANSYS的接触穿透计算
ANSYS的接触刚度系数(FKN)与穿透公差系数(Ftoln)用来计算有限元分析时实际采用的接触刚度ＫN。

若接触单元的附着层为实体单元，单元厚度h，则允许的最大穿透Tmax为：
Tmax=Ftoln*h/2 （1）
若接触单元的附着层为壳或面单元，单元实常数(厚度)为h，则允许的最大穿透Tmax为：
Tmax=Ftoln*10*(4h)/2 （2）
式中，4h为单元四个节点的实常数厚度和，在此假设四个节点为等厚度h。

ANSYS计算采用的实际接触刚度ＫN为：
KN＝E*FKN/Tmax（3）
E为附着层单元的弹性模量。

若定义的接触单元总面积为S，接触力F，对实体单元平均穿透值T为：
T=F/(KN*S)=F*Ftoln*h/(2*E*FKN*S) （4）
h为附着层单元厚度，对于实体单元实际等于单元划分长度。

若附着层单元为壳或面单元，平均穿透值T为：
T=F/(KN*S)=F*Ftoln*10*4h/(2*E*FKN*S) （5）
h为单元实常数厚度。

从公式（4）、（5）看出，在同样的作用力下，刚度系数FKN越小，穿透越大；接触面积S越小，穿透越大；材料越软（既弹性模量E越小），穿透越大。

对实体单元，单元划分越稀疏，穿透越大；对面或壳单元，实常数厚度越大，则穿透越大。

从公式中看出，对于壳或面单元，穿透值与单元划分长度（疏密程度）没有直接关系，与厚度h有关。

但是大量工程计算表明，随着单元划分长度的减小使穿透值出现增大的现象，后续弯曲模型计算论述了这一结论。

图2为接触刚度系数0.3，单元划分长度50mm时的计算结果，图3是穿透值与位移的关系曲线。

接触处于闭合状态后，穿透值、接触力与C点位移呈近似线形关系。

当C点位移为10mm时，D点位移为1.235mm，接触力为5.23x10^7N，穿透值范围为0.008~0.0116mm，平均穿透值在0.009mm左右。

C，D点的位置见图1。

当接触刚度系数为0.3，单元划分尺度设定为30mm时，接触力为5.23x10^7N，穿透值范围为0.0049~0.0076mm，平均穿透值在0.0055mm左右，见图4。

接触刚度系数为0.3，单元划分长50mm时，由公式（5）计算出的平均穿透值：
T=F*Ftoln*h/(2*E*FKN*S)=5.23x10^7*0.1*50/(2*210000*0.3*800*300)=0.0087(mm)
接触刚度系数为0.3，单元划分长30mm时，由公式（5）计算出的平均穿透值：
T=F*Ftoln*h/(2*E*FKN*S)=5.23x10^7*0.1*30/(2*210000*0.3*800*300)=0.0052(mm)
有限元计算值与公式计算值有很好的一致性。

刚度系数和单元划分密度会影响接触穿透值的大小，间接对接触力等主要计算结果产生影响，在拉压问题中，这种影响是可以忽略的，但弯曲问题穿透过大造成的影响往往不容忽视。

由前面的公式可以计算出，在5.23x10^7N力作用下，当单元尺寸为50mm，刚度系数为0.0026时，接触穿透达到1mm，但与刚度系数0.3时的有限元计算结果相比，A，B板的应力和变形几乎没有什么影响。

但对于弯曲问题，接触过大的穿透会严重影响计算精度，单元划分密度会影响最大穿透值的测量基准位置，间接对穿透值的大小产生影响。

3 接触穿透对有限元计算精度的影响
对于纯拉压问题，接触穿透对计算结果的影响可以忽略。

但对于弯曲问题，接触过大的穿透会严重影响计算精度。

图5是弯曲接触问题计算模型，两种梁的刚度等效，左端梁在中间处沿垂直方向有位移约束，右端梁在相同位置采用柱面支撑，在端部力作用下，梁与柱发生接触。

在小变形条件下，接触近似于线接触，因此两种模型的弯曲变形应近似相等。

非接触等效结构梁支撑处反力为：左端752N，中间1252N，最大变形为13.151mm。

接触刚度系数为0.3，单元划分尺寸2mm时，接触结构梁支撑处反力为：左端736N，中间1236N。

最大变形为13.521mm，接触面的最大接触穿透为0.0785mm。

两种梁在同样载荷作用下，最大变形相差3%左右，由于接触穿透的影响使接触梁的变形大于非接触梁变形。

采用ANSYS的APDL编程，选取多种刚度系数和单元划分长度进行计算，得到图7~图10的结果。

图7是选用的接触刚度系数0.3，计算得到的接触穿透值随单元划分尺寸关系曲线，随着单元划分尺寸的增大，接触穿透越来越小，这一结论对接触结构分析至关重要，从公式（5）中看出，对于壳或面单元，穿透值与单元划分长度没有直接关系，与厚度h有关。

但是，由于单元划分长度的减小，弯曲结构真实的接触面积S减小，从而使穿透值出现增大现象。

图8是接触刚度系数0.3时两种结构的变形随单元划分尺寸的变化关系，非接触模型的变形基本不随单元尺寸变化的影响，接触模型的变形随单元尺寸的增大逐渐增大，但总体变化幅度不大，与非接触模型的结果相比，计算精度是可以接受的。

图9是单元划分尺寸
2mm时穿透值与刚度系数的关系曲线，随着接触刚度增加，穿透值逐渐减小。

但增大刚度系数可能会使计算发散，同时会加大求解时间，因此应综合考虑计算精度和计算时间的平衡关系。

图10是单元划分尺寸2mm时接触模型梁变形与刚度系数的关系曲线，随着接触刚度增加，接触梁变形逐渐逼近非接触梁变形值。

因此如果计算机性能允许，尽可能选择较大的接触刚度系数、较密的网格划分计算，一般接触刚度系数取0.3~0.7就可得到很高的计算精度。

4 结论
基于ANSYS接触计算理论，得到了平均穿透值的计算公式，通过对拉压及弯曲接触模型的有限元计算验证了公式的准确性。

分析表明：刚度系数和单元划分密度会影响接触穿透值的大小，间接对接触力等主要计算结果产生影响，在拉压问题中，这种影响是可以忽略的，但对于弯曲问题，接触过大的穿透会严重影响计算精度，单元划分密度会影响最大穿透值的测量基准位置，间接对穿透值的大小产生影响。

随着单元划分尺寸的增大，接触穿透越来越小，这一结论对接触结构分析至关重要，从公式（5）中看出，对于壳或面单元，穿透值与单元划分长度没有直接关系，与厚度h有关。

但是，由于单元划分长度的减小，弯曲结构真实的接触面积S减小，从而使穿透值出现增大现象。

图8是接触刚度系数0.3时两种结构的
变形随单元划分尺寸的变化关系，非接触模型的变形基本不随单元尺寸变化的影响，接触模型的变形随单元尺寸的增大逐渐增大，但总体变化幅度不大，与非接触模型的结果相比，计算精度是可以接受的。

图9是单元划分尺寸2mm时穿透值与刚度系数的关系曲线，随着接触刚度增加，穿透值逐渐减小。

但增大刚度系数可能会使计算发散，同时会加大求解时间，因此应综合考虑计算精度和计算时间的平衡关系。

图10是单元划分尺寸2mm时接触模型梁变形与刚度系数的关系曲线，随着接触刚度增加，接触梁变形逐渐逼近非接触梁变形值。

因此如果计算机性能允许，尽可能选择较大的接触刚度系数、较密的网格划分计算，一般接触刚度系数取0.3~0.7就可得到很高的计算精度。