CN114357717A - 一种应用于双色币压印成形仿真中的基于改进接触算法的物质点法 - Google Patents

Abstract

本发明提供了一种应用于双色币压印成形仿真中的基于改进接触算法的物质点法，该方法包括以下步骤：对双色币压印成形工艺过程建立动力学物理模型；对压印成形过程动力学模型进行物质点积分；对成形模具进行有限元三角形/四边形单元离散；根据坯饼数据模型建立背景网格，并对双色币坯饼进行四面体/六面体单元离散、建立铜、铝材料区域的质点集；经过计算质点的应变增量和旋率增量，根据本构模型求解质点的应力及弹塑性相关变量，更新所有质点的密度。本发明通过物质点法来避免双色币压印成形有限元模拟中材料自接触判断问题。

Description

一种应用于双色币压印成形仿真中的基于改进接触算法的物 质点法

技术领域

本发明涉及压印技术领域，具体涉及一种应用于双色币压印成形仿真中的基于改进接触算法的物质点法

背景技术

双色金属币制作是当今硬币防伪的一个很重要的技术，有利于提高居民们的使用金属币的防伪性能。同时，设计精美的双色币观赏性强，极具收藏价值。随着科学技术不断革新，对金属币的制造也提出了更高的要求，特别是对于制作这种拥有明显优势并且技术含量很高的双色币。双色币在实际压印过程中，两种不同材料之间相互接触挤压，容易产生应力集中问题，进而导致微裂纹甚至裂纹的产生，严重影响产品质量。在实际生产过程中，需要不断的试模才能克服应力集中问题。这会极大增加生产周期和生产成本。采用仿真算法对该压印过程进行分析，可以辅助工程师为解决应力集中问题提供理论依据，同时有效地降低生产周期和生产成本。有限元法在很多工程领域得到了广泛的应用和认可，但是在压印成形大变形问题中同样存在网格畸变问题，导致计算终止或者仿真结果偏离实际。同时，针对双色币材料的压印仿真分析，有限元法需要额外判断不同材料界面的接触穿透问题，极大地增加了仿真时间和难度。而物质点法是一种采用拉格朗日质点的欧拉网格双重描述的、非常适合于分析涉及大变形和接触问题的数值计算方法。它弥补了拉格朗日法在网格划分问题中发生网格畸变导致参元插值与积分精度下降、结果不收敛的问题，在接触和冲击等涉及大变形问题上具有明显优势。

基于物质点法在大变形数值模拟领域中的明显优势，本研究发明拟采用物质点法对双色币压印成形进行仿真分析。但是，传统的物质点法存在虚接触的缺陷，导致动量不守恒和界面缺陷。本发明提出改进接触算法的物质点法，避免了虚拟接触，提高了物质点法接触判断的准确性。此外，通过物质点法来避免双色币压印成形有限元模拟中材料自接触判断问题，极大降低了仿真中接触算法难度和计算量。

发明内容

本发明的目的是针对物质点法中点对点接触会造成虚接触的问题，而提供了一种应用于双色币压印成形仿真中的基于改进的接触算法的物质点法。本发明的目的是这样实现的：包括如下步骤：

步骤1：对双色币压印成形工艺过程建立动力学物理模型；

步骤2：对压印成形过程动力学模型进行物质点积分；

步骤3：对成形模具进行有限元三角形/四边形单元离散；

步骤4：根据坯饼数据模型建立背景网格，并对双色币坯饼进行四面体/六面体单元离散、建立铜、铝材料区域的质点集；

步骤5：计算背景网格节点载荷；

步骤6：通过插值形函数将质点携带的物理量映射到背景网格节点；

步骤7：根据双色币不同材料特性求解系统自适应时间步长；

步骤8：求解动量方程并施加边界条件；

步骤9：更新质点的物理量；

步骤10：更新背景网格点的物理量；

步骤11：更新模具空间位置；

步骤12：判断物质点的穿透状态及求解穿透距离；

步骤13：求解接触力；

步骤14：再次更新质点的物理量；

步骤15：施加边界条件；

步骤16：计算质点的应变增量和旋率增量，根据本构模型求解质点的应力及弹塑性相关变量，更新所有物质点的密度。

本发明还包括这样的一些特征：

1.步骤1具体为：

从力学角度出发，压印成形过程是一个拟静态的低速运动过程。一般采用隐式积分方法求解。但是，考虑到纪念币表面图案复杂，需要数以百万计的网格来描述，导致线性方程组求解会占用大量内存。同时，压印过程伴随复杂的材料非线性。采用牛顿-拉夫森迭代求解可能会导致不收敛的发生。因此，本研究拟采用动力显式算法来描述及求解压印成形过程。设坯饼为被研究物体，其发生塑性变形时应满足基本平衡方程：

其中，σ_ij分别为柯西应力张量，b_i为体积力，ρ为材料密度，γ为阻尼系数，

和

分别是物体内任一点的速度和加速度，i和j分别为坐标方向。根据虚位移原理、相关应力边界条件和高斯散度定理，可推导系统的虚功率方程为：

式中，V为物体区域，

是虚速度，

是对应于柯西应力σ_ij的虚应变速率。其中S_p为已知外力p_i的表面；S_c为与另一物体接触的表面，记接触表面力为q_i；惯性力和阻尼力功率项反映了物体系统的惯性效应和物理阻尼效应。

