一种SLM成形过程温度场数值模拟方法与流程

本发明属于增材制造计算机数值模拟技术领域，具体涉及一种slm成形过程温度场数值模拟方法。

背景技术：

slm成形过程受工艺参数影响，成形件容易发生卷曲、球化、热变形和热裂纹等缺陷。slm成形产品结构与材料需求多变，现有研究主要通过大量反复实验确定工艺参数，普适性较差。目前技术人员要求高、准备时间长以及制造失败率高的现状，使得slm技术难以大范围普及应用。

通过计算机数值模拟获取slm成形过程的温度场，特别是熔池温度，成为指导slm成形的重要手段。在成形过程中能否保持稳定且尺寸适中的熔池，是成形件质量是否过关的重要因素。此外，温度场数值模拟所得温度数据还为模拟应力提供了数据基础。

slm成形过程热源激光速度较快，温度梯度大。为准确模拟成形过程的温度变化，现有温度场数值模拟方法在模拟过程中常采用小尺寸的离散单元(边长0.1mm左右)，则尺寸仅为10mm*10mm*10mm的样件离散就将会得到10⁶个离散单元。现有温度场数值模拟方法要求较高的计算机配置，模拟小尺寸样件(如10mm*10mm*10mm)仍需要大量计算时间(一周左右)。样件尺寸越大，则计算时间成本越大。这一缺陷限制了slm成形过程数值模拟的发展与应用。有方法通过整体增大离散单元的尺寸缩减计算量，但又带来了计算精度不足的问题，因此在保证计算精度的情况下降低计算量是十分必要的。

技术实现要素：

本发明提供了一种slm成形过程温度场数值模拟方法，能够减小slm温度场模拟计算的计算量，为厘米级大尺寸样件模拟提供了新的解决办法。

为达到上述目的，本发明所述一种slm成形过程温度场数值模拟方法，包括以下步骤：

步骤1，建立仿真模型：在三维建模软件中建立基板和样件的三维模型，对样件和基板分别赋予材料属性；

步骤2，离散有限差分单元：将步骤1得到的基板和样件的三维模型按照设定的尺寸离散得到若干正六面体有限差分单元el单元，将设定区域内的el单元细化分割得到比el单元尺寸小的em单元，将部分em单元细化分割得到比em单元尺寸小的es单元；el单元、em单元和es单元组成嵌套单元；

步骤3，建立传热计算模型，确定最大允许计算时间步长δtj：建立同类型单元之间的换热方程，根据步骤2确定的el单元、em单元及em单元的尺寸，结合基板材料属性、样件材料属性以及粉末的热物性参数计算同类型单元之间传热的最大允许计算时间步长δtj，以保证传热计算收敛；

步骤4，确定传热策略：根据激光的扫描速度确定激光在步骤2所得es单元之间的移动时间，并根据其与步骤3得到的各单元之间的最大允许计算时间步长的大小关系，确定传热计算策略，并对步骤3确定的最大允许计算时间步长δtj进行修正；

步骤5，计算温度场：按照步骤4确定的传热计算策略开始模拟slm成形过程，即可得到slm成形过程温度场的数值模拟结果，计算结果包含整个成形过程中各类型单元的温度变化历程。

进一步的，在步骤2中，将基板的三维模型和样件的三维模型按照长为δx，宽为δy，高为δz的正六面体el单元离散，el单元各方向尺寸为热源激光光斑直径的2-10倍；定义em单元是长为δx，宽为δy，高为δz的正六面体单元；定义es单元是长为δx，宽为δy，高为δz的正六面体单元，δx与δy大小与热源激光光斑直径相同；定义n是正六面体el单元与正六面体es单元的尺寸比，即n为正整数。

进一步的，在模拟slm成形时，所述设定区域为：所有z方向坐标等于激光所在el单元坐标zc以及等于zc-1的el单元；将连同激光所处el单元在内的4×2×2个el单元1所包含的em单元再次细化分割为若干es单元。

进一步的，步骤2中，被划分为es单元的区域跟随激光移动而改变。

进一步的，步骤2中，δz尺寸与实际打印铺粉层的厚相同。

进一步的，步骤3中，slm成形过程中的空间内任一单元位置用其几何中心的点坐标表示，其与周围单元的有限差分换热公式表示为：

