CN115455776B - 增材制造残余应力与变形预测的基于物理的固有应变法 - Google Patents

Abstract

本发明公开了一种增材制造残余应力与变形预测的基于物理的固有应变法，属于增材制造技术领域。包括：提出了更能反映增材制造物理过程的固有应变计算方法，将固有应变进行分解，大幅提高了残余应力的计算精度；提出了附加塑性应变的概念和估算方法，当进行多熔覆层合并计算时，大幅消减了合并层之间不真实的应力震荡；提出了将塑性应变分解为局部相关和结构相关——两部分的思路，局部相关塑性应变反映局部热循环及局部约束的影响，而结构相关塑性应变反映热源路径及结构形状、尺寸的影响。本发明可以广泛应用于工业级大尺寸金属构件的增材过程残余应力与变形模拟，在大幅提高计算效率、保证变形计算精度的同时，得到合理的残余应力计算结果。

Description

增材制造残余应力与变形预测的基于物理的固有应变法

技术领域

本发明属于增材制造技术领域，更具体地，涉及一种增材制造残余应力与变形预测的基于物理的固有应变法。

背景技术

金属增材制造技术是以金属粉末或焊丝为原料，通过激光或电弧等热源熔化、快速凝固逐层沉积“生长制造”，基于构件CAD模型一步完成全致密、高性能金属构件的“近净成形制造”。其一经提出就成为受到广泛关注的“高性能材料制备与复杂构件直接近净成形一体化”新技术，在航空、航天、兵器、能源、化工、机械和医疗等领域有广阔的应用前景。

但现阶段，仍存在制约大型金属增材制造技术广泛工业化应用的技术瓶颈。构件在增材成形过程中经历的激光或电弧热源的长期周期性剧烈加热与冷却、移动熔池在池底强约束下的快速凝固收缩等物理现象，在构件内部产生很大、极其复杂的残余应力以及应力集中，导致构件严重变形与开裂。

由于成形过程中热、力等演化规律极其复杂，几乎不可能求出控制方程的解析解。而如果依靠大量试验、完全凭经验摸索最优的成形控制参数，无法明确其中的热、力学机理，且存在成本高、效率低的问题，特别是航空航天、武器装备等领域中的大型金属构件，单个构件制造成本动辄几十万，甚至几百万，不允许进行“试错式”的研制。因此，需要发展高效的数值模拟方法。

金属增材制造过程是一个强非线性的多物理场耦合问题，大型构件的应力、变形的数值模拟方法主要有“热弹塑性有限元法”及“固有应变法”。但是，这两种方法目前存在明显局限。

热弹塑性有限元法计算精度较高，但是计算效率很低，计算成本很高。中国发明专利CN110348072A公开说明一种提高电弧增材构件热力分析有限元数值计算效率的方法，该方法提出了“模型迭代-分割计算”方法，取得了较好的效果，计算效率较传统方法提高了一倍，但由于没有量级上的提高，因此对于大型构件而言，计算效率仍然远远不够。

固有应变方法计算效率高，但是现有固有应变法由于过强的主观假设，缺乏完善的理论支撑，计算精度和可靠性存在较大问题，尤其是无法合理预测大型构件的残余应力。中国发明专利CN108399280A公开了一种预测激光选区熔化成形件变形的有限元仿真方法，中国发明专利CN113976920A公开了选区激光熔化成形结构残余变形的跨尺度控制方法及系统，二者都是基于固有应变法思想，针对激光选区熔化成形，但主要涉及变形的预测，没有涉及残余应力的预测。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种增材制造残余应力与变形预测的基于物理的固有应变法，其为一种保证计算精度的同时，大幅提高计算效率的金属增材制造残余应力与变形的数值模拟方法。其基于对增材制造物理过程的观察和规律总结，因此被称为基于物理的固有应变方法。本发明提出了更能反映增材制造物理过程的固有应变计算方法，即把每一熔覆层内的固有应变视为多个部分的叠加，分别反映熔覆和多次重熔的影响，大幅提高了残余应力的计算精度；提出了附加塑性应变的概念和估算方法，当进行多熔覆层合并计算时，将附加塑性应变引入到所施加的固有应变后，保证了熔覆层与被覆盖层之间位移的连续性，大幅消减了合并层之间不真实的应力震荡。另外，在固有应变方法的理论基础方面，提出了将塑性应变分解为局部相关和结构相关——两部分的思路，局部相关塑性应变反映局部热循环(熔覆、冷却、重熔)及局部约束的影响，作为固有应变(也即“载荷”)施加到构件上，而结构相关塑性应变反映热源路径及结构形状、尺寸的影响，是“载荷”作用的结果。由此，使固有应变法的物理意义更加明确，即由于局部相关塑性应变只与局部热循环(熔覆、冷却、重熔)及局部约束相关，因而可以通过代表性小尺寸模型求解得到，并可以作为固有应变施加到任意全尺寸构件模型中，以求得构件的变形和残余应力。与传统固有应变法相比，本发明不仅可以高效计算增材制造构件的变形，还可以得到残余应力的合理结果；与热弹塑性有限元法相比，本发明计算效率可以提高两个数量级。