2.步骤2具体为：

物质点法中，每一个质点占据一定的体积空间和携带一定质量以及其他相关物理量。任意空间点x的材料密度可以表示为：

其中n_p为质点数，m_p为质点质量，δ()为Dirac函数，x_p为质点p的坐标。背景网格中任意质点p的位移u_ip可通过背景网格结点形函数插值获取：

其中，N_Ip为质点p相关的结点I的形函数,u_iI为背景网格结点I的第i方向的速度。将上述密度和位移公式带入步骤1中的系统虚功率方程，同时考虑压印成形过程中无体积力b_i和外力p_i，可得：

其中，

为结点I质量，

和

分别为结点I第i个方向的速度和加速度。

为结点I第i个方向的外力，

为结点I第i个方向的内力。

3.步骤4具体为：

首先，对坯饼所占区域划分规则六面体背景网格(沿x,y,z方向的网格数分别为Gx,Gy,Gz)，并求出每个背景网格的积分点。然后，对铜、铝两种不同材料区域分别进行四面体网格划分。如图1所示，为简单的二维示例(背景网格为矩形，坯料网格为三角形。。实际三维坯饼尺寸如图2所示)。每一个小矩形表示背景网格单元，实心圆代表网格节点，小实心三角形代表铜物质点，小实心正方形代表铝物质点。接着，在铜材料区域找出每一个四面体单元的背景网格单元集合I_Cu。求该集合的具体方法如下：将背景网格单元沿x,y,z方向从1开始编号至G＝Gx*Gy*Gz。根据四面体单元的四个节点坐标，可以找出一个空间长方体来刚好包围该四面体。假设背景网格最左前点坐标(x_min，y_min，z_min)和最右后点坐标(x_max，y_max，z_max)。依据空间任意一点x_p＝(x，y，z)的背景网格编号

找到该长方体所包含的所有背景网格集合I_Cu。将该集合中所有单元的高斯点组成的集合定义为点集II_Cu。根据铜材料的表面边界单元(由网格剖分系统自动生成并编写程序读取该表面单元)，将点集II_Cu中铜材料外的高斯点剔除，剩下的点形成的集合则为铜材料域的质点集合III_cu，如图1实心小三角形所示。同理，对铝材料区域进行类似的步骤，获得铝材料的质点集合III_Al，如图1实心小正方形所示。最终，将两集合III_Cu和III_Al合并为坯饼的物质点集III，并为每个物质点赋予相对应的材料属性，如密度，体积，泊松比，屈服应力、硬化指数、强化系数等参数。为了准确描述硬币初始轮廓,将坯饼表面各三角形单元的高斯点集添加到坯饼质点集合中，组成新的坯饼质点集。如图3所示为其四分之一模型被切除，以方便查看内部质点。将坯饼表面质点集分为三部分：上表面、侧面和下表面质点集，它们分别与上模，中圈和下模构成接触对。在仿真过程中，此三对接触对分别进行接触判断处理。

4.步骤7具体为：

在通常情况下，针对单一材料压印成形仿真问题，我们采用中心差分动力显式算法进行求解。该算法为条件稳定，系统时间步长Δt必须小于临界步长Δt_cr，即

Δt＝αΔt_cr

式子中，α是一个常数，取值范围为：0.8≤α≤0.98

由于CFL(Courant Friedrichs Lewy，柯朗-弗里德里希斯-列维)条件，为了保证显式时间积分的稳定性，采用有限元法进行动态显式分析时，一个时间步长内波传输的距离必须不超过一个单元。而弹性材料的声速为：

式子中，E为弹性模量，v为泊松比，ρ为材料密度。假定单元的名义长度为l_e，则临界时间步长为：

随变形的增加，单元名义尺寸逐渐减小，临界时间步长也逐渐减小。当网格发生严重畸变时，所求解的最小名义长度趋近于0，导致时间步长亦趋近于0，从而仿真分析无法进行。

物质点法与有限元法不同，所有时间积分都是需要在固定的背景网格上进行。在确定临界时间步长时，不仅需要考虑材料的声速，同时需要考虑质点的速度。特别地，在超高速撞击问题中，质点速度较大，与材料声速在同一数量级甚至超过材料声速。因此，质点速度对仿真计算的影响不可忽略。在压印成形模拟过程中，基于动力显式中心差分算法的物质点法求解，在每一步需要计算临界时间步长使得数值求解稳定：

式子中，c为质点的声速，u为质点速度。在使用均匀背景网格时，l_g又可以取为d_c，即背景网格尺寸。本发明中坯饼有两种不同材料，即内铜外铝。时间步长的大小由系统的属性和模型材料参数所决定。由于涉及两种不同材料，需要分别计算两种物质区域内的时间步长并取其较小值。铜质材料中时间步长为：

式子中，

和

分别为质点集III_Cu中声速和任意质点p₁质点速度。

同理，铝质材料中时间步长为：

式子中，

和

分别为质点III_Al中声速和任意质点p₂的质点速度。

我们需要计算出两种材料的波速，以及压印过程中两种材料物质点的速度，进而比较两者的临界时间步长，最终确定系统的临界时间步长为两者最小值。因此，系统的临界时间步长为：

式子中，