上式中分别是i、i+1时刻坐标为(x,y,z)的单元的温度，是i时刻坐标为(x+1,y,z)的单元的温度，是i时刻坐标为(x-1,y,z)的单元的温度，是i时刻坐标为(x,y+1,z)的单元的温度，是i时刻坐标为(x,y-1,z)的单元的温度，是i时刻坐标为(x,y,z+1)的单元的温度，是i时刻坐标为(x,y,z-1)的单元的温度；λ1是坐标为(x,y,z)的单元与坐标为(x+1,y,z)的单元材料的平均导热系数，λ2是坐标为(x,y,z)的单元与坐标为(x-1,y,z)的单元材料的平均导热系数，λ3是坐标为(x,y,z)的单元与坐标为(x,y+1,z)的单元材料的平均导热系数，λ4是坐标为(x,y,z)的单元与坐标为(x,y-1,z)的单元材料的平均导热系数，λ5是坐标为(x,y,z)的单元与坐标为(x,y,z+1)的单元材料的平均导热系数，λ6是坐标为(x,y,z)的单元与坐标为(x,y,z+1)的单元材料的平均导热系数，δtj为最大允许计算时间步长，ρ是坐标为(x,y,z)的单元的密度，c是坐标为(x,y,z)的单元的比热容，δx、δy、δz分别表示空间内在x方向上相邻、y方向上相邻、z方向上相邻的任意2个单元中心在空间内沿x方向、y方向和z方向的距离，这与单元在x方向、y方向、z方向上的单元尺寸相同。

进一步的，步骤3中，空间内的某一单元a单元的最大允许计算时间步长δtj满足下式：

上式中，j＝1,2...，δtj是最大允许计算时间步长，ρ是a单元材料的密度，c是a单元材料的比热容，λa是a单元材料的导热系数，λb是b单元材料的导热系数，λc是c单元材料的导热系数，λd是d单元材料的导热系数，λe是e单元材料的导热系数，λf是f单元材料的导热系数，λg是g单元材料的导热系数，δxj、δyj、δzj是a单元的尺寸；其中b单元和c单元为在x方向上与a单元相邻的单元，其中d单元和e单元为在y方向上与a单元相邻的单元，其中z单元和g单元为在x方向上与a单元相邻的单元。

进一步的，步骤4中，定义δtv＝δx/v，v是激光扫描速度，δtv是激光在es单元之间的移动时间；em-em单元之间以及es-es单元之间的传热时间均为δt1，el-el单元之间的传热时间为δt2；

在速度v逐渐增大过程中，各时间步长存在三种关系:δt1<δt2<δtv、δt1<δtv<δt2或δtv<δt1<δt2；

当δt1<δt2<δtv成立时：

定义n1＝δtv/δt2，进行δt2数值修正，使n1取正整数，修正方法为将n1向大值取整,δt2＝δtv/n1；定义n2＝δt2/δt1，进行δt1值修正，使n2取正整数，修正方法为将n2向大值取整，δt1＝δt2/n2；

当δt1<δtv<δt2成立时：

定义m1＝δt2/δtv，进行δt2数值修正，使m1取正整数，修正方法为将m1向小值取整，δt2＝m1·δtv；定义m2＝δtv/δt1，进行δt1数值修正，使m2取正整数，修正方法为将m2向大值取整，δt1＝δtv/m2；

当δtv<δt1<δt2成立时：

保持δtv不变，减小δt1使得δtv＝δt1，n1＝δtv/δt2，进行δt2数值修正，使n1取正整数，修正方法为将n1向大值取整,δt2＝δtv/n1。

与现有技术相比，本发明至少具有以下有益的技术效果：

1)基于有限差分原理模拟slm成形过程，简化了建模难度；

2)考虑材料的可变热物性参数(热导率、密度和比热容)使得计算结果更加准确；

3)大幅减少了模拟所需计算机内存和存储空间；在保证计算精度的同时加快了计算速度，使大尺寸样件模拟成为可能。

本方法基于有限差分原理，采用单元嵌套方法进行温度场模拟，突破了有限元方法、其他有限差分方法模拟尺寸小和计算时间长的限制。相较于有限体积法，本方法仅在温度梯度大的熔池周围采用小尺寸离散单元，而在远离熔池且温度梯度小的区域使用大尺寸离散单元，使用单元嵌套的方法缩减了参与计算的离散单元数目。更进一步的，通过计算时间步长的嵌套，能够减少单元与单元之间的计算次数，从而能够减小slm温度场模拟计算的计算量，且能保持较高的计算精度。为大尺寸样件的温度场数值模拟提供新的解决办法。

