一种模拟压制过程中粉末转化为连续体的方法

本发明涉及粉末成型与数值模拟，具体涉及一种模拟压制过程中粉末转化为连续体的方法。

背景技术：

1、粉末的压制过程是在高压下把造型粉压制成连续体的过程，涉及的造型粉尺寸为几十微米至毫米量级，连续体尺寸为厘米至米量级，压制压力为兆帕至百兆帕量级，压制时长为分钟至小时量级。

2、造型粉在加载过程中，离散颗粒体系逐渐出现运动、变形、压密及粘接现象，从而实现从离散体到连续体的渐进转变。有限元法及离散元法是两类主要的压制成型数值模拟方法。有限元法重点关注压制过程中宏观力学行为的变化，如载荷-位移曲线的演化、空间密度的演化等，该方法的计算精度主要依赖于宏观本构的精度，但宏观本构的建立及本构参数的获取则主要依赖于试验，由于宏观试验数据往往有限，因此无法对压制过程中造型粉的宏观演化进行精确表征，更无法直接根据造型粉的细观参数对压制过程中的宏观力学行为进行直接预测。离散元法重点关注造型粉间的细观行为，通过精细建立造型粉的细观模型并引入细观参数，可对压制过程中造型粉的运动、变形及挤密行为进行精准模拟；但受限于离散元的计算量，该类算法仅适用于小尺度的压制机理分析。

3、总体而言，当前已有数值计算方法或只关注连续部分的计算（如有限元法），或只进行离散颗粒的描述（如离散元法），缺乏定量分析压制过程中造型粉从离散态到连续态转变的计算方法。

技术实现思路

1、针对上述背景技术中存在的技术问题，本发明提出了一种模拟压制过程中粉末转化为连续体的方法，其可真实复现粉末压制过程中离散-连续渐进转变过程，具有物理一致性好、仿真精度高、等优点，可以解决现有技术中，细观仿真计算效率低、宏观仿真无法考虑细观结构影响的技术问题。

2、为解决上述技术问题，本发明提供的一种模拟压制过程中粉末转化为连续体的方法，其包括如下步骤：

3、（1）将粉末压制区域离散为离散元及有限元两个区域，两个区域完全重叠；

4、（2）选取空间映射模式，将离散元区域中的离散元颗粒映射到对应的有限元区域中的有限元单元上；

5、（3）在离散元的颗粒间引入接触本构，在有限元区域的有限元的单元上引入粘弹塑性本构；

6、（4）选定离散元转化为有限元的转化模型，并设定转化参数；

7、（5）在粉末四周创建模具及上下压板；

8、（6）根据设定的加载模式，选择在设定的压板上施加准静态速度载荷；

9、（7）计算离散元与有限元的耦合转化，根据叠合区域内离散元颗粒的转化因子计算有限元单元的转化因子；

10、（8）将粉末压制区域作为一个模型，根据离散元颗粒的转化因子及有限元单元的转化因子，计算压制区域的整体受力。

11、作为本发明的一种优选方案，所述有限元单元的最小网格尺寸≥离散元颗粒的最大颗粒直径，具体表示为：

12、；

13、上式中，为第个单元的特征尺寸；为第个颗粒的直径。

14、作为本发明的一种优选方案：所述步骤（2）中的空间映射模式包括体心映射模式和面积映射模式两类；

15、所述空间映射模式为体心映射模式时，若颗粒体心位于某个有限元单元的内部，该颗粒即和该单元形成绑定关系；设有限元单元由个面组成，第个平面的方程为：

16、；

17、上式中，为平面的单位外法向量，为常数；

18、若满足下式，则表示颗粒体心位于某个有限元单元的内部：

19、；

20、所述空间映射模式为面积映射模式时，则需计算离散元颗粒与每个有限元单元的交集，若两者的交集不为空，则该离散元颗粒即和该有限元单元形成绑定关系；设离散元颗粒占据的空间区域为，有限元单元占据的空间区域为；若，则定义离散元颗粒与每个有限元单元的交集体积为，即：

21、；

22、上式中，为离散元颗粒的指示函数，为有限元单元的指示函数；当交集为空时，该积分值为0。