为Cu的时间步长，

为Al的时间步长。

5.步骤12为物质点穿透状态判断及求解穿透距离，主要包含以下内容：如图4所示，根据传统物质点法理论，模具质点j′与坯饼质点j对背景网格结点I的速度均有贡献(因为它们各自所在的背景网格共享结点I)，因而判定模具与该坯饼质点接触。而实际上它们并没有接触，从而产生虚假接触。本发明采用点面接触算法替代传统物质点的点-点接触算法以避免虚拟接触，提高仿真计算准确性。判断坯饼表面每个质点与上模、下模和中圈之间的接触状态和计算接触值(接触力、穿透深度和摩擦力)时，在每一个时间步长里都会造成很大的计算量。为了减少接触判断的计算时间，采用了全局搜索与局部搜索相结合的策略。全局搜索是为每个质点找到与之相接触的候选模具单元，过程如下所示。

构建包含特定模具网格的块：块的大小根据各模具的最大和最小空间坐标来定义。块的左下角和右下角的坐标分别为(x_min，y_min，z_min)和(x_max，y_max，z_max)。然后将该块区域划分为子单元格，每个方向上的单元格数目分别为N_x，N_y，N_z，并分别为每个单元格沿x,y,z三个方向进行编号。注意这个块随模具网格一起运动。为了区分这种单元格和背景单元网格，我们将这种单元格命名为空间单元格。

搜索每个单元格所包含的候选的模具网格单元。如图5所示为某二维物体的单元格示意图。共有N＝N_x×N_y＝8×5＝40个空间单元格，图中只标注了四个角的单元格编号。编号为21的单元格包含五个模具网格单元，其节点用实心圆表示。查找每个单元格所包含的候选单元一般包括两个步骤。首先，将某模具空间格域沿三个方向扩展，扩展的距离为子空间格尺寸的一半，并循环每个模具单元以获取包含它的单元格，最终得到每个单元格包含的候选模具单元。对于给定的粒子x_p＝(x，y，z)，其定位的单元编号为

其中[]是一个四舍五入的运算操作。针对质点x_p，可以定位到包含该质点的子空间格。根据上述的全局搜索结果可以得到该质点的候选模具单元集合。接下来在该候选集中进行局部搜索并进行接触判断。为了判断网格之间的接触状态和质点x_p，我们只需要找到质点是否接触候选单元或者穿透它。执行以下步骤可以实现此目标。

求潜在接触单元：如图6所示，质点x_p，S₁，S₂，S₃，S₄有四个候选单元。找到离x_p最近的节点m_s，它与x_p不重合。如果满足下列不等式，则质点将与S_i单元接触

(C_i×g₂)·(C_i×C_i+1)>0，(C_i×g₂)·(g₂×C_i+1)>0

这里的i＝1,2,3,4为待判断的第i个模具单元。C_i和C_i+1是两个单位向量，它们有共同的起始点，沿着S_i单元的两条相邻边，g₁从m_s开始，到x_p结束。g₂是在单元S_i在向量g₁上的投影，写成：

g₂＝g₁-(g₁·m)m

这里

如图7所示，为了计算x_c的空间坐标，在单元上建立参数公式：

r(ξ，η)＝f₁(ξ，η)i₁+f₂(ξ，η)i₂+f₃(ξ，η)i₃

其中

N_j(ξ，η)是壳体单元的双线性函数，

为节点j的第i个坐标，(ξ，η)为单元内任意一点的局部坐标，i₁，i₂，i₃分别为空间坐标下的单位向量。t是从原始坐标开始，到x_p结束的向量。因为向量

是垂直于接触单元，接触点x的局部坐标(ξ_c，η_c)必须满足：

判断接触状态和穿透深度d的计算：

其中n_c为在点(ξ_c，η_c)接触单元的单位向外的法向量，计算方法为：

附图说明

图1是物质点离散的示意图

图2是双色币截面示意图

图3是质点潜在接触对的示意图

图4是物质点虚接触示意图

图5是背景网格单元区域划分的示意图

图6是质点间的接触状态判断示意图

图7是计算质点接触位置示意图

图8是计算质点接触力的示意图

具体实施方式

为了阐明本发明的目的、技术方案和优点，现结合所附结构和实施例对本发明进行更详细的说明。本发明是一种应用于双色币压印成形仿真的基于改进的接触算法的物质点法，包括以下步骤：

步骤1：对双色币压印成形工艺过程建立动力学物理模型；

从力学角度出发，压印成形过程是一个拟静态的低速运动过程。一般采用隐式积分方法求解。但是，考虑到纪念币表面图案复杂，需要数以百万计的网格来描述，导致线性方程组求解会占用大量内存。同时，压印过程伴随复杂的材料非线性。采用牛顿-拉夫森迭代求解可能会导致不收敛的发生。因此，本研究拟采用动力显式算法来描述及求解压印成形过程。设坯饼为被研究物体，其发生塑性变形时应满足基本平衡方程:

其中，σ_ij，j分别为柯西应力张量，b_i为体积力，ρ为材料密度，γ为阻尼系数，

和

分别是物体内任一点的速度和加速度，i和j分别为坐标方向。根据虚位移原理、相关应力边界条件和高斯散度定理，可推导系统的虚功率方程为：

式中，V为物体区域，