进一步的，步骤2中，在成形过程中，细化分割区域跟随激光移动而改变，始终保证熔池附近的热量传递以小尺寸的es单元计算。

进一步的，步骤2中，δz尺寸与实际打印铺粉层厚相同，以便模拟铺粉行为。

附图说明

图1是本发明的流程图；

图2是各种类型单元从属关系图；

图3是传热单元接触示意图；

图4是温度监测点位置示意图；

图5是示例温度模拟结果图。

附图中：1-el单元，2-em单元，3-es单元，4-温度监测点，5-样件平面。

具体实施方式

为了使本发明的目的和技术方案更加清晰和便于理解。以下结合附图和实施例，对本发明进行进一步的详细说明，此处所描述的具体实施例仅用于解释本发明，并非用于限定本发明。

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

参照图1，一种slm成形过程温度场数值模拟方法，包括以下步骤：

步骤1，建立仿真模型。在三维建模软件中建立基板和样件的三维模型，对样件和基板分别赋予相应的材料属性——热导率、密度和比热容。

步骤2，离散有限差分单元。将步骤1得到的基板和样件的三维模型均按照一定尺寸离散得到若干正六面体有限差分单元el单元(el单元尺寸为热源激光光斑直径的2-10倍)。将熔池附近区域的el单元细化分割得到比el单元尺寸小的em单元，将连同激光所处el单元在内的4×2×2个el单元1(x方向数量为4、y方向数量为2、z方向数量为2)所包含的em单元再次细化分割为若干es单元。el单元、em单元和es单元组成层层嵌套的“嵌套单元”。通过使用嵌套单元的方法，使得熔池附近区域的传热计算以尺寸相对较小的es单元和em单元进行，而在其他区域传热计算使用尺寸较大的el单元进行。

步骤3，建立传热计算模型，确定计算时间步长。基板、成形件和粉末三者之间的传热主要通过接触热传导进行。考虑slm成形过程中热源形式以及铺粉行为，基于有限差分原理，建立同类型单元之间的换热方程。根据步骤2确定的el单元、em单元及em单元的尺寸，结合基板材料、样件材料(以及粉末)的热物性参数计算同类型单元之间传热的最大允许计算时间步长以保证传热计算收敛。

步骤4，确定传热策略。根据激光的扫描速度大小确定激光在步骤2所得es单元之间的移动时间，并比较其与步骤3得到的各单元之间的最大允许计算时间步长的大小关系，确定传热计算策略，并对计算时间步长进行修正。

步骤5，计算温度场。按照步骤4确定的传热计算策略开始模拟slm成形过程，即可得到slm成形过程温度场的数值模拟结果，计算结果包含整个成形过程中各类型单元的温度变化历程，也可根据研究需求可提取具体单元的温度变化历程。

在所述步骤2中，将基板的三维模型和样件的三维模型按照长为δx，宽为δy，高为δz的正六面体el单元离散，el单元各方向尺寸为热源激光光斑直径的2-10倍。定义em单元是长为δx，宽为δy，高为δz的正六面体单元。定义es单元是长为δx，宽为δy，高为δz的正六面体单元，δx与δy尺寸应与热源激光光斑直径相同，δz尺寸与实际打印铺粉层厚相同，以便模拟铺粉行为。定义n是正六面体el单元与正六面体es单元的尺寸比，即n始终为正整数。其中δx、δy、δz是el单元的边长尺寸，δx、δy、δz是es单元的边长尺寸。定义在模拟slm成形时，将所有z方向坐标等于激光所在el单元坐标zc以及等于zc-1的el单元1细化分割为若干em单元2。应当理解,在z方向上可以将更多层el单元细化分割，但细化分割单元数越多，参与传热计算的单元数目越多，将影响计算速度。el单元细化分割层数应为2-5层。将在激光形成的熔池附近的el单元1——连同激光所处el单元在内的4×2×2个el单元1(x方向数量为4、y方向数量为2、z方向数量为2)所包含的em2单元再次细化分割为若干es单元3。应当理解，“熔池附近”所指区域亦可以增大或减小，此区域越大，则计算量越大，计算结果越可靠；此区域越小，则计算量越小，计算速度越快。熔池附近区域至少应包含2×2×2个el单元,该区域包含的el单元数在x方向、y方向上el单元不能超过样件模型离散得到的x方向el单元数、y方向el单元数。上述三种单元及其嵌套关系如图2所示。在成形过程中，细化分割区域跟随激光移动而改变，始终保证熔池附近的热量传递以小尺寸的es单元3计算。