23、作为本发明的一种优选方案：所述步骤（3）中的接触本构为线性接触模型、hertz接触模型和dpc接触模型中的任意一种；所述步骤（3）中的粘弹塑性本构为dpc模型、mc模型和burgers模型中的任意一种。

24、作为本发明的一种优选方案，所述步骤（4）中离散元转化为有限元的转化模型为体应力转化模型及体积模量转化模型中的任意一种；

25、在体应力转化模型中，离散元颗粒转化因子的计算公式表示为：

26、；

27、上式中，为开始转变时的体应力，为终止转变时的体应力；

28、当某个离散元颗粒所受到的体应力≥，开始转化；当某个离散元颗粒所受到的体应力≥，转化完毕；

29、在体积模量转化模型中，离散元颗粒转化因子的计算公式表示为：

30、；

31、上式中，为开始转变时的离散元颗粒所在区域等效体积模量，为终止转变时的离散元颗粒所在区域等效体积模量；

32、当某个离散元颗粒所在区域的等效体积模量≥，开始转化；当某个离散元颗粒所在区域的等效体积模量≥，转化完毕。

33、作为本发明的一种优选方案，所述步骤（7）中离散元的计算包括离散元颗粒接触力的计算、离散元颗粒的受力计算和离散元颗粒的运动计算；

34、所述离散元颗粒的接触力计算的具体过程为：

35、颗粒间接触力计算包含法向力计算及切向摩擦计算两部分，采用线性接触模型时的法向力计算公式为：

36、；

37、上式中，、为法向接触刚度及法向阻尼系数，、为嵌入量及嵌入速率；

38、切向力计算时采用增量-更新模式，每个时间步内先计算试探切向力，再根据库仑摩擦定律修正；试探应力计算公式为：

39、；

40、上式中，为当前时步试探切向力，为上一时步切向力，为切向刚度，为切向位移差的增量；

41、接着计算当前时步最大静摩擦力，为：

42、；

43、上式中，为最大静摩擦力，为摩擦系数；

44、接着进行试探切向力的修正，计算公式为：

45、；

46、所述离散元颗粒的受力计算的具体过程为：

47、计算离散元颗粒所受合力及合力矩，计算公式为：

48、；

49、上式中，及分别依次为颗粒所受的合力及合力矩；及分别依次为施加于颗粒的重力及力矩；及分别依次为颗粒之间的接触力及连接力矩；为颗粒与刚性面之间的接触力；为颗粒与颗粒之间的接触作用点相对于颗粒质心的位置矢量；为颗粒与刚性面间的作用点相对于颗粒质心的位置矢量；

50、所述离散元颗粒的运动计算，包含颗粒平动计算及颗粒转动计算：

51、颗粒平动计算的计算公式为：

52、；

53、上式中，、、分别依次为加速度、速度及位移矢量，为颗粒质量，为平动局部阻尼系数；

54、颗粒转动计算的计算公式为：

55、；

56、上式中，、、分别依次为转动角加速度、角速度及角位移；为转动局部阻尼系数；为转动惯量；用转动速度的四元数表征量。

57、作为本发明的一种优选方案，所述步骤（7）中有限元的计算包含有限元单元变形受力的计算及有限元节点受力运动的计算；

58、所述有限元单元变形受力的计算过程如下：

59、根据节点位移增量计算某个单元某个高斯点的应变增量，计算公式为：

60、；

61、上式中，为某个单元某个高斯点的应变增量矩阵，为高斯点的应变矩阵，为单元应变增量矩阵。

62、根据应变增量计算当前时刻的应力增量，计算公式为：

63、；

64、上式中，为材料矩阵，由该固体单元所采用的本构模型确定。

65、根据当前时刻的应力增量计算应力全量，计算公式为：

66、；

67、上式中，为上一时步的应力全量；

68、所述有限元节点受力运动的计算过程如下：

69、首先，计算有限元单元每个节点的合力，计算公式为：

70、；

71、上式中，为有限元节点所受到的合外力，为有限元单元变形对该节点产生的力，为该节点上施加的外力；