是虚速度，

是对应于柯西应力σ_ij的虚应变速率。

其中S_p为已知外力p_i的表面；S_c为与另一物体接触的表面，记接触表面力为q_i；惯性力和阻尼力功率项反映了物体系统的惯性效应和物理阻尼效应。

步骤2：对压印成形过程动力学模型进行物质点积分；

物质点法中，每一个质点占据一定的体积空间和携带一定质量以及其他相关物理量。任意空间点x的材料密度可以表示为：

其中n_p为质点数，m_p为质点质量，δ()为Dirac函数，x_p为质点p的坐标。背景网格中任意质点p的位移u_ip可通过背景网格结点形函数插值获取：

其中，N_Ip为质点p相关的结点I的形函数,u_iI为背景网格结点I的第i方向的速度。将上述密度和位移公式带入步骤1中的系统虚功率方程，同时考虑压印成形过程中无体积力b_i和外力p_i，可得：

其中，

为结点I质量，

和

分别为结点I第i个方向的速度和加速度。

为结点I第i个方向的外力，

为结点I第i个方向的内力。

步骤3：对成形模具进行有限元三角形/四边形单元离散；

使用开源软件对双色币进行四面体网格划分并获取表面三边形网格离散，建立只包含双色币材料区域的背景网格。

步骤4：根据坯饼数据模型建立背景网格，并对双色币坯饼进行四面体/六面体单元离散、建立铜、铝材料区域的质点集；

坯饼进行质点离散后，其质点密度可以表示为：

其中n_p为质点数，m_p为质点质量，δ()为Dirac函数，x_p为质点p的坐标。文中物质点法(Material Point Method,MPM)将材料域即初始双色币坯料模型离散为一组粒子，如图1所示。首先，对坯饼所占区域划分规则六面体背景网格(沿x,y,z方向的网格数分别为Gx,Gy,Gz)，并求出每个背景网格的积分点。然后，对铜、铝两种不同材料区域分别进行四面体网格划分。如图1所示，为简单的二维示例(背景网格为矩形，坯料网格为三角形。实际三维坯饼尺寸如图2所示)。每一个小矩形表示背景网格单元，实心圆代表网格节点，小实心三角形代表铜物质点，小实心正方形代表铝物质点。接着，在铜材料区域找出每一个四面体单元的背景网格单元集合I_Cu。求该集合的具体方法如下：将背景网格单元沿x,y,z方向从1开始编号至G＝Gx*Gy*Gz。根据四面体单元的四个节点坐标，可以找出一个空间长方体来刚好包围该四面体。假设背景网格最左前点坐标(x_min，y_min，z_min)和最右后点坐标(x_max，y_max，z_max)。依据空间任意一点x_p＝(x，y，z)的背景网格编号

找到该长方体所包含的所有背景网格集合I_Cu。将该集合中所有单元的高斯点组成的集合定义为点集II_Cu。根据铜材料的表面边界单元(由网格剖分系统自动生成并编写程序读取该表面单元)，将点集II_Cu中铜材料外的高斯点剔除，剩下的点形成的集合则为铜材料域的质点集合III_cu，如图1实心小三角形所示。同理，对铝材料区域进行类似的步骤，获得铝材料的质点集合III_Al，如图1实心小正方形所示。最终，将两集合III_Cu和III_Al合并为坯饼的物质点集III，并为每个物质点赋予相对应的材料属性，如密度，体积，泊松比，屈服应力、硬化指数、强化系数等参数。为了准确描述硬币初始轮廓,将坯饼表面各三角形单元的高斯点集添加到坯饼质点集合中，组成新的坯饼质点集。如图3所示为其四分之一模型被切除，以方便查看内部质点。将坯饼表面质点集分为三部分：上表面、侧面和下表面质点集，它们分别与上模，中圈和下模构成接触对。在仿真过程中，此三对接触对分别进行接触判断处理。

步骤5：计算背景网格节点载荷；

节点内力矢量为：

式中，

是节点I内力第i个分量，m_p为质点p的质量，σ_ij为质点p的应力，

为在质点p处节点I的形函数的第j个偏导数，

为质点p在n时刻的位置，ρ_p为质点p的密度。

节点外力矢量为：

式中，

是节点外力，b_i是体力，h为边界层的厚度。针对压印成形问题，在不考虑接触的情形下，此节点外力为0。

n时刻节点I处的合力为：

式子中，

是n时刻节点I处的合力。

步骤6：通过插值形函数将质点携带的物理量映射到背景网格节点；

节点I质量为：

节点动量以及节点速度为：

式子中，

是节点I的形函数第i个分量，

是n-1/2时刻第p个质点的速度第i个分量，

是n+1/2时刻节点I的动量第i个分量，

是n+1/2时刻节点I的速度第i个分量，

表示节点I在n时刻的质量。

步骤7：根据双色币不同材料特性求解系统自适应时间步长；

在通常情况下，针对单一材料压印成形仿真问题，我们采用中心差分动力显式算法进行求解。该算法为条件稳定，系统时间步长Δt必须小于临界步长Δt_cr，即

Δt＝αΔt_cr (13)

式子中，α是一个常数，取值范围为：0.8≤α≤0.98