在所述步骤3中，固体的导热系数远大于气体，主要考虑slm成形过程中基板内部、基板与成形件之间、基板与粉末之间、成形件内部、成形件与粉末之间、粉末内部的热传导效应。slm成形过程中的空间内任一单元位置用其几何中心的点坐标表示，其与周围单元的有限差分换热公式可表示为：

上式中分别是i、i+1时刻坐标为(x,y,z)的单元的温度，是i时刻坐标为(x+1,y,z)的单元的温度，是i时刻坐标为(x-1,y,z)的单元的温度，是i时刻坐标为(x,y+1,z)的单元的温度，是i时刻坐标为(x,y-1,z)的单元的温度，是i时刻坐标为(x,y,z+1)的单元的温度，是i时刻坐标为(x,y,z-1)的单元的温度。λ1是坐标为(x,y,z)的单元与坐标为(x+1,y,z)的单元材料的平均导热系数，λ2是坐标为(x,y,z)的单元与坐标为(x-1,y,z)的单元材料的平均导热系数，λ3是坐标为(x,y,z)的单元与坐标为(x,y+1,z)的单元材料的平均导热系数，λ4是坐标为(x,y,z)的单元与坐标为(x,y-1,z)的单元材料的平均导热系数，λ5是坐标为(x,y,z)的单元与坐标为(x,y,z+1)的单元材料的平均导热系数，λ6是坐标为(x,y,z)的单元与坐标为(x,y,z+1)的单元材料的平均导热系数，δtj是计算时间步长，ρ是坐标为(x,y,z)的单元的密度，c是坐标为(x,y,z)的单元的比热容，δx、δy、δz分别表示空间内在x方向上相邻、y方向上相邻、z方向上相邻的任意2个单元中心在空间内沿x方向、y方向和z方向的距离，这与单元在x方向、y方向、z方向上的单元尺寸相同。

为保证显式有限差分方程计算收敛，如图3所示，空间内任意尺寸的单元a的计算时间步长δtj应满足式(2)：

上式中，j＝1,2...。δtj是计算时间步长，ρ是a单元材料的密度，c是a单元材料的比热容，λa是a单元材料的导热系数，λb是b单元材料的导热系数，λc是c单元材料的导热系数，λd是d单元材料的导热系数，λe是e单元材料的导热系数，λf是f单元材料的导热系数，λg是g单元材料的导热系数，δxj、δyj、δzj是a单元的尺寸。

计算中所采用的时间步长δtj小于或等于由式(2)得到的极大值。材料热物性参数受温度影响，如图3所示的a单元及其周围单位温度发生改变时，会影响计算时间步长的极值。通过计算a单元及其周围单元温度的多种组合(如a单元为50℃，相邻单元为30℃、a单元为100℃，相邻单元为200℃等)，即可得到多个最大允许时间步长。应取尽可能多的各种温度组合所确定的δtj最小值为最大允许时间步长，以使得所有可能的温度组合下传热计算都收敛。结合单元尺寸，根据式(2)计算es-es单元之间的最大允许计算传热时间步长δt1，计算el-el单元之间的最大允许计算传热时间步长δt2。

下文称最大允许计算传热时间步长δt1为计算时间步长δt1，称最大允许计算传热时间步长δt2为计算时间步长δt2。

在所述步骤4中，定义δtv＝δx/v，v是激光扫描速度，δtv是激光在步骤2所得es单元之间的移动时间。在速度v逐渐增大过程中，各时间步长存在三种关系:

①δt1<δt2<δtv②δt1<δtv<δt2③δtv<δt1<δt2

关系①成立时，定义n1＝δtv/δt2。进行δt2数值修正，使n1取正整数，修正方法为将n1向大值取整,δt2＝δtv/n1；定义n2＝δt2/δt1，进行δt1值修正，使n2取正整数，修正方法为将n2向大值取整，δt1＝δt2/n2。修正方法能确保计算时间步长均满足式(2)。经过数值修正，相邻的es-es单元间、相邻的em-em单元间传热计算n2次后，进行1次相邻的el-el单元间传热计算，el-el单元间传热计算n1次后，模拟激光光源移动1次。