72、然后根据每个节点的合力计算节点的加速度；若节点当前速度大于零，则加速度；若节点当前速度小于等于零，则加速度；其中，为局部阻尼系数，为节点质量；

73、接着根据加速度更新速度，再根据速度更新本步的位移增量和位移，计算公式为：

74、；

75、上式中，为节点加速度矢量；为节点速度矢量；为节点位移增量矢量；为节点位移矢量；为计算时步。

76、作为本发明的一种优选方案：所述准静态速度载荷可在上压板上施加，也可在上下压板上同时施加，准静态载荷的范围设定为10-7~10-9m/step。

77、作为本发明的一种优选方案，所述有限元单元的转化因子的计算公式表示为：

78、；

79、上式中，为与某个有限元单元有绑定关系的离散元颗粒总数，为第个离散元颗粒的体积，为第个离散元颗粒的转化因子。

80、作为本发明的一种优选方案，所述计算模型整体受力的公式为：

81、；

82、上式中，为某个有限元单元节点所在区域受到的合力，为直接由有限元计算获得的节点力，为直接由离散元颗粒计算得到的节点力，为第i个离散元颗粒在第j个有限元节点上的权系数；

83、计算时，首先根据有限元单元中各高斯点的应力计算有限元单元各节点的节点力矩阵，计算公式为：

84、；

85、上式中，为高斯点的应变矩阵，为当前时步的应力全量，及为高斯点的权系数及雅可比行列式的值；

86、接着，在节点力矩阵中找到节点对应的节点力向量，计算公式为：

87、；

88、颗粒所受到的其他颗粒给其的接触力的计算公式为：

89、；

90、上式中，为与颗粒接触的其他颗粒的数量，为第个与颗粒接触的颗粒作用到颗粒上的力。

91、作为本发明的一种优选方案：所述粉末转化为连续体的方法可进行二维问题仿真或三维问题仿真；所述数值方法进行二维仿真时，所述离散元颗粒为圆、椭圆或多边形，所述有限元单元为三角形或四边形；所述数值方法进行三维仿真时，所述离散元颗粒为球、椭球或多面体，所述有限元单元为四面体、金字塔、三棱柱或六面体。

92、采用上述技术方案，本发明具有如下有益效果：

93、本发明将离散元方法与有限元方法进行耦合，有效克服了传统数值模拟中无法实现细宏观一体化计算的局限性；借助离散元向有限元的渐进转化策略，确保了转化过程自然、无震荡、鲁棒性强，可真实复现从松散颗粒到致密连续体的全过程力学响应。本发明通过构建离散元与有限元在空间上的叠合覆盖模型，以体应力及体积模量作为转化判据，采用转化因子动态分配离散元与有限元对空间区域的贡献程度，实现松散粉末介质到致密连续介质的渐进平滑转化。

94、本发明的主要特点和优点主要体现在以下几方面：：

95、（1）实现了细观-宏观一体化高保真模拟：本发明创新性地将离散元(dem)与有限元(fem)在空间上完全叠合，并通过转化因子实现二者的动态耦合。该方法既能利用离散元法精准捕捉压制初期粉末颗粒的运动、重排及挤压破碎等细观力学行为，又能通过渐进转化为有限元法高效模拟压制后期致密连续体的宏观应力-应变响应，突破了传统单一数值方法无法同时兼顾细观机理与宏观性能的局限。

96、（2）转化过程物理意义明确且数值稳定性高：本发明提出的基于体应力或体积模量的转化模型，具有明确的物理内涵。转化因子在0到1之间平滑过渡，避免了离散元与有限元之间切换时的非物理震荡或数值突变。同时，通过转化因子对节点合力进行加权平均，确保了转化过程中系统力学响应的连续性，算法鲁棒性强，能够真实复现从松散粉末到致密连续体的全过程力学演化。

97、（3）适用性强，可模拟复杂几何与多元材料体系：本发明提供了两种灵活的空间映射模式，并支持多种离散元接触本构与有限元本构的组合。该方法支持二维和三维问题分析，允许粉末颗粒为球形、椭球、多边形或多面体，能够灵活适应不同形状粉末、复杂模具几何以及不同材料体系的压制过程模拟。