由于CFL(Courant Friedrichs Lewy，柯朗-弗里德里希斯-列维)条件，为了保证显式时间积分的稳定性，采用有限元法进行动态显式分析时，一个时间步长内波传输的距离必须不超过一个单元。而弹性材料的声速为：

式子中，E为弹性模量，v为泊松比，ρ为材料密度。假定单元的名义长度为l_e，则临界时间步长为：

随变形的增加，单元名义尺寸逐渐减小，临界时间步长也逐渐减小。当网格发生严重畸变时，所求解的最小名义长度趋近于0，导致时间步长亦趋近于0，从而仿真分析无法进行。

物质点法与有限元法不同，所有时间积分都是需要在固定的背景网格上进行。在确定临界时间步长时，不仅需要考虑材料的声速，同时需要考虑质点的速度。特别地，在超高速撞击问题中，质点速度较大，与材料声速在同一数量级甚至超过材料声速。因此，质点速度对仿真计算的影响不可忽略。在压印成形模拟过程中，基于动力显式中心差分算法的物质点法求解，在每一步需要计算临界时间步长使得数值求解稳定：

式子中，c为质点的声速，u为质点速度。在使用均匀背景网格时，l_g又可以取为d_c，即背景网格尺寸。本发明中坯饼有两种不同材料，即内铜外铝。时间步长的大小由系统的属性和模型材料参数所决定。由于涉及两种不同材料，需要分别计算两种物质区域内的时间步长并取其较小值。铜质材料中时间步长为：

式子中，

和

分别为质点集III_Cu中声速和任意质点p₁质点速度。

同理，铝质材料中时间步长为：

式子中，

和

分别为质点III_Al中声速和任意质点p₂的质点速度。

我们需要计算出两种材料的波速，以及压印过程中两种材料物质点的速度，进而比较两者的临界时间步长，最终确定系统的临界时间步长为两者的最小值。因此，系统的临界时间步长为：

式子中，

为Cu的时间步长，

为Al的时间步长。

步骤8：求解动量方程并施加边界条件；

积分背景网格节点动量方程：

式中，

和

分别为节点I在n+1/2和n-1/2时刻的第i个方向的动量。

步骤9：更新质点的物理量；

更新物质点速度：

表示在p质点在n+1/2时刻的速度第i个分量。

更新物质点位置：

式子中，

表示p质点在n+1时刻的坐标第i个分量，

表示p质点在n时刻的坐标第i个分量。

步骤10：更新背景网格点的物理量；

更新背景网格节点速度

式子中，

表示节点I在n+1/2时刻的速度。

步骤11：更新模具空间位置；

更新模具的位置：

d_punch＝d_punch+Δd_punch (23)

d_punch表示上模的当前总位移，Δd_punch表示上模的当前时刻的位移。

步骤12：判断物质点的穿透状态及求解穿透距离；

如图4所示，根据传统物质点法理论，模具质点j′与坯饼质点j对背景网格结点I的速度均有贡献(因为它们各自所在的背景网格共享结点I)，因而判定模具与该坯饼质点接触。而实际上它们并没有接触，从而产生虚假接触。本发明采用点面接触算法替代传统物质点的点-点接触算法以避免虚拟接触，提高仿真计算准确性。判断坯饼表面每个质点与上模、下模和中圈之间的接触状态和计算接触值(接触力、穿透深度和摩擦力)时，在每一个时间步长里都会造成很大的计算量。为了减少接触判断的计算时间，采用了全局搜索与局部搜索相结合的策略。全局搜索是为每个质点找到与之相接触的候选模具单元，过程如下所示。

构建包含特定模具网格的块：块的大小根据各模具的最大和最小空间坐标来定义。块的左下角和右下角的坐标分别为(x_min，y_min，z_min)和(x_max，y_max，z_max)。然后将该块区域划分为子单元格，每个方向上的单元格数目分别为N_x，N_y，N_z，并分别为每个单元格沿x,y,z三个方向进行编号。注意这个块随模具网格一起运动。为了区分这种单元格和背景单元网格，我们将这种单元格命名为空间单元格。

搜索每个单元格所包含的候选的模具网格单元。如图5所示为某二维物体的单元格示意图。共有N＝N_x×N_y＝8×5＝40个空间单元格，图中只标注了四个角的单元格编号。编号为21的单元格包含五个模具网格单元，其节点用实心圆表示。查找每个单元格所包含的候选单元一般包括两个步骤。首先，将某模具空间格域沿三个方向扩展，扩展的距离为子空间格尺寸的一半，并循环每个模具单元以获取包含它的单元格，最终得到每个单元格包含的候选模具单元。对于给定的粒子x_p＝(x，y，z)，其定位的单元编号为

其中[]是一个四舍五入的运算操作。针对质点x_p，可以定位到包含该质点的子空间格。根据上述的全局搜索结果可以得到该质点的候选模具单元集合。接下来在该候选集中进行局部搜索并进行接触判断。为了判断网格之间的接触状态和质点x_p，我们只需要找到质点是否接触候选单元或者穿透它。执行以下步骤可以实现此目标。

求潜在接触单元：如图6所示，质点x_p，S₁，S₂，S₃，S₄有四个候选单元。找到离x_p最近的节点m_s，它与x_p不重合。如果满足下列不等式，则质点将与S_i单元接触

