一种基于SPH的筒形件旋压成形过程分析方法

一种基于sph的筒形件旋压成形过程分析方法
技术领域
1.本发明涉及计算力学领域，尤其涉及一种基于sph的筒形件旋压成形过程分析方法。


背景技术：

2.传统的以网格为主的数值方法（如有限元法、有限差分法等）虽然可以对筒形件旋压成形过程进行数值模拟，但当面临大变形及裂纹扩展时却容易产生网格畸变而造成计算误差甚至计算中止，且在建立多旋轮动态运动分析对称模型中存在困难。因此，基于sph构建一种新型的筒形件旋压成形过程力学分析数值方法，能够为旋压过程中面临的复杂动态力学分析提供科学依据，同时为该领域的力学分析提供新的计算方法。
3.光滑粒子流体动力学方法（smoothed particle hydrodynamics，简称sph）是近年来逐渐发展起来的一种纯拉格朗日无网格计算方法，主要用于天体物理学、爆炸力学、流体力学和固体力学等领域，然而在面临金属成型系统领域的研究还处于起步阶段。本发明将sph方法引入筒形件旋压成形过程力学分析领域，基本思想是将毛坯、芯模、旋轮部件离散为有限个无网格连接的离散点，各部件的物理、力学及运动参数，如密度、质量、体积、位置、速度、加速度、应力等都由离散点携带，并通过光滑核函数及质量守恒、动量守恒和本构方程建立离散点间的联系。由于离散点间无网格连接，所以不会产生网格畸变，特别适合求解金属成型的大变形动态力学行为。同时，为使sph方法满足旋压过程工艺特点、减小计算成本，一方面，本发明将毛坯视为变形体，将芯模和旋轮视为刚性体，并将毛坯的旋转等效到旋轮上，即使毛坯和芯模不旋转，而旋轮绕毛坯和芯模做旋转及纵向的进给运动；另一方面，建立了适用于sph在旋压成型过程分析中应用的轴对称分析方法，最终实现了该高性能无网格算法在该领域内的成功应用。


技术实现要素：