为达到上述目的，本发明采用以下技术方案予以实现：

增材制造残余应力与变形预测的基于物理的固有应变法，包括以下步骤：

步骤(1)、代表性小尺寸模型-提取固有应变：根据目标成形构件的工艺参数、材料参数，建立小尺寸模型，进行热弹塑性有限元模拟，提取某层熔覆时热影响区内各层的局部相关塑性应变，用以计算跟熔覆及各次重熔相对应的局部相关塑性应变，作为真实固有应变；

步骤(2)、全尺寸构件模型-施加固有应变，求解目标成形构件残余应力和变形：根据目标成形构件的结构尺寸图，建立成形构件和基板全尺寸实体模型，对全尺寸模型划分网格，得到网格划分后的全尺寸有限元模型，每次激活N层熔覆层，计算相应的附加塑性应变，叠加真实固有应变后，作为总固有应变，施加到全尺寸模型进行弹塑性计算。

进一步地，所述步骤(1)中，热弹塑性有限元模拟采用顺序热力耦合的方式，隐式算法求解。

进一步地，所述步骤(1)中，熔覆为k层，每层在熔覆过程中最多受到t次重熔，分别提取第k-t层熔覆并冷却到室温后中间区域的平均塑性应变、t次重熔并分别冷却到室温后中间区域的平均塑性应变，作为局部相关塑性应变，其中k>t。

进一步地，所述步骤(2)中，总固有应变等于局部相关塑性应变与附加塑性应变之和，局部相关塑性应变通过代表性小尺寸模型的热弹塑性分析计算得到，而附加塑性应变近似等于局部相关塑性应变与弹性应变的逐层累积。

进一步地，所述步骤(2)中，若N＝1，即采用逐层激活，此时附加塑性应变为零，所施加的固有应变只包含局部相关塑性应变；或者为了进一步提高计算效率，采用多层同时激活，即N>1，此时为了保证激活层与被覆盖层之间位移连续，所施加的固有应变还应包含附加塑性应变。

进一步地，所述步骤(2)中，考虑到材料熔覆后的多次重熔，第m层被激活时，被激活层及其下面附近s层施加固有应变，分别表示第m层的熔覆、第m-1至第m-s层的重熔；激活每一层时，被激活层及其下面附近s层内施加的固有应变分布相同，且s＝t。

进一步地，所述步骤(2)中，利用生死单元法每次激活N层熔覆层，且被激活前沉积件的单元为“死单元”，采用静力学隐式分析步，需要设置的材料参数和材料模型包括常温下的弹性常数、塑性加载函数、硬化模型和流动法则。

本发明与现有技术相比，有益效果为：

本发明基于对增材制造物理过程的观察和规律总结以及理论建模，夯实了固有应变法的理论基础，提出了新的固有应变法的实施过程，使固有应变的求解与施加更加符合增材成形的实际物理过程，进而提高了固有应变法预测残余应力与变形的可靠性和精度。

固有应变求解方面：根据固有应变的产生机制，提出了更能反映增材制造物理过程的固有应变计算方法，即把每一熔覆层内的固有应变视为多个部分的叠加，分别反映熔覆和多次重熔的影响，并且给出了各个应变分量的解析表达式，避免了过强的主观假设，大幅提高了残余应力的计算精度。