(C_i×g₂)·(C_i×C_i+1)>0，(C_i×g₂)·(g₂×C_i+1)>0 (25)

这里的i＝1,2,3,4为待判断的第i个模具单元。C_i和C_i+1是两个单位向量，它们有共同的起始点，沿着S_i单元的两条相邻边，g₁从m_s开始，到x_p结束。g₂是在单元S_i在向量g₁上的投影，写成：

g₂＝g₁-(g₁·m)m (26)

这里

如图7所示，为了计算x_c的空间坐标，在单元上建立参数公式：

r(ξ，η)＝f₁(ξ，η)i₁+f₂(ξ，η)i₂+f₃(ξ，η)i₃ (28)

其中

N_j(ξ，η)是壳体单元的双线性函数，

为节点j的第i个坐标，(ξ，η)为单元内任意一点的局部坐标，i₁，i₂，i₃分别为空间坐标下的单位向量。t是从原始坐标开始，到x_p结束的向量。因为向量

是垂直于接触单元，接触点x的局部坐标(ξ_c，η_c)必须满足：

判断接触状态和穿透深度d的计算：

其中n_c为在点(ξ_c，η_c)接触单元的单位向外的法向量，计算方法为：

步骤13：求解接触力；

如图8所示，阻力F_R是关于d的一个函数，其表达式为：

F_R＝-km_pdn_c/(△tⁿ)² (33)

其中m_p为质点质量，k为界面刚度。

施加在质点上法向和切向接触力的表达式为：

式子中，

为质点p的法向接触力，

为质点p的切向接触力，

为阻力，

为均为接触点的法向量。总接触力可表示为：

其中τ_ip为接触点的单位切向矢量，μ为摩擦系数，

为质点p的总接触力。

步骤14：再次更新质点的物理量：

质点x_p的更新坐标可以写成：

式子中，

是p质点在n+1/2时刻的速度，

为质点p的总接触力，m_p为质点p的质量，

是p质点在n+1/2时刻的坐标。

最后更新背景网格动量：

式子中，

为节点I在n+1/2时刻的动量，n_p为质点个数，N_Ip质点p的形函数，

是p质点在n+1/2时刻的速度。

步骤15：施加边界条件：

对于固定边上的节点，令

步骤16：计算质点的应变增量和旋率增量，根据本构模型求解质点的应力及弹塑性相关变量，更新所有质点的密度。计算物质点的应变增量和旋率增量：

式子中，

是应变增量，

是旋率增量，

是n+1/2时刻第I个节点的速度,

是n+1/2时刻第I个节点的速度。

更新物质点的密度：

式子中，

是n+1时刻第p个质点的密度，

是n时刻第p个质点的密度。综上，本发明公开了一种应用于双色币压印成形仿真的基于改进的接触算法的物质点法，该方法包括以下步骤：

步骤1：对双色币压印成形工艺过程建立动力学物理模型；

步骤2：对压印成形过程动力学模型进行物质点积分；

步骤3：对成形模具进行有限元三角形/四边形单元离散；

步骤4：根据坯饼数据模型建立背景网格，并对双色币坯饼进行四面体/六面体单元离散、建立铜、铝材料区域的质点集；

步骤5：计算背景网格节点载荷；

步骤6：通过插值形函数将质点携带的物理量映射到背景网格节点；

步骤7：根据双色币不同材料特性求解系统自适应时间步长；

步骤8：求解动量方程并施加边界条件；

步骤9：更新质点的物理量；

步骤10：更新背景网格点的物理量；

步骤11：更新模具空间位置；

步骤12：判断物质点的穿透状态及求解穿透距离；

步骤13：求解接触力；

步骤14：再次更新质点的物理量；

步骤15：施加边界条件；

步骤16：计算质点的应变增量和旋率增量，根据本构模型求解质点的应力及弹塑性相关变量，更新所有物质点的密度。

传统有限元算法在处理实际压印成形或者类似有模具的加工问题时，材料界面接触算法需要额外建立材料自接触判断处理。本发明提出改进点-面接触算法的物质点法，不仅避免了虚拟接触以提高物质点法接触判断的准确性，同时，不需要进行额外的材料自接触判断算法，为此类涉及材料自接触问题提供了有效途径。

Claims (6)

1.一种应用于双色币压印成形仿真的基于改进接触算法的物质点法，其特征在于：包括如下步骤：

步骤1：对双色币压印成形工艺过程建立动力学物理模型；

步骤2：对压印成形过程动力学模型进行物质点积分；

步骤3：对成形模具进行有限元三角形/四边形单元离散；

步骤4：根据坯饼数据模型建立背景网格，并对双色币坯饼进行四面体/六面体单元离散、建立铜、铝材料区域的质点集；

步骤5：计算背景网格节点载荷；

步骤6：通过插值形函数将质点携带的物理量映射到背景网格节点；

步骤7：根据双色币不同材料特性求解系统自适应时间步长；

步骤8：求解动量方程并施加边界条件；

步骤9：更新质点的物理量；

步骤10：更新背景网格点的物理量；

步骤11：更新模具空间位置；

步骤12：判断物质点的穿透状态及求解穿透距离；

步骤13：求解接触力；

步骤14：再次更新质点的物理量；

步骤15：施加边界条件；