当关系②成立时，定义m1＝δt2/δtv。进行δt2数值修正，使m1取正整数，修正方法为将m1向小值取整，δt2＝m1·δtv。定义m2＝δtv/δt1。进行δt1数值修正，使m2取正整数，修正方法为将m2向大值取整，δt1＝δtv/m2。修正方法能确保计算时间步长均满足式(2)。相邻的es-es单元间、相邻的em-em单元间传热计算m1次后，模拟光源移动进行1次；模拟光源移动m2次后，相邻的el-el单元间传热计算1次。

当关系③成立时，为完整模拟激光移动过程，捕捉激光位置变动，保持δtv不变，减小δt1使得δtv＝δt1。如前文所述即n1＝δtv/δt2，进行δt2数值修正，使n1取正整数，修正方法为将n1向大值取整,δt2＝δtv/n1。经过es-es单元间、em-em单元间传热计算1次后，模拟激光移动1次；模拟激光移动n1次后，进行1次el-el单元间传热计算。

参与计算时，遵循能量守恒原则，在本示例中每个el单元温度为其所包含所有em单元的平均值，每个em单元温度为其所包含所有es单元的平均值。

实施例1

考虑基于本方法开发一种计算软件实现温度场模拟功能。

(1)建立尺寸为20mm*20mm*10mm的基板模型，建立尺寸为10mm*10mm*1.5mm的样件模型，基板材料为316l不锈钢，成形粉末为316l不锈钢粉末。

(2)将基板模型和样件模型按照1mm*1mm*0.5mm的el单元离散，共得到9200个el单元，9200个el单元包含了在slm铺粉过程中形成的粉末单元。其中基板始终存在，样件随成形过程逐渐生成。将正在成形生成的样件层el单元及其下层el单元细化分割得到尺寸为1mm*1mm*0.05mm的8000个片层型em单元，使得在熔池附近将4×2×2个el单元包含的所有em单元细化分割得到16000个尺寸为0.1mm*0.1mm*0.05mm的es单元，此时的铺粉层厚为0.05mm。

而若模型全部采用0.1mm*0.1mm*0.05mm尺寸进行离散，将得到9200000个es单元。与之相较，本方法大幅缩减了数据量，模型尺寸越大，其优化作用越明显。

(3)某单元与其周围单元的换热计算公式为：

根据公式

考虑多种温度下各种材料的比热容后计算得到es-es单元之间的最大传热时间步长为δt1＝5.86242×10^-5s，em-em单元之间传热时间同为δt1。el-el单元之间的传热时间为δt2＝5.86242×10^-3。

(4)设定激光的扫描速度为3500mm/s，则δtv＝2.85714×10^-5，此时满足δtv<δt1<δt2。

(5)保持δtv不变，δt1向小值修正，则有δt1＝δtv＝2.85714×10^-5s。定义为充分保证计算收敛，δt2向小值修正，取n2＝204，δt2＝δt1×204＝5.8286×10^-3s。即传热计算策略为，成形过程每经过2.85714×10^-5s，es-es单元之间进行一次传热计算，随后立即进行em-em单元间传热计算，并进行一次热源模拟移动。es-es单元间计算204次后，所有el-el单元之间进行传热计算一次。

需要说明的是，参与计算时，遵循能量守恒原则，在本示例中每个el单元温度为其所包含所有em单元的平均值，每个em单元温度为其所包含所有es单元的平均值。但此方式并不用于限定本方法。当激光移动次数等于边长之比n时，需要移动细化分割区域以保证在熔池附近始终以小尺寸网格进行传热计算。在本示例中，在空间三个方向均设定激光移动n＝10次后重新规划细化区域。

采用本方法所述计算传热方法，使在远离熔池的温度变化缓慢区域以较长时间步长计算传热，在减少所需计算单元数目的基础上，更进一步减小了计算量。需要注意的是，各单元尺寸相差越大，计算量优化效果越好，但尺寸差不宜过大，以免影响计算精度。其中es单元尺寸应尽量与激光光斑大小相同。

(6)开始传热计算。定义一个温度监测点4来描述在特定点在slm成形过程中的温度模拟数据，其位于样件平面5的几何中心，如图4所示。图5给出了该点的温度模拟结果。slm成形过程中，熔池是否稳定是成形关键因素。使用本方法能正确模拟当热源靠近该点(即该点处于熔池附近)时的温度变化。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。