固有应变施加方面：考虑到增材过程的逐层熔覆，多次反复重熔、冷却的特点，将固有应变的施加过程进行分解，分别表示材料首次熔覆和后续的重熔，充分考虑到增材的实际物理过程；在进行多熔覆层合并施加固有应变时，引入附加塑性应变，保证了熔覆层与被覆盖层之间位移的连续性，进而大幅消减了合并层之间不真实的应力震荡。

本发明既具有与热弹塑性有限元法相当的求解精度，又具有固有应变法的仿真效率，适用于大尺寸构件变形与残余应力的高精度、高效率仿真。

附图说明

图1为增材过程和逐层施加固有应变的示意图；

图2为合并6层施加真实固有应变的示意图；

图3为本发明方法的流程图；

图4为实施例中的全尺寸构件模型的几何尺寸图；

图5为实施例中的代表性小尺寸模型的几何尺寸图；

图6为实施例中的全尺寸构件模型的热弹塑性有限元法的残余应力计算结果；

图7为实施例中的全尺寸构件模型的基于物理的固有应变法的残余应力计算结果(逐层激活)；

图8为实施例中的全尺寸构件模型的基于物理的固有应变法的残余应力计算结果(合并6层激活)；

图9为实施例中的全尺寸构件模型第20层熔覆层的基于物理的固有应变法与热弹塑性有限元法残余应力计算结果对比；

图10为实施例中的全尺寸构件模型第30层熔覆层的基于物理的固有应变法与热弹塑性有限元法残余应力计算结果对比；

图11为实施例中的全尺寸构件模型第40层熔覆层的基于物理的固有应变法与热弹塑性有限元法残余应力计算结果对比；

图12为实施例中的全尺寸构件模型的基板底面中线变形计算结果对比；

图13为参考文献中模型1的尺寸图；

图14为参考文献中模型2的尺寸图；

图15为基于物理的固有应变法的位移计算结果与参考文献中模型1的实验结果对比；

图16为基于物理的固有应变法的残余应力计算结果与参考文献中模型2的实验结果对比。

具体实施方式

下面结合具体实例对本发明详细说明。结合图3，增材制造残余应力与变形预测的基于物理的固有应变法，包括以下步骤：

步骤(1)、代表性小尺寸模型-提取固有应变。建立代表性小模型的有限元模型，进行热弹塑性有限元模拟，计算结束后，提取固有应变。具体包括：

(1.1)代表性小模型的热弹塑性有限元模拟：

根据目标成形构件的工艺参数、材料参数，建立代表性小尺寸几何模型，划分网格后得到有限元模型。在移动热源作用下，采用隐式算法，进行热力顺序耦合模拟，即通过热分析得到温度场，热分析的结果作为温度载荷施加到力分析模型，得到塑性应变场。总共熔覆k层；计算结束后，提取真实塑性应变。

(1.2)真实塑性应变的分解：

将真实塑性应变分解为“局部相关塑性应变”和“结构相关塑性应变”，即：

其中，为真实塑性应变；/>为局部相关塑性应变，反映了熔覆与重熔的影响，与结构形状、尺寸无关；/>为结构相关塑性应变；X⁽ⁱ⁾＝{X₁,X₂,X₃}为全局坐标系下第i熔覆层材料点空间坐标，如图4所示，1方向为热源移动方向，2方向为熔覆层宽度方向，3方向为沉积件高度方向。

(1.3)局部相关塑性应变的提取：

将局部相关塑性应变进一步分解为首次熔覆造成的塑性应变与后续t次重熔造成的塑性应变；熔覆k层，分别计算第k-t层熔覆并冷却到室温后中间区域的平均塑性应变、t次重熔并分别冷却到室温后中间区域的平均塑性应变，即局部相关塑性应变，其中k>t，计算公式如下：

其中，x⁰为层内局部坐标，并且全局坐标X⁽ⁱ⁾＝x⁰+{0,0,(i-1)h}，其中h为熔覆层的厚度；表示首次熔覆造成的塑性应变，/>(l＝1...t)表示第l次重熔造成的塑性应变。

将此局部相关塑性应变作为总固有应变的一部分，作为全尺寸构件模型的输入参数：“载荷”。

(1.4)结构相关塑性应变：

结构相关塑性应变反映了结构形状、尺寸的影响，视为总固有应变“载荷”的作用结果。

优选的，步骤(1)中，取k＝6。