步骤16：计算质点的应变增量和旋率增量，根据本构模型求解质点的应力及弹塑性相关变量，更新所有物质点的密度。

2.根据权利要求1所述的一种应用于双色币压印成形模拟中的基于改进的接触算法的物质点法，其特征在于：步骤1的主要内容为建立双色币压印成形工艺过程建立动力学物理模型，具体内容如下：

从力学角度出发，压印成形过程是一个拟静态的低速运动过程，采用动力显式算法来描述及求解压印成形过程，设坯饼为被研究物体，其发生塑性变形时应满足基本平衡方程：

其中，σ_ij，j分别为柯西应力张量，b_i为体积力，ρ为材料密度，γ为阻尼系数，

和

分别是物体内任一点的速度和加速度，i和j分别为坐标方向，根据虚位移原理、相关应力边界条件和高斯散度定理，可推导系统的虚功率方程为：

式中，V为物体区域，

是虚速度，

是对应于柯西应力σ_ij的虚应变速率。其中S_p为已知外力p_i的表面；S_c为与另一物体接触的表面，记接触表面力为q_i；惯性力和阻尼力功率项反映了物体系统的惯性效应和物理阻尼效应。

3.根据权利要求1所述的一种应用于双色币压印成形模拟中的基于改进的接触算法的物质点法，其特征在于：步骤2为压印成形过程动力学方程的物质点积分，具体步骤如下：

物质点法中，每一个质点占据一定的体积空间和携带一定质量以及其他相关物理量，任意空间点x的材料密度可以表示为：

其中n_p为质点数，m_p为质点质量，δ()为Dirac函数，x_p为质点p的坐标，背景网格中任意质点p的位移u_ip可通过背景网格结点形函数插值获取：

其中，N_Ip为质点p相关的结点I的形函数，u_iI为背景网格结点I的第i方向的速度，将上述密度和位移公式带入步骤1中的系统虚功率方程，同时考虑压印成形过程中无体积力b_i和外力p_i，可得：

其中，

为结点I质量，

和

分别为结点I第i个方向的速度和加速度，

为结点I第i个方向的外力，

为结点I第i个方向的内力。

4.根据权利要求1所述的一种应用于双色币压印成形模拟中的基于改进的接触算法的物质点法，其特征在于：步骤4的主要内容为提取铜、铝两种不同材料区域(铜为内芯，铝为外环)的物质点集，具体内容如下：

首先，对坯饼所占区域划分规则六面体背景网格，沿x，y，z方向的网格数分别为Gx，Gy，Gz，并求出每个背景网格的积分点，然后，对铜、铝两种不同材料区域分别进行四面体网格划分，每一个小矩形表示背景网格单元，实心圆代表网格节点，小实心三角形代表铜物质点，小实心正方形代表铝物质点；接着，在铜材料区域找出每一个四面体单元的背景网格单元集合I_Cu，求该集合的具体方法如下：将背景网格单元沿x，y，z方向从1开始编号至G＝Gx*Gy*Gz，根据四面体单元的四个节点坐标，可以找出一个空间长方体来刚好包围该四面体，假设背景网格最左前点坐标(x_min，y_min，z_min)和最右后点坐标(x_max，y_max，z_max)，依据空间任意一点x_p＝(x，y，z)的背景网格编号

找到该长方体所包含的所有背景网格集合I_Cu，将该集合中所有单元的高斯点组成的集合定义为点集II_Cu，根据铜材料的表面边界单元，由网格剖分系统自动生成并编写程序读取该表面单元，将点集II_Cu中铜材料外的高斯点剔除，剩下的点形成的集合则为铜材料域的质点集合III_cu，同理，对铝材料区域进行类似的步骤，获得铝材料的质点集合III_Al，最终，将两集合III_Cu和III_Al合并为坯饼的物质点集III，并为每个物质点赋予相对应的材料属性，如密度，体积，泊松比，屈服应力、硬化指数、强化系数等参数，为了准确描述硬币初始轮廓，将坯饼表面各三角形单元的高斯点集添加到坯饼质点集合中，组成新的坯饼质点集，将坯饼表面质点集分为三部分：上表面、侧面和下表面质点集，它们分别与上模，中圈和下模构成接触对；在仿真过程中，此三对接触对分别进行接触判断处理。

5.根据权利要求1所述的一种应用于双色币压印成形仿真中的基于改进的接触算法的物质点法，其特征在于：步骤7为动力显式中心差分算法的时间积分步长确定，具体包含以下内容：

在通常情况下，针对单一材料压印成形仿真问题，采用中心差分动力显式算法进行求解，该算法为条件稳定，系统时间步长Δt必须小于临界步长Δt_cr，即：

Δt＝αΔt_cr

式子中，α是一个常数，取值范围为：0.8≤α≤0.98；

由于CFL(Courant Friedrichs Lewy，柯朗-弗里德里希斯-列维)条件，为了保证显式时间积分的稳定性，采用有限元法进行动态显式分析时，一个时间步长内波传输的距离必须不超过一个单元，而弹性材料的声速为：

式子中，E为弹性模量，v为泊松比，ρ为材料密度，假定单元的名义长度为l_e，则临界时间步长为：

随变形的增加，单元名义尺寸逐渐减小，临界时间步长也逐渐减小，当网格发生严重畸变时，所求解的最小名义长度趋近于0，导致时间步长亦趋近于0，从而仿真分析无法进行；

