一种圆柱滚子轴承简化有限元模型的快速建模方法
简介
简介
目前，传动系统仍然处于持续开发和快速迭代的阶段，以满足未来对更大载荷、更长使用寿命以及降低噪音和重量的需求。为了有效地实现这些目标，采用诸如有限元方法（FEM）等计算机辅助工程（CAE）技术，是较为常见的手段。
可靠的有限元分析依赖于对传动系统在运行过程中承受的载荷有深入了解，同时必须正确模拟这些载荷对传统系统的影响。由于轴承在传动系统中起到传递载荷的作用，因此轴承是分析传动系统时必须考虑的关键部件。
尽管建立具有实体滚子和考虑接触行为的轴承精细有限元模型完全可行，但对于包含多个轴承的传动系统，或是具有大量滚子的大型轴承，从计算成本的角度来看，采用这种方式并非最佳选择。这是因为，对于球形滚子或圆柱形滚子，滚子与滚道之间的弹性接触可近似为点接触或线接触。由于接触面积相比于轴承的尺寸非常小，为了准确模拟滚子与滚道之间的接触行为，必须采用非常细密的网格来划分，这将严重增加计算成本。
在开展传动系统的有限元分析时，如果关注区域并非轴承本身，而是与轴承相关的结构。此时，如果能够建立复杂度适中、数值稳定、经过充分测试和验证的简化轴承模型，使得简化轴承模型能够提供正确的结构整体刚度以及载荷分布，预计可以获得足够准确的计算结果，并节约计算成本。
为此，本研究将针对圆柱滚子轴承，讨论其简化有限元模型的建模方法，并基于Abaqus二次开发，给出一种快速建模方法，能够根据有限的参数快速生成模型，以进一步提高建模效率。
圆柱滚子轴承简化有限元模型建模方法
圆柱滚子轴承简化有限元模型建模方法
模型分类
目前，在建立圆柱滚子轴承的简化有限元模型时，普遍的做法为采用实体单元建立轴承的外圈和内圈，并通过具有非线性刚度的弹簧来代替滚子。这种简化建模方式的优点为可以较好地反映轴承外圈和内圈的刚度，通过调整弹簧的布置方式和刚度，可以获得与真实轴承接近的载荷-变形行为。
根据弹簧布置方式的不同，上述简化模型可以分为如图1所示的两类模型。在图中，左侧的轴承模型可以称为分布载荷（Distributed Load）模型，右侧的轴承模型可以称为离散载荷（Discrete Load）模型。
图1 两类简化轴承模型
在分布载荷模型中，外圈滚道和内圈滚道的所有节点均通过弹簧单元连接；而在离散载荷模型中，只有滚子所在区域的节点通过弹簧连接。对于分布载荷模型，当轴承承受载荷作用时，载荷将会被均匀地分配到内圈和外圈上，从而获得更平滑的应力分布，但这也意味着此类模型无法研究滚子和滚道之间的作用载荷。对于离散载荷模型，由于只有滚子所在区域的节点采用弹簧连接，因此相比之下更接近真实轴承的承载情况，并且能够获得外载作用下滚子和滚道之间的作用载荷的分布特性。
对于圆柱滚子轴承，由于滚子与滚道之间的接触行为可以近似为线接触，因此对于每一个滚子，外圈和内圈滚道沿轴向方向上的节点均通过一排弹簧单元连接，如图2所示。通过布置多个弹簧，可以反映单个滚子在轴向方向的载荷分布情况。
图2 单个滚子的节点连接示意图
弹簧非线性刚度计算
在上述简化轴承模型中，弹簧刚度可以通过滚子的载荷-变形行为计算得到。由于随着滚子承受的载荷增大，滚子与滚道之间的接触面积也会随之增大。因此，即使是线弹性模型，滚子的载荷-变形行为也是非线性的。大量的研究人员通过开展试验或仿真，得到对于单个滚子与滚道之间的载荷-变形关系可采用如下所示的表达式近似描述：
式中，Q为滚子载荷，δ为变形量，K为载荷-变形系数，n为一个指数，对于滚子轴承可近似取值为1.1。
由于内圈和外圈均会产生变形，因此内圈滚道相对外圈滚道的总的径向变形量δtot可表示为：
式中，δi和δo分别为内圈和外圈的变形量。
由于滚子受到内圈和外圈的作用力而处于平衡状态，因此有：
则对于一个与内圈和外圈滚道同时接触的滚子，其载荷-变形关系可表示为：
式中，Kn可表示为
式中，Ki和Ko分别为内圈和外圈的载荷-变形系数。载荷-变形系数取决于滚子的有效长度L，也同样受到滚子直径D的影响。因此，滚子-滚道的载荷-变形可写为如下所示的一般形式：
式中，L0、D0和δ0用于对参数进行归一化，当参数以毫米表示时，其取值均为1。n1、n2和n3为常数，基于部分文献的研究成果，其取值如表1所示。
表1 计算滚子-滚道接触的载荷-变形系数使用的常数
需要说明的是，只有滚道与滚子处于受压状态时，滚道与滚子才会处于接触状态；如果滚道与滚子处于受拉状态时，则滚道将与滚子脱开。这意味着，代表滚子的弹簧只能受压，而不能受拉。上述载荷-变形行为可以通过将弹簧处于拉伸状态时的刚度设置为零来实现，此时弹簧的载荷-变形行为如图3所示。在图中，左侧的载荷-变形曲线代表轴承的游隙为零时的情况，右侧的载荷-变形曲线代表轴承的游隙不为零时的情况。此时，可以通过将左侧的载荷-变形曲线向左侧偏移游隙值来获得。其基本原理为当代表滚子的弹簧受压时，当弹簧产生的压缩量小于游隙时，弹簧不会受到载荷的作用；只有当变形量超过游隙值后，弹簧才会承受载荷的作用。此外，对于圆柱滚子轴承，在沿着滚子的轴向方向上，滚子边缘的接触应力将大于滚子中心的接触应力。为了获得较为均匀的接触应力分布，圆柱滚子表现为沿着轴向方向上滚子直径发生改变。对于单个滚子中沿轴向方向上的弹簧，通过调整每个弹簧所代表的轴承游隙值，可以考虑滚子直径的变化。
图3 非线性弹簧的载荷-变形行为
此外，从图3中可以看到，当弹簧受拉时，弹簧的刚度值并不为零，而是具有一个非常微小的刚度值。通过该方法，可以消除轴承内圈变形时的刚体 位移，从而避免在静力分析时出现收敛性问题。需要特别说明的是，在定义非线性弹簧的载荷-变形行为时，需要在较宽的变形范围内定义。这是因为，对于超出数据定义范围内的变形量，有限元软件Abaqus并不会对数据进行线性外推处理，而是假设载荷保持不变，这可能会带来收敛性问题。
另一个需要特别说明的问题为非线性弹簧刚度的标定。这是因为采用上述理论公式计算得到的弹簧刚度包含了滚道的变形。由于有限元模型中已经通过建立内圈和外圈的实体模型考虑了滚道的变形，因此直接使用理论公式计算的刚度会带来额外的误差。此外，由于弹簧是直接通过连接外圈和内圈滚道的节点创建的，而真实的滚子与滚道在变形时存在一定的接触面积，因此同样会存在误差。一些研究人员通过计算滚子与滚道的接触宽度，在接触区域创建一层非常薄的单元层来考虑真实的接触面积。然而，这种方法 会引入长宽比差异较大的单元，同样可能会带来一定的收敛性问题。综上所述，为了获得接近真实轴承的载荷-变形行为，弹簧的刚度应该略大于理论计算值。一种较为可行的方法为通过建立包含单个实体滚子的轴承有限元模型或试验的方式来获得真实轴承的载荷-变形行为，并基于该载荷-变形行为来标定弹簧的刚度。由于本文仅讨论圆柱滚子轴承简化有限元模型的建模方法，因此仍使用理论公式计算得到的弹簧刚度。关于简化轴承有限元模型中弹簧刚度对于轴承滚子载荷分布以及周围结构的影响，仍有待进一步研究。
圆柱滚子轴承简化有限元模型快速建模方法
圆柱滚子轴承简化有限元模型快速建模方法
对于圆柱滚子轴承，通过将滚子所在部位的外圈和内圈节点采用具有非线性刚度的弹簧单元依次连接，即可建立圆柱滚子轴承简化有限元模型。然而，对于具有较多滚子数量的轴承结构，如果采用手动连接的方式，建模过程将非常繁琐。为此，本文基于有限元软件Abaqus的二次开发功能，通过批量建立非线性弹簧单元，提出一种圆柱滚子轴承简化有限元模型的快速建模方法。
该方法的基本原理为将外圈滚道和内圈滚道的所有表面节点分别定义两个节点集。通过定义变量来指定两个节点集的名称，如变量Outer_Ring_Node_Set_Name和变量Inner_Ring_Node_Set_Name，可以编写脚本获得两个节点集中所有节点对应的节点编号以及坐标值。然后，通过指定轴承滚子的数量和第一个轴承滚子所在位置（指定滚子中心位置的X和Y坐标，假定轴承轴向为Z轴），脚本程序将会依照节点的极坐标值确定所有轴承滚子对应区域的外圈和内圈滚道表面节点，并依照滚子的编号存放在不同的节点集中。例如，对于第一个滚子，其所在区域的外圈滚道表面节点将会存放在节点集Outer_Ring_Roller1_Node_Set中，内圈滚道表面节点将会存放在节点集Inner_Ring_Roller1_Node_Set中。其余滚子的节点集名称可以依次类推。
在建立每个滚子所在区域节点的节点集后，通过将单个滚子的外圈节点集中的全部节点与内圈节点集中的全部节点通过弹簧单元连接即可建立单个滚子的弹簧单元。需要指出的是，由于节点集中的节点并非沿着轴向顺序排列，因此需要首先判断节点的连接顺序。本文给出的解决方法为首先获取节点集中节点的轴向坐标（Z坐标）按从小到大的顺序排列时，其在节点集中的索引值（即节点对象在节点集对象中的排列顺序，如第一个节点的索引值为0）。然后，根据该索引值来连接外圈和内圈节点，即可确保节点连接顺序无误。
在通过节点集中节点的索引值获得待连接的节点后，只需在两个节点之间创建弹簧单元即可。需要说明的是，对于有限元软件Abaqus，在图形界面中创建的单元仅支持恒定弹簧刚度，如果需要创建非线性刚度的弹簧单元，则必须手动修改inp文件。并且，修改之后的inp文件只能直接进行求解，如果再次导入到Abaqus中，则弹簧的非线性刚度将会失效。此外，由于需要修改inp文件，因此这种方法不利于基于Abaqus二次开发进行交互处理。
为此，本文采用连接器单元来创建具有非线性刚度的弹簧。有限元软件Abaqus提供的连接器单元支持创建各种复杂的连接关系，如铰接、衬套、万向节等。在创建类似于弹簧的连接器单元时，需要在Interaction模块中将连接器截面（Connector Section）指定为轴向（Axial），如图4所示。
图4 具有弹簧变形行为的连接器单元
由该截面定义的连接器单元仅有轴向自由度，而轴向方向由连接器单元两个端点构成的矢量来定义。由于该单元没有额外的转动自由度，因此该单元在本质上可以认为是Abaqus中的桁架（Truss）单元。两者的区别为桁架单元的变形由弹性模量E定义，而这里的连接器单元的变形行为由连接器的刚度定义。
在连接器截面的定义菜单中继续点击Continue即可定义连接器的变形行为。在定义时，选择行为选项（Behavior Options）为弹性（Elasticity），并设置弹性定义为非线性（Nonlinear），如图5。
图5 连接器非线性刚度定义
此时，连接器的非线性刚度由连接器单元的载荷（F or M）-变形（U or UR）关系来指定。注意，连接器单元的变形以受拉为正，以受压为负，这与弹簧单元的定义正好相反。通过定义如图3所示的载荷-变形关系，即可建立具有非线性刚度的连接器单元。需要说明的是，由于这里定义了非线性刚度，因此在开展静力分析时，必须打开几何非线性（Nlgeom）选项，否则可能会得到错误的变形行为。
在完成连接器单元的截面定义后，只需通过Interaction模块的Connector Builder选项定义连接器两个端点的节点，并赋予该截面，即可完成非线性刚度的连接器单元创建。在本文中，上述过程将基于Abaqus的Python脚本自动实现，并批量创建所有滚子对应的连接器单元。
此外，脚本程序将针对每一个连接器单元自动创建时间历程输出，以输出连接器作用的轴向载荷（Connector Total Force，CTF1）。这里的轴向载荷即为轴承滚子承受的载荷。在计算完成后，基于本文提供的后处理脚本，可以批量获取每一个滚子所在区域的连接器载荷，从而得到轴承滚子的载荷分布情况。
关于圆柱滚子简化有限元模型快速建模方法的具体操作流程，可参见本文在附录提供的建模和后处理脚本，以及操作视频。
计算示例
计算示例
考虑如图6所示的一个圆柱滚子轴承，假设滚子数量为50。
图6 圆柱滚子轴承几何尺寸示意图
基于上述简化方法，将滚子采用连接器单元代替，因此仅划分内圈和外圈实体网格，滚子局部网格如图7所示。
图7 圆柱滚子轴承局部网格示意图
为了便于加载，将轴承外圈的外表面节点全部施加固定约束，在轴承中心创建参考点，将参考点与轴承内圈的表面节点耦合。约束参考点除Y向平动自由度以外的全部自由度，并在参考点上施加沿Y向的集中载荷F，取值为1000 kN。采用上述方法批量建立连接器单元，得到圆柱滚子轴承简化有限元模型如图8所示。
图8 圆柱滚子轴承简化有限元模型示意图
计算完成后，采用上述后处理脚本提取每个连接器上的载荷，并将同一个滚子上所有连接器的载荷累加，即可得到作用于单个滚子上的载荷。通过统计所有滚子上作用的载荷，可以得到圆柱滚子轴承的载荷分布如图9所示。
图9 圆柱滚子轴承的载荷分布
在图9中，蓝色曲线代表将轴承外圈节点全部固定时的轴承载荷分布。作为对比，图中的橙色曲线代表将轴承外圈节点部分固定时的轴承载荷分布。可以看到，在径向载荷的作用下，轴承上部的滚子并不会受到载荷的作用，载荷将主要集中在下半部分。在这里，部分固定和完全固定轴承外圈节点可视为轴承周围结构的局部刚度对轴承载荷分布的影响。不难看出，当轴承连接的局部结构的刚度较小时，轴承载荷的分布范围更广，并且载荷峰值将会略微降低。通过本文提出的方法，建立考虑完整结构和简化轴承的有限元模型，可以考虑大载荷作用时结构刚度对于轴承载荷分布的影响。
模型局限性
模型局限性
现有模型在模拟真实轴承的载荷-变形行为时仍存在一定的局限性。如上述模型尚无法考虑轴承承受轴向力时的情况。在此类情况下，需要在外圈和内圈滚道之间建立梁单元，以模拟滚子的轴向刚度。此外，现有模型将外圈和内圈滚道节点直接通过具有非线性刚度的连接器单元连接。因此，不能实现轴承外圈和内圈的相对转动。一种可能的解决方法为将内圈节点与孤立节点连接，然后建立孤立节点与外圈节点之间的接触关系，由该方法可以在考虑滚子刚度的情况下，模拟轴承外圈和内圈之间的转动。
附录
附录
圆柱滚子轴承简化有限元模型建模脚本（Python）
#-* - coding:UTF-8 -*-
#------------------------------------------------------------------
#文件名: Bearing_Connector_Generator
#程序用于批量生成轴承简化有限元模型所需要连接器单元
#当前程序假设轴承轴向方向为Z向
#作者: Tian Xu, Southwest Jiaotong University
#------------------------------------------------------------------
from abaqus import *
from abaqusConstants import *
from caeModules import *
import math
#--------------------------程序输入参数定义--------------------------
#模型名称
Model_Name='Main_Bearing'
#存放外圈滚道所有节点的集 合名
Outer_Ring_Node_Set_Name='Outer_Ring_Node'
#存放内圈滚道所有节点的集 合名
Inner_Ring_Node_Set_Name='Inner_Ring_Node'
#定义连接器载荷历程输出的分析步名称
Step_Name='Step-1'
#轴承滚子数量
Roller_Num=50
#第一个滚子的圆心坐标
Roller1_Pos_Coord_X=0.0
Roller1_Pos_Coord_Y=-653.0
#滚子刚度定义数据
Roller_Stiffness_Data_U=[-1.0, 0.0, 1.0] #滚子刚度定义数据-变形值
Roller_Stiffness_Data_F=[-7.0e5,0.0,7.0e5] #滚子刚度定义数据-载荷值
#-------------------------------------------------------------------
#-------------------生成连接器截面信息(滚子非线性刚度)-----------------
My_Model=mdb.models[Model_Name]
#创建连接器截面(连接器类型:Axial 活动自由度:U1)
My_Connector_SEC=My_Model.ConnectorSection(name='Bearing_Srping_SEC',
translationalType=AXIAL)
#设置连接器非线性刚度
Stiffness_Table=tuple([(Roller_Stiffness_Data_F[i],
Roller_Stiffness_Data_U[i]) for i in range(len(Roller_Stiffness_Data_U))])
My_Stiffness_Data=connectorBehavior.ConnectorElasticity(components=(1, ),
behavior=NONLINEAR, table=Stiffness_Table)
My_Connector_SEC.setValues(behaviorOptions=(My_Stiffness_Data,))
#-------------------------------------------------------------------
#-------------------输出轴承内圈和外圈滚道所有节点的信息---------------
#获取存放外圈滚道所有节点的节点集对象
My_Outer_Ring_Node_Set=My_Model.rootAssembly.sets[Outer_Ring_Node_Set_Name]
#获取存放内圈滚道所有节点的节点集对象
My_Inner_Ring_Node_Set=My_Model.rootAssembly.sets[Inner_Ring_Node_Set_Name]
#创建列表存放节点坐标信息
Outer_Ring_Node_Set_Coord_X=[]
Outer_Ring_Node_Set_Coord_Y=[]
Outer_Ring_Node_Set_Coord_Z=[]
Inner_Ring_Node_Set_Coord_X=[]
Inner_Ring_Node_Set_Coord_Y=[]
Inner_Ring_Node_Set_Coord_Z=[]
#读取外圈滚道所有节点的节点号
Outer_Ring_Node_ID_List=[]
for node in My_Outer_Ring_Node_Set.nodes:
Outer_Ring_Node_ID_List.append(node.label)
coords=node.coordinates
Outer_Ring_Node_Set_Coord_X.append(coords[0])
Outer_Ring_Node_Set_Coord_Y.append(coords[1])
Outer_Ring_Node_Set_Coord_Z.append(coords[2])
#读取内圈滚道所有节点的节点号
Inner_Ring_Node_ID_List=[]
for node in My_Inner_Ring_Node_Set.nodes:
Inner_Ring_Node_ID_List.append(node.label)
coords=node.coordinates
Inner_Ring_Node_Set_Coord_X.append(coords[0])
Inner_Ring_Node_Set_Coord_Y.append(coords[1])
Inner_Ring_Node_Set_Coord_Z.append(coords[2])
#输出外圈和内圈滚道所有节点的节点号至屏幕
#print('Node Set '+Outer_Ring_Node_Set_Name+ ':')
#print(Outer_Ring_Node_ID_List)
#print('Node Set '+Inner_Ring_Node_Set_Name+ ':')
#print(Inner_Ring_Node_ID_List)
#-------------------------------------------------------------------
#------------------------建立滚子所在节点的节点集---------------------
#创建列表存放节点坐标信息(极坐标形式)
Outer_Ring_Node_Set_PolarCoord_R=[]
Outer_Ring_Node_Set_PolarCoord_Theta=[]
Outer_Ring_Node_Set_PolarCoord_Z=[]
Inner_Ring_Node_Set_PolarCoord_R=[]
Inner_Ring_Node_Set_PolarCoord_Theta=[]
Inner_Ring_Node_Set_PolarCoord_Z=[]
#计算外圈和内圈滚道所有节点的极坐标值
for i in range(len(Outer_Ring_Node_ID_List)):
X=Outer_Ring_Node_Set_Coord_X[i]
Y=Outer_Ring_Node_Set_Coord_Y[i]
Z=Outer_Ring_Node_Set_Coord_Z[i]
#计算极径
R=math.sqrt(X**2.0+Y**2.0)
#计算角度(从-pi到pi)
Theta=math.atan2(Y,X)
#变换角度(从0到2pi)
if0.0<=Theta<=math.pi:
Theta=Theta
else:
Theta=Theta+2.0*math.pi
#存放计算结果
Outer_Ring_Node_Set_PolarCoord_R.append(R)
Outer_Ring_Node_Set_PolarCoord_Theta.append(Theta)
Outer_Ring_Node_Set_PolarCoord_Z.append(Z)
for i in range(len(Inner_Ring_Node_ID_List)):
X=Inner_Ring_Node_Set_Coord_X[i]
Y=Inner_Ring_Node_Set_Coord_Y[i]
Z=Inner_Ring_Node_Set_Coord_Z[i]
#计算极径
R=math.sqrt(X**2.0+Y**2.0)
#计算角度(从-pi到pi)
Theta=math.atan2(Y,X)
#变换角度(从0到2pi)
if0.0<=Theta<=math.pi:
Theta=Theta
else:
Theta=Theta+2.0*math.pi
#存放计算结果
Inner_Ring_Node_Set_PolarCoord_R.append(R)
Inner_Ring_Node_Set_PolarCoord_Theta.append(Theta)
Inner_Ring_Node_Set_PolarCoord_Z.append(Z)
#格式化输出外圈和内圈全部节点的笛卡尔和极坐标信息至外部文件(程序调试用)
Outer_Ring_Node_Set_Info_File=open('Outer_Ring_Node_Set_Info.txt','w')
for index, cx, cy, cz, pr, ptheta, pz in zip(
Outer_Ring_Node_ID_List,
Outer_Ring_Node_Set_Coord_X,
Outer_Ring_Node_Set_Coord_Y,
Outer_Ring_Node_Set_Coord_Z,
Outer_Ring_Node_Set_PolarCoord_R,
Outer_Ring_Node_Set_PolarCoord_Theta,
Outer_Ring_Node_Set_PolarCoord_Z):
Outer_Ring_Node_Set_Info_File.write(
'%6d %.6f %.6f %.6f %.6f %.6f %.6f\n' %(index, cx, cy, cz, pr, ptheta, pz))
Outer_Ring_Node_Set_Info_File.close()
Inner_Ring_Node_Set_Info_File=open('Inner_Ring_Node_Set_Info.txt','w')
for index, cx, cy, cz, pr, ptheta, pz in zip(
Inner_Ring_Node_ID_List,
Inner_Ring_Node_Set_Coord_X,
Inner_Ring_Node_Set_Coord_Y,
Inner_Ring_Node_Set_Coord_Z,
Inner_Ring_Node_Set_PolarCoord_R,
Inner_Ring_Node_Set_PolarCoord_Theta,
Inner_Ring_Node_Set_PolarCoord_Z):
Inner_Ring_Node_Set_Info_File.write(
'%6d %.6f %.6f %.6f %.6f %.6f %.6f\n' %(index, cx, cy, cz, pr, ptheta, pz))
Inner_Ring_Node_Set_Info_File.close()
#筛选外圈滚子所在区域的节点
#计算第一个滚子所在位置的极角(从-pi到pi)
Roller1_Pos_Coord_Theta=math.atan2(Roller1_Pos_Coord_Y,Roller1_Pos_Coord_X)
#变换第一个滚子所在位置的极角(从0到2pi)
if0.0<=Roller1_Pos_Coord_Theta<=math.pi:
Roller1_Pos_Coord_Theta=Roller1_Pos_Coord_Theta
else:
Roller1_Pos_Coord_Theta=Roller1_Pos_Coord_Theta+2.0*math.pi
#计算筛选滚子节点所用的角度容差(滚子数量的100倍)
Pi_Val=math.pi
Angle_Tol=2.0*Pi_Val/(Roller_Num*1000.0)
for Roller_Index in range(1,Roller_Num+1):
#计算当前滚子所在位置的准确角度值
Roller_Angle=0.0+2.0*Pi_Val/Roller_Num*(Roller_Index-1)
#初始化存放滚子所在区域节点的列表
Outer_Ring_Roller_Node=[]
#计算考虑角度容差时角度筛选的上限和下限值
Roller_Angle_Max_Lim=Roller_Angle+Angle_Tol
Roller_Angle_Min_Lim=Roller_Angle-Angle_Tol
#基于角度筛选符合要求的节点
for i in range(len(Outer_Ring_Node_ID_List)):
#计算将节点变换到以第一个滚子为起始位置时的角度(局部坐标)
Local_CSYS_PolarCoord_Theta=\
Outer_Ring_Node_Set_PolarCoord_Theta[i]+2.0*math.pi-Roller1_Pos_Coord_Theta
if Local_CSYS_PolarCoord_Theta>=2.0*math.pi:
Local_CSYS_PolarCoord_Theta=Local_CSYS_PolarCoord_Theta-2.0*math.pi
#筛选获得滚子所在区域的节点
if Roller_Angle_Min_Lim<=Local_CSYS_PolarCoord_Theta<=Roller_Angle_Max_Lim:
Outer_Ring_Roller_Node.append(My_Outer_Ring_Node_Set.nodes[i])
#将节点列表转换为MeshNodeArray对象
Outer_Ring_Roller_Node=mesh.MeshNodeArray(Outer_Ring_Roller_Node)
#创建新的节点集存放外圈滚子的节点
My_Model.rootAssembly.Set(nodes=Outer_Ring_Roller_Node,
name='Outer_Ring_Roller'+str(Roller_Index)+'_Node_Set')
#筛选内圈滚子所在区域的节点
for Roller_Index in range(1,Roller_Num+1):
#计算当前滚子所在位置的准确角度值
Roller_Angle=0.0+2.0*Pi_Val/Roller_Num*(Roller_Index-1)
#初始化存放滚子所在区域节点的列表
Inner_Ring_Roller_Node=[]
#计算考虑角度容差时角度筛选的上限和下限值
Roller_Angle_Max_Lim=Roller_Angle+Angle_Tol
Roller_Angle_Min_Lim=Roller_Angle-Angle_Tol
#基于角度筛选符合要求的节点
for i in range(len(Inner_Ring_Node_ID_List)):
#计算将节点变换到以第一个滚子为起始位置时的角度(局部坐标)
Local_CSYS_PolarCoord_Theta=\
Inner_Ring_Node_Set_PolarCoord_Theta[i]+2.0*math.pi-Roller1_Pos_Coord_Theta
if Local_CSYS_PolarCoord_Theta>=2.0*math.pi:
Local_CSYS_PolarCoord_Theta=Local_CSYS_PolarCoord_Theta-2.0*math.pi
#筛选获得滚子所在区域的节点
if Roller_Angle_Min_Lim<=Local_CSYS_PolarCoord_Theta<=Roller_Angle_Max_Lim:
Inner_Ring_Roller_Node.append(My_Inner_Ring_Node_Set.nodes[i])
#将节点列表转换为MeshNodeArray对象
Inner_Ring_Roller_Node=mesh.MeshNodeArray(Inner_Ring_Roller_Node)
#创建新的节点集存放外圈滚子的节点
My_Model.rootAssembly.Set(nodes=Inner_Ring_Roller_Node,
name='Inner_Ring_Roller'+str(Roller_Index)+'_Node_Set')
#-------------------------------------------------------------------
#--------------------批量创建各个滚子的所有连接器单元-----------------
for Roller_Index in range(1,Roller_Num+1):
#处理滚子的外圈节点的索引
#读取待处理的节点集
Processed_Set_Name='Outer_Ring_Roller'+str(Roller_Index)+'_Node_Set'
My_Processed_Set=My_Model.rootAssembly.sets[Processed_Set_Name]
#初始化用于存放滚子外圈节点索引排序索引的列表
#获取节点列表
My_Processed_Set_Nodes=My_Processed_Set.nodes
#构建索引和坐标z对应的列表
temp_list=[]
for i in range(len(My_Processed_Set_Nodes)):
z_coord=My_Processed_Set_Nodes[i].coordinates[2]
temp_list.append((z_coord,i))
#基于x坐标进行排序
temp_list.sort()
#提取按x坐标排序后的索引
Outer_Ring_Arrange_Index_List=[item[1] for item in temp_list]
#处理滚子的内圈节点的索引
#读取待处理的节点集
Processed_Set_Name='Inner_Ring_Roller'+str(Roller_Index)+'_Node_Set'
My_Processed_Set=My_Model.rootAssembly.sets[Processed_Set_Name]
#初始化用于存放滚子外圈节点索引排序索引的列表
#获取节点列表
My_Processed_Set_Nodes=My_Processed_Set.nodes
#构建索引和坐标z对应的列表
temp_list=[]
for i in range(len(My_Processed_Set_Nodes)):
z_coord=My_Processed_Set_Nodes[i].coordinates[2]
temp_list.append((z_coord,i))
#基于x坐标进行排序
temp_list.sort()
#提取按x坐标排序后的索引
Inner_Ring_Arrange_Index_List=[item[1] for item in temp_list]
#批量生成当前滚子对应的所有连接器单元
for j in range(len(Inner_Ring_Arrange_Index_List)):
#读取当前处理的滚子对应的外圈和内圈滚动节点
Processed_Set_Name='Outer_Ring_Roller'+str(Roller_Index)+'_Node_Set'
Process_Outer_Node_Set=My_Model.rootAssembly.sets[Processed_Set_Name]
Processed_Set_Name='Inner_Ring_Roller'+str(Roller_Index)+'_Node_Set'
Process_Inner_Node_Set=My_Model.rootAssembly.sets[Processed_Set_Name]
#创建基准坐标系(坐标原点位于外圈节点 x轴指向内圈节点)
#origin和point1的输入必须是一个网格节点对象
My_Datum_CSYS=My_Model.rootAssembly.DatumCsysByThreePoints(
origin=Process_Outer_Node_Set.nodes[Outer_Ring_Arrange_Index_List[j]],
point1=Process_Inner_Node_Set.nodes[Inner_Ring_Arrange_Index_List[j]],
coordSysType=CARTESIAN)
My_Datum_CSYS_ID=My_Model.rootAssembly.datums[My_Datum_CSYS.id]
#创建连接线
My_Wire=My_Model.rootAssembly.WirePolyLine(
points=((Process_Outer_Node_Set.nodes[Outer_Ring_Arrange_Index_List[j]],
Process_Inner_Node_Set.nodes[Inner_Ring_Arrange_Index_List[j]]), ),
mergeType=IMPRINT,
meshable=False)
#更改连接线名称
My_Model.rootAssembly.features.changeKey(fromName=My_Wire.name,
toName='Spring-'+str(Roller_Index)+'-'+str(j+1))
#存放连接线至几何集 合
My_Edges=My_Model.rootAssembly.edges[0:1]
My_Model.rootAssembly.Set(edges=My_Edges,
name='Spring-'+str(Roller_Index)+'-'+str(j+1))
#赋予连接器截面至连接线
My_Region=My_Model.rootAssembly.sets['Spring-'+str(Roller_Index)+'-'+str(j+1)]
My_CSA=My_Model.rootAssembly.SectionAssignment(
sectionName='Bearing_Srping_SEC', region=My_Region)
#赋予连接器方向
My_Model.rootAssembly.ConnectorOrientation(
region=My_CSA.getSet(), localCsys1=My_Datum_CSYS_ID)
#--------------------------------------------------------------------
#-----------------------创建连接器单元作用力的时间历程输出--------------
for Roller_Index in range(1,Roller_Num+1):
for j in range(len(Inner_Ring_Arrange_Index_List)):
#定义时间历程输出的节点集
Hist_Output_Region_Name='Spring-'+str(Roller_Index)+'-'+str(j+1)
Hist_Output_Region=My_Model.rootAssembly.sets[Hist_Output_Region_Name]
#创建各个连接器单元作用力(轴力)时间历程数据输出
My_Model.HistoryOutputRequest(
name=Hist_Output_Region_Name,
createStepName=Step_Name,variables=('CTF1', ),
region=Hist_Output_Region,sectionPoints=DEFAULT, rebar=EXCLUDE)
#--------------------------------------------------------------------
圆柱滚子轴承简化有限元模型连接器载荷提取脚本（Python）
#-* - coding:UTF-8 -*-
#------------------------------------------------------------------
#文件名: Bearing_Load_Extract
#程序用于批量提取轴承简化有限元模型中的轴承载荷(连接器载荷)
#程序与Bearing_Connector_Generator程序配合使用
#作者: Tian Xu, Southwest Jiaotong University
#------------------------------------------------------------------
from abaqus import *
from abaqusConstants import *
#from caeModules import *
from odbAccess import *
from abaqusConstants import *
#--------------------------程序输入参数定义--------------------------
#数据库文件名(.odb)
ODB_File_Name='Main_Bearing.odb'
#读取数据所在分析步名称
Step_Name='Step-1'
#轴承滚子数量
Roller_Num=50
#单个滚子使用的连接器数量
Roller_Spring_Num=14
#------------------------------------------------------------------
#---------------------------批量提取滚子载荷-------------------------
#打开数据库文件(以session对象的形式)
My_Odb=session.openOdb(name=ODB_File_Name)
#打开外部文件用于存放轴承载荷数据
Bearing_Load_File=open('Bearing_Load.txt','w')
#初始化存放轴承载荷数据的列表
Bearing_Load_List=[]
#批量读取连接器载荷
#初始化连接器编号
Connector_Index=0
for Roller_Index in range(1,Roller_Num+1):
for j in range(1,Roller_Spring_Num+1):
#累加连接器编号
Connector_Index=Connector_Index+1
#构造当前连接器的时间历程输出名称
Output_Hist_Name='Connector element total force: CTF1 PI: rootAssembly Element '+\
str(Connector_Index)+' in ELSET SPRING-'+str(Roller_Index)+'-'+str(j)
#构造XYData的显示名称(不影响数据读取)
XYData_Name='Bearing Load-'+str(Roller_Index)+'-'+str(j)
#通过XYData的形式读取时间历程输出数据
My_Hist_Data=session.XYDataFromHistory(
name=XYData_Name,odb=My_Odb,
outputVariableName=Output_Hist_Name,steps=(Step_Name,))
#储存连接器载荷数据(CTF1)(分析步最后时刻的载荷)
Bearing_Load_List.append(My_Hist_Data[-1][1])
#写入连接器载荷数据至文件
Bearing_Load_File.write('%6d %6d %.6f\n' %(Roller_Index,j,My_Hist_Data[-1][1]))
#关闭文件
Bearing_Load_File.close()
#-------------------------------------------------------------------
圆柱滚子轴承非线性载荷-变形关系计算程序（MATLAB）
clear,clc
%--------------------------------------------------------------------------
%程序用于计算圆柱滚子轴承中滚子-滚道接触的载荷-变形关系
%--------------------------------------------------------------------------
%------------------------------参数输入-------------------------------------
L=113; %滚子有效长度(mm)
D=63; %滚子直径(mm)
Play_Val=0; %轴承游隙(mm)
%载荷-变形系数计算使用的常数(Brändlein et al. (1999))
C=2.65e4;
n1=0.9189;
n2=0;
n3=1.0811;
%--------------------------------------------------------------------------
%----------------------------计算载荷-变形关系------------------------------
%设置压缩变形量(mm)(零游隙)
Delta_Negative=linspace(10,0,21)';
%计算压缩变形量对应的滚子载荷(N)
Q_Negative=C*(L.^n1).*(D.^n2).*(Delta_Negative.^n3);
Delta_Negative=-1*Delta_Negative;
Q_Negative=-1*Q_Negative;
%设置拉伸变形量(mm)(零游隙)
Delta_Positive=10;
%计算拉伸变形量对应的滚子载荷(N)
Q_Positive=0;
%组装滚子的完整载荷-变形关系(零游隙)
Delta=[Delta_Negative;Delta_Positive];
Q=[Q_Negative;Q_Positive];
%计算滚子的完整载荷-变形关系(考虑游隙)
Delta=Delta-Play_Val;
%输出滚子的完整载荷-变形关系(考虑游隙)
Roller_Load_Disp_Data=[Q Delta];
%绘制滚子的载荷-变形关系
plot(Delta,Q,'-o');
hold on;scatter(0,0);
%--------------------------------------------------------------------------
轴承载荷分布绘制（MATLAB）
clear,clc
%--------------------------------------------------------------------------
%程序用于读取轴承载荷数据并绘制轴承载荷分布
%轴承载荷数据文件格式:滚子编号-单个滚子的连接器编号-连接器载荷
%--------------------------------------------------------------------------
%------------------------------参数输入-------------------------------------
%轴承滚子数量
Roller_Num=50;
%--------------------------------------------------------------------------
%------------------------读取并处理轴承载荷数据------------------------------
%读取载荷数据文件
Bearing_Load_Data=importdata('Bearing_Load.txt');
%获取单个轴承滚子使用的连接器数量
Roller_Spring_Num=size(Bearing_Load_Data,1)/Roller_Num;
%初始化数组存放轴承滚子全局载荷(单个滚子所有连接器载荷之和)
Roller_Global_Load=zeros(Roller_Num,1);
%循环读取轴承载荷数据
Data_Pair_Num=0;
for Roller_Index=1:Roller_Num
k=0;
temp_Load_Data=zeros(Roller_Spring_Num,1);
for Spring_Index=1:Roller_Spring_Num
k=k+1;
Data_Pair_Num=Data_Pair_Num+1;
temp_Load_Data(k,1)=Bearing_Load_Data(Data_Pair_Num,3);
end
Roller_Global_Load(Roller_Index)=sum(temp_Load_Data);
field_name = ['Roller', num2str(Roller_Index)];
Bearing_Data_Struct.(field_name)=temp_Load_Data;
end
%绘制滚子全局载荷分布
Theta=linspace(0,2*pi,Roller_Num)';
plot(Theta,Roller_Global_Load,'-o');
%--------------------------------------------------------------------------