优选的，步骤(1)中，取t＝4。重熔总次数t的取值与具体增材形式、工艺参数、材料等有关，可以根据热源的作用深度确定其取值。如图1所示，在移动热源作用下，材料逐层熔覆，热源作用深度为5层熔覆层厚度，即每层熔覆层沉积之后最多会受到4次重熔作用。

步骤(2)、全尺寸构件模型-施加固有应变。建立全尺寸构件有限元模型，设定等效合并层数，逐层激活等效合并层，施加固有应变，进行弹塑性求解，得到残余应力和变形。具体包括：

(2.1)根据目标成形构件的结构尺寸图，建立全尺寸构件几何模型，划分网格后得到全尺寸有限元模型，需要设置的材料参数和材料模型包括常温下的弹性常数、塑性加载函数、硬化模型和流动法则。设定等效合并层数r，每个等效合并层包含N层实际熔覆层，N为大于等于1的自然数。

(2.2)若N＝1，采用逐层激活的方式：

将“局部相关塑性应变”作为固有应变。如图1所示(其中IS0表示熔覆所造成的真实固有应变，IS1～IS4表示重熔造成的真实固有应变)，第m层被激活时，被激活层及其下面附近s层施加固有应变，分别表示第m层的熔覆、第m-1至第m-s层的重熔；激活每一层时，被激活层及其下面附近s层内施加的固有应变分布相同，其中m为非零自然数。激活第m层时，对第i层施加的固有应变为：

其中，s＝t。

(2.3)若N>1，采用合并多层激活的方式：

将“局部相关塑性应变”叠加“附加塑性应变”作为总固有应变，即：

其中，ε^tin(X⁽ⁱ⁾,r)，ε^trin(X⁽ⁱ⁾,r)和ε^ain(X⁽ⁱ⁾,r)分别为第r个等效合并层中的第i熔覆层的总固有应变，由局部相关塑性应变构成的真实固有应变，以及附加固有(塑性)应变。

真实固有应变的计算公式为：

施加真实固有应变的示意图如图2所示(其中IS0表示熔覆所造成的真实固有应变，IS1～IS4表示重熔造成的真实固有应变)。

附加固有(塑性)应变的计算公式为：

其中，括号<·>表示保留面内分量，其余分量为零；热源移动方向弹性应变分量熔覆道宽度方向弹性应变分量/>其余弹性应变分量为零，σ_s是常温屈服极限，E和μ分是常温杨氏模量和泊松比。

(2.4)采用静力学隐式分析步，利用生死单元法逐层激活等效合并层(每次激活N层熔覆层)，且被激活前沉积件的单元为“死单元”，进行弹塑性求解，得到残余应力、变形。

优选的，步骤(2)中，取N<10。

实施例：

以电弧熔丝增材316L不锈钢60层“高墙”模型为例说明本发明的具体实施方法。分别采用高精度的热弹塑性有限元法和基于物理的固有应变法进行模拟。

全尺寸构件几何模型如图4所示。基板几何模型尺寸为250mm*100mm*20mm，沉积件几何模型的尺寸为180mm*10mm*90mm，总共沉积60层。由于模型的对称性，为了提高计算速度，采用半模型。

代表性小尺寸几何模型如图5所示。基板几何模型尺寸为140mm*50mm*12mm，沉积件几何模型的尺寸为90mm*10mm*9mm，总共沉积6层。同样采用半模型进行计算。

高精度的热弹塑性有限元法计算得到的全尺寸模型应力结果如图6所示。

基于物理的固有应变法，若采用如图1所示的逐层激活的方式，全尺寸模型的应力计算结果如图7所示；若采用如图2所示的合并6层激活的方式，全尺寸模型的应力计算结果如图8所示。

图9，图10，图11分别为全尺寸模型沉积件第20、30、40层熔覆层处，热弹塑性有限元法和基于物理的固有应变法(逐层激活、合并6层激活)计算得到的Mises应力结果对比；图12为全尺寸模型基板底面中线3方向变形结果对比。

可见，本发明的变形计算结果与热弹塑性有限元的计算结果吻合很好，应力计算结果能够反映应力的分布规律，并能较好地预测增材制造构件内应力较大的位置，即容易产生破坏的位置。与传统固有应变法相比，本发明不仅可以高效计算增材制造构件的变形，还可以得到残余应力的合理结果。