4.本发明提供了一种基于sph的筒形件旋压成形过程分析方法，其特征在于,包括如下步骤：步骤一：根据待研究对象的几何及运动信息确定分析模型的对称性，并根据毛坯的几何信息将毛坯离散为三维实粒子模型，将毛坯实粒子定义为ⅰ类型粒子并进行初始状态的密度、质量、弹性模量、泊松比、位置、速度、应力、光滑长度的参数赋值；步骤二：设置计算旋压成型力学分析过程中需要采用的数值处理技术及参数，并设置时间步长及时间步数；步骤三：根据分析模型对称性建立毛坯对称边界虚粒子模型，根据芯模、旋轮的几何信息将芯模、旋轮离散为三维虚粒子模型，并将毛坯边界虚粒子、芯模虚粒子、旋轮虚粒子分别定义为ⅱ、ⅲ、ⅳ类型虚粒子并进行密度、质量、位置、速度、应力、光滑长度的参数赋值；步骤四：分别进行
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型、
ⅰ‑ⅲ
类型、
ⅰ‑ⅳ
类型粒子的相邻粒子对的搜
索；步骤五：通过
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型粒子对进行密度求解，得到ⅰ类型粒子密度变化率；步骤六：通过
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型粒子进行应力应变的求解，得到ⅰ类型粒子应力变化率，并由应力应变结果得到由内力引起的ⅰ类型粒子内力速度变化率；步骤七：对
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型、
ⅰ‑ⅲ
类型粒子对进行人工粘度求解，得到ⅰ类型粒子的人工粘度速度变化率；步骤八：通过
ⅰ‑ⅲ
、
ⅰ‑ⅳ
类型粒子对进行外力作用下引起的外力速度变化率的求解，得到ⅰ类型粒子的外力速度变化率；步骤九：根据ⅰ类型粒子的密度变化率、应力变化率、内力速度变化率、人工粘度速度变化率、外力速度变化率之和更新ⅰ类型粒子的质点密度、质量、应力、位置、速度信息；步骤十：判定时间步是否达到预设的时间步，如果未达到则进入下一个循环，即重复步骤三到步骤八，如果达到则计算终止。
5.进一步的，步骤一中，为减少计算量和提高计算效率，会根据待研究对象的几何及运动信息建立不同的对称分析模型，对双旋轮且旋轮间夹角为180
°
建立1/2对称模型，对三旋轮且旋轮间夹角为120
°
建立1/3对称模型，对单旋轮建立全分析模型。
6.进一步的，步骤二中，所述数值处理技术包括密度计算方法、光滑核函数、粒子搜索方法、刚体模拟方法、时间积分方法、质量放大方法、对称模型方法、人工粘度方法。
7.进一步的，步骤三中，ⅱ类型粒子仅在采用对称分析模型中才存在，其位置处在毛坯对称边界处，且在周向设置的ⅱ类型粒子数量列数不应该少于ⅰ类型粒子的光滑半径与圆周向初始间距的倍数，即：其中， nⅱ代表ⅱ类型粒子在圆周方向的列数，h为光滑长度， r为与光滑核函数相关的参数，h
×
r为光滑半径，r0代表ⅰ类型粒子在圆周向的初始间距。且ⅱ类型粒子的密度、质量、光滑长度、应力值、位置、速度与对称位置的ⅰ类型粒子的数值对称。
8.进一步的，步骤一及步骤三进行位置和速度参数设置时，将毛坯与芯模的旋转过程等效到旋轮上来，即采用芯模和毛坯固定不动，旋轮做进给运动同时绕毛坯旋转的方法进行模拟。同时，当旋轮旋转至对称边界时，其超过对称边界的部分会自动切换至另一端对称边界，并根据模型对称特性及平面旋转矩阵对称其运动参数。
9.进一步的，步骤六中，所述毛坯采用分段弹性模型模拟弹塑性过程来进行应力应变的求解。
10.本发明的有益效果为：与目前常用的数值模拟方法，如有限元法、有限差分法等方法相比，本发明提供的一种基于sph的筒形件旋压成形过程分析方法是一种纯拉格朗日无网格方法，因此不会在产生大变形时发生网格畸变。同时，通过提供的对称模型建立方法并设置质量放大系数，可以显著的减少计算量，为筒形件旋压成形过程的动态力学分析提供科学依据，同时为该领域的力学分析提供了一种新的计算方法，本发明可以同时为其他旋压成型问题提供有益的
帮助。
附图说明
11.图1为本发明所述方法计算流程图。
12.图2为当旋轮移动至对称边界时，某旋压成型过程的1/2模型及各类型粒子示意图。
13.图3为采用本发明方法对某旋压成型过程进行分析计算时，计算时间为1s时模型的位置及其mises应力分布。
14.图4为采用本发明方法对某旋压成型过程进行分析计算时，计算时间为2s时模型的位置及其mises应力分布。
15.图5为采用本发明方法对某旋压成型过程进行分析计算时，计算时间为3s时模型的位置及其mises应力分布。
16.图6为采用本发明方法对某旋压成型过程进行分析计算时，计算时间为4s时模型的位置及其mises应力分布。
17.图7为采用有限元计算为某旋压成型过程进行分析计算时，计算时间为1s时模型的位置及其mises应力分布。
18.图8为采用有限元计算为某旋压成型过程进行分析计算时，计算时间为2s时模型的位置及其mises应力分布。
19.图9为采用有限元计算为某旋压成型过程进行分析计算时，计算时间为3s时模型的位置及其mises应力分布。
20.图10为采用有限元计算为某旋压成型过程进行分析计算时，计算时间为4s时模型的位置及其mises应力分布。
21.图11为某旋压加工模型塑性阶段应力应变数据图。
具体实施方式
22.为使本发明实施方式的目的、技术方案和优点更加清楚，下面将结合本发明实施方式中的实施案例和附图，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。因此，以下对在附图中提供的本发明的实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施方式。基于本发明中的实施方式，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
23.如图1，一种基于sph的筒形件旋压成形过程分析方法，包括步骤如下：步骤一：根据待研究对象的几何及运动信息确定分析模型的对称性，并根据毛坯的几何信息将毛坯离散为三维实粒子模型，将毛坯实粒子定义为ⅰ类型粒子并进行初始状态的密度、质量、弹性模量、泊松比、位置、速度、应力、光滑长度的参数赋值；步骤二：设置计算旋压成型力学分析过程中需要采用的数值处理技术及参数，并设置时间步长及时间步数；
步骤三：根据分析模型对称性建立毛坯对称边界虚粒子模型，根据芯模、旋轮的几何信息将芯模、旋轮离散为三维虚粒子模型，并将毛坯边界虚粒子、芯模虚粒子、旋轮虚粒子分别定义为ⅱ、ⅲ、ⅳ类型虚粒子并进行密度、质量、位置、速度、应力、光滑长度的参数赋值；步骤四：分别进行
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型、
ⅰ‑ⅲ
类型、
ⅰ‑ⅳ
类型粒子的相邻粒子对的搜索；步骤五：通过
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型粒子对进行密度求解，得到ⅰ类型粒子密度变化率；步骤六：通过
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型粒子进行应力应变的求解，得到ⅰ类型粒子应力变化率，并由应力应变结果得到由内力引起的ⅰ类型粒子内力速度变化率；步骤七：对
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型、
ⅰ‑ⅲ
类型粒子对进行人工粘度求解，得到ⅰ类型粒子的人工粘度速度变化率；步骤八：通过
ⅰ‑ⅲ
、
ⅰ‑ⅳ
类型粒子对进行外力作用下引起的外力速度变化率的求解，得到ⅰ类型粒子的外力速度变化率；步骤九：根据ⅰ类型粒子的密度变化率、应力变化率、内力速度变化率、人工粘度速度变化率、外力速度变化率之和更新ⅰ类型粒子的质点密度、质量、应力、位置、速度信息；步骤十：判定时间步是否达到预设的时间步，如果未达到则进入下一个循环，即重复步骤三到步骤八，如果达到则计算终止。
24.进一步的，步骤一中，为减少计算量和提高计算效率，会根据待研究对象的几何及运动信息建立不同的对称分析模型，对双旋轮且旋轮间夹角为180
°
建立1/2对称模型，对三旋轮且旋轮间夹角为120
°
建立1/3对称模型，对单旋轮建立全分析模型。
25.进一步的，步骤二中，所述数值处理技术包括密度计算方法、光滑核函数、粒子搜索方法、刚体模拟方法、时间积分方法、质量放大方法、对称模型方法、人工粘度方法、平均速度方法。
26.进一步的，步骤三中，ⅱ类型粒子仅在采用对称分析模型中才存在，其位置处在毛坯对称边界处，且在周向设置的ⅱ类型粒子数量列数不应该少于ⅰ类型粒子的光滑半径与圆周向初始间距的倍数，即：其中， nⅱ代表ⅱ类型粒子在圆周方向的列数，h为光滑长度， r为与光滑核函数相关的参数，h
×
r为光滑半径，r0代表ⅰ类型粒子在圆周向的初始间距。且ⅱ类型粒子的密度、质量、光滑长度、应力值、位置、速度与对称位置的ⅰ类型粒子的数值对称。
27.进一步的，步骤一及步骤三进行位置和速度参数设置时，将毛坯与芯模的旋转过程等效到旋轮上来，即采用芯模和毛坯固定不动，旋轮做进给运动同时绕毛坯旋转的方法进行模拟。同时，当旋轮旋转至对称边界时，其超过对称边界的部分会自动切换至另一端对称边界，并根据模型对称特性及平面旋转矩阵对称其运动参数，如图2所示。
28.进一步的，步骤六中，所述毛坯采用分段弹性模型模拟弹塑性过程来进行应力应变的求解。
29.下面以具体数值例子对本发明方法进行说明。
30.实施例1：以某旋压加工模型为例，采用本发明的数值仿真方法进行仿真模拟。该旋压模型由两个旋轮组成，两个旋轮间夹角为180
°
，旋压过程采用反旋转。其中，旋轮入旋角为20
°
，退出角为30
°
，圆角半径为5mm，毛坯厚度为10mm，内径为120mm，长度为30mm，芯模直径120mm，长度32mm，减薄率为23%，进给速度为0.01m/s，主轴转速为50rad/s，毛坯密度为7800kg/m3，弹性模量为20gpa，泊松比为0.3，旋轮塑性阶段应力应变数据表如图11（图11为某旋压加工模型塑性阶段应力应变数据图）所示。具体实施过程如下：步骤一：将毛坯模型离散为有限个实粒子。由于该模型为轴对称结构，为节约计算时间，因此建立1/2模型对模型进行分析。其中毛坯实粒子378个，粒子初始间轴向间距为5mm，厚度方向为10/3mm，圆周向方向为每圈40个粒子，并将弹性模量、泊松比、密度、质量等物理参数分别给粒子进行赋值，光滑长度选择为毛坯圆周外径粒子初始平均间距的1.2倍，即0.0066m。
31.步骤二：设置时间步长为4
×
108，共计算106个时间步，选择核函数为三次样条核函数，选择密度计算方法为连续密度法、粒子搜索方法为全配对搜索法、刚体模拟方法采用固壁边界处理、时间积分方法为跳蛙法，质量放大系数选择为104，采用轴对称方法进行模拟。
32.步骤三：根据芯模、旋轮的几何信息建立芯模和旋轮的三维虚粒子模型，其中芯模虚粒子共780个，为简化模型，仅对与毛坯相接触位置附近的旋轮粒子进行建模，共782个，虚粒子间距为毛坯实粒子的1/2。芯模位置固定，速度为0，而旋轮则根据主轴转速其位置和速度在每个时间步进行更新。由于是轴对称模型且核函数选择为三次样条核函数，则毛坯虚粒子数量设置为2排，共84个，位置在对称边界附近与实粒子相邻，且其位置、速度与应力毛坯实粒子轴对称，其弹性模量、密度、质量、光滑长度、应力与毛坯实粒子相等，并将毛坯虚粒子、芯模、旋轮分别定义为ⅱ、ⅲ、ⅳ类型虚粒子。离散模型的初始状态如图3所示。
33.步骤四：进行粒子对的搜索，分别得到每个时间步
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型、
ⅰ‑ⅲ
类型、
ⅰ‑ⅳ
类型粒子对的个数和对应粒子。
34.步骤五：通过
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型粒子对进行密度求解，得到ⅰ类型粒子密度变化率。
35.步骤六：通过
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型粒子进行应力应变的求解，得到ⅰ类型粒子应力变化率，并由应力应变结果得到由内力引起的ⅰ类型粒子内力速度变化率。
36.步骤七：对
ⅰ‑ⅰ
类型、
ⅰ‑ⅱ
类型、
ⅰ‑ⅲ
类型粒子对进行人工粘度求解，得到ⅰ类型粒子的人工粘度速度变化率；步骤八：通过
ⅰ‑ⅲ
、
ⅰ‑ⅳ
类型粒子对进行外力作用下引起的外力速度变化率的求解，得到ⅰ类型粒子的外力速度变化率。
37.步骤九：根据ⅰ类型粒子的密度变化率、应力变化率、内力速度变化率、人工粘度速度变化率、外力速度变化率之和更新ⅰ类型粒子的质点密度、质量、应力、位置、速度信息。
38.步骤十：判定时间步是否达到预设的时间步，如果未达到则进入下一个循环，即重复步骤三到步骤八，如果达到则计算终止。
39.通过以上步骤，计算得到了该筒形件旋压成型的运动及力学过程。其在1s、2s、3s、4s时的位置及mises应力状态如图3-图6所示。为观察筒形件旋压后的成型质量，采用绝对
百分比误差（mape）及均方根误差（rmse）两种指标对预测模型进行精度评价,其中，mape表达了计算值与目标值的平均偏离程度，其越接近0表示效果越好；rmse为计算值与目标值之间的偏差，越接近0表示与目标值越吻合。
40.经计算，该次旋压结果的mape值为1.0075%，rmse值为0.012043。mape值和rmse值可以来评价该旋压条件下的旋压效果，而如果要得到更好的旋压效果或提高旋压效率，则可在调整不同的旋压参数后，通过数值计算对比其计算结果，从而对旋压参数进行优化。
41.同时，为验证本发明方法的正确性，通过有限元方法对该例进行了计算，其在1s、2s、3s、4s时的位置及mises应力状态如图7-图10所示，可见采用本发明的sph方法与采用有限元法在进行该旋压成型模型计算时具有较高的一致性，证明了本发明算法可以较好的模拟筒形件旋压成型的力学过程。