物质点法与有限元法不同，所有时间积分都是需要在固定的背景网格上进行，在确定临界时间步长时，不仅需要考虑材料的声速，同时需要考虑质点的速度；特别地，在超高速撞击问题中，质点速度较大，与材料声速在同一数量级甚至超过材料声速；因此，质点速度对仿真计算的影响不可忽略，在压印成形模拟过程中，基于动力显式中心差分算法的物质点法求解，在每一步需要计算临界时间步长使得数值求解稳定：

式子中，c为质点的声速，u为质点速度，在使用均匀背景网格时，l_g又可以取为d_c，即背景网格尺寸；本发明中坯饼有两种不同材料，即内铜外铝，时间步长的大小由系统的属性和模型材料参数所决定，由于涉及两种不同材料，需要分别计算两种物质区域内的时间步长并取其较小值，铜质材料中时间步长为：

式子中，

和

分别为质点集III_Cu中声速和任意质点p₁质点速度；

同理，铝质材料中时间步长为：

式子中，

和

分别为质点III_Al中声速和任意质点p₂的质点速度；

需要计算出两种材料的波速，以及压印过程中两种材料物质点的速度，进而比较两者的临界时间步长，最终确定系统的临界时间步长为两者最小值，因此，系统的临界时间步长为：

式子中，

为Cu的时间步长，

为Al的时间步长。

6.根据权利要求1所述的一种应用于双色币压印成形仿真中的基于改进的接触算法的物质点法，其特征在于：步骤12为物质点穿透状态判断及求解接触力，具体包含以下内容：

根据传统物质点法理论，模具质点j′与坯饼质点j对背景网格结点I的速度均有贡献(因为它们各自所在的背景网格共享结点I)，因而判定模具与该坯饼质点接触，而实际上它们并没有接触，从而产生虚假接触，本发明采用点面接触算法替代传统物质点的点-点接触算法以避免虚拟接触，提高仿真计算准确性，判断坯饼表面每个质点与上模、下模和中圈之间的接触状态和计算接触值(接触力、穿透深度和摩擦力)时，在每一个时间步长里都会造成很大的计算量，为了减少接触判断的计算时间，采用了全局搜索与局部搜索相结合的策略，全局搜索是为每个质点找到与之相接触的候选模具单元，过程如下所示：

构建包含特定模具网格的块：块的大小根据各模具的最大和最小空间坐标来定义，块的左下角和右下角的坐标分别为(x_min，y_min，z_min)和(x_max，y_max，z_max)，然后将该块区域划分为子单元格，每个方向上的单元格数目分别为N_x，N_y，N_z，并分别为每个单元格沿x，y，z三个方向进行编号；注意这个块随模具网格一起运动，为了区分这种单元格和背景单元网格，将这种单元格命名为空间单元格；

搜索每个单元格所包含的候选的模具网格单元，共有N＝N_x×N_y＝8×5＝40个空间单元格，图中只标注了四个角的单元格编号；编号为21的单元格包含五个模具网格单元，其节点用实心圆表示，查找每个单元格所包含的候选单元一般包括两个步骤；首先，将某模具空间格域沿三个方向扩展，扩展的距离为子空间格尺寸的一半，并循环每个模具单元以获取包含它的单元格，最终得到每个单元格包含的候选模具单元；对于给定的粒子x_p＝(x，y，z)，其定位的单元编号为

其中[]是一个四舍五入的运算操作，针对质点x_p，可以定位到包含该质点的子空间格，根据上述的全局搜索结果可以得到该质点的候选模具单元集合；接下来在该候选集中进行局部搜索并进行接触判断，为了判断网格之间的接触状态和质点x_p，只需要找到质点是否接触候选单元或者穿透它，执行以下步骤可以实现此目标；

求潜在接触单元：质点x_p，S₁，S₂，S₃，S₄有四个候选单元，找到离x_p最近的节点m_s，它与x_p不重合，如果满足下列不等式，则质点将与S_i单元接触

(C_i×g₂)·(C_i×C_i+1)＞0，(C_i×g₂)·(g₂×C_i+1)＞0

这里的i＝1，2，3，4为待判断的第i个模具单元，C_i和C_i+1是两个单位向量，它们有共同的起始点，沿着S_i单元的两条相邻边，g₁从m_s开始，到x_p结束，g₂是在单元S_i在向量g₁上的投影，写成：

g₂＝g₁-(g₁·m)m

这里

为了计算x_c的空间坐标，在单元上建立参数公式：

r(ξ，η)＝f₁(ξ，η)i₁+f₂(ξ，η)i₂+f₃(ξ，η)i₃

其中

N_j(ξ，η)是壳体单元的双线性函数，

为节点j的第i个坐标，(ξ，η)为单元内任意一点的局部坐标，i₁，i₂，i₃分别为空间坐标下的单位向量；t是从原始坐标开始，到x_p结束的向量，因为向量

是垂直于接触单元，接触点x的局部坐标(ξ_c，η_c)必须满足：

判断接触状态和穿透深度d的计算：

d＝n_c·[t-r(ξ_c，η_c)]

其中n_c为在点(ξ_c，η_c)接触单元的单位向外的法向量，计算方法为：

Images