计算时长如下表所示：

可见，基于物理的固有应变法在保证残余应力、变形计算精度的同时，大大提高了仿真效率。与热弹塑性有限元法相比，本发明计算效率可以提高两个数量级。而且，因为通过小尺寸模型得到的固有应变可用于相同材料和工艺参数的任意大尺寸构件的变形与残余应力预测，所以构件尺寸越大，效率提高会越明显。

图13，图14，图15，图16为本发明的基于物理的固有应变方法与文献(Ding J,Colegrove P,Mehnen J,et al.Thermo-mechanical analysis of Wire and ArcAdditive Layer Manufacturing process on large multi-layer parts[J].Computational Materials Science,2011,50(12):3315-3322.)中的实验结果对比。图13为论文中测量位移的构件，其基板尺寸为500mm*60mm*12mm，沉积件尺寸为500mm*5mm*8mm，包含4层熔覆层。打印结束后释放基板约束，测量基板底面中线的位移。图14为论文中测量残余应力的构件，其基板尺寸为500mm*60mm*12mm，沉积件尺寸为500mm*5mm*40mm，包含20层熔覆层。打印结束后释放基板约束，测量沉积件中间的纵向残余应力。基于物理的固有应变法的位移预测结果与实验结果的对比如图15所示；基于物理的固有应变法的残余应力预测结果与实验结果的对比如图16所示。可见，本发明的残余应力、位移的预测结果与实验结果吻合较好。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种增材制造残余应力与变形预测的基于物理的固有应变法，其特征在于，包括以下步骤：

步骤(1)、代表性小尺寸模型-提取固有应变：根据目标成形构件的工艺参数、材料参数，建立小尺寸模型，进行热弹塑性有限元模拟，提取某层熔覆时热影响区内各层的局部相关塑性应变，用以计算跟熔覆及各次重熔相对应的局部相关塑性应变，作为真实固有应变；热弹塑性有限元模拟采用顺序热力耦合的方式，隐式算法求解；

步骤(2)、全尺寸构件模型-施加固有应变，求解目标成形构件残余应力和变形：根据目标成形构件的结构尺寸图，建立成形构件和基板全尺寸实体模型，对全尺寸模型划分网格，得到网格划分后的全尺寸有限元模型，每次激活N层熔覆层，计算相应的附加塑性应变，叠加真实固有应变后，作为总固有应变，施加到全尺寸模型进行弹塑性计算；

所述步骤(2)中，若N＝1，即采用逐层激活，此时附加塑性应变为零，所施加的固有应变只包含局部相关塑性应变；或者为了进一步提高计算效率，采用多层同时激活，即N>1，此时为了保证激活层与被覆盖层之间位移连续，所施加的固有应变还应包含附加塑性应变；利用生死单元法每次激活N层熔覆层，且被激活前沉积件的单元为“死单元”，采用静力学隐式分析步，需要设置的材料参数和材料模型包括常温下的弹性常数、塑性加载函数、硬化模型和流动法则。

2.根据权利要求1所述的增材制造残余应力与变形预测的基于物理的固有应变法，其特征在于，所述步骤(1)中，熔覆为k层，每层在熔覆过程中最多受到t次重熔，分别提取第k–t层熔覆并冷却到室温后中间区域的平均塑性应变、t次重熔并分别冷却到室温后中间区域的平均塑性应变，作为局部相关塑性应变，其中k>t。

3.根据权利要求1所述的增材制造残余应力与变形预测的基于物理的固有应变法，其特征在于，所述步骤(2)中，总固有应变等于局部相关塑性应变与附加塑性应变之和，局部相关塑性应变通过代表性小尺寸模型的热弹塑性分析计算得到，而附加塑性应变近似等于局部相关塑性应变与弹性应变的逐层累积。

4.根据权利要求1所述的增材制造残余应力与变形预测的基于物理的固有应变法，其特征在于，所述步骤(2)中，考虑到材料熔覆后的多次重熔，第m层被激活时，被激活层及其下面附近s层施加固有应变，分别表示第m层的熔覆、第m-1至第m–s层的重熔；激活每一层时，被激活层及其下面附近s层内施加的固有应变分布相同，且s＝t。