CN112115634A

CN112115634A - 金属液单向凝固过程中晶粒组织三维数值预测方法

Info

Publication number: CN112115634A
Application number: CN202010995830.8A
Authority: CN
Inventors: 刘东戎; 浦震鹏; 赵思聪
Original assignee: Harbin University of Science and Technology
Current assignee: Harbin University of Science and Technology
Priority date: 2020-09-21
Filing date: 2020-09-21
Publication date: 2020-12-22
Anticipated expiration: 2040-09-21
Also published as: CN112115634B

Abstract

一种金属液单向凝固过程中晶粒组织三维数值预测方法，属于金属材料加工领域。本发明针对采用三维元胞自动机模型模拟金属液凝固过程并对晶粒组织进行三维表征，目前已有模拟计算方法存在三维数值预测计算量大，并且未考虑固液界面能对晶粒生长的影响的问题。本发明包括：对计算域划分网格并确定网格标识；对每个网格赋予状态变量；基于每个剖分步长在所有时间步长内计算获得的网格欧拉角计算结果形成一个剖分步长确认文件；计算获得满足要求的最大的剖分步长；针对实际铸件进行凝固过程的预测，将每个时间步长内具有相同欧拉角的网格作为同一晶粒。本发明提高了实际凝固过程数值预测的计算速度。

Description

金属液单向凝固过程中晶粒组织三维数值预测方法

技术领域

本发明涉及金属液单向凝固过程中晶粒组织三维数值预测方法，属于金属材料加工领域。

背景技术

金属液的凝固产生于型腔，因为铸型相比于金属液温度低，可有效将金属液的热量传导出去，因而可促进金属液凝固。金属液在型腔中发生的凝固具有多方向散热特点，因此通常形成具有较大晶界取向差的柱状晶组织。具有较大晶界取向差的柱状晶组织会引起铸件力学性能的不均匀分布，在受到单向载荷时，某些与载荷垂直的晶界极易成为裂纹源。

单向凝固技术可避免形成具有较大晶界取向差的柱状晶组织。其特点是型腔内没有多向散热，而是建立某一特定方向的温度梯度，高温金属液只沿着一个方向散热，使柱状晶粒之间具有相近的结晶取向，因此形成具有较小晶界取向差的柱状晶组织。具有较小晶界取向差的柱状晶组织没有横向晶界，纵向晶界之间基本相互平行，力学性能分布均匀，使其在经受单向载荷时晶界不容易被撕裂，这大大提高了铸件的抗热冲击性和抗蠕变性，从而延长了铸件的疲劳寿命。单向凝固技术目前已被广泛应用于飞机发动机叶片的制造过程。

在单向凝固技术中，控制其凝固过程的组织形成有助于提高铸件力学性能。可通过竞争生长尽量减小柱状晶晶界夹角，例如，柱状晶晶界夹角小于15°要好于小于25°；另外，单晶组织的力学性能优于多晶组织的力学性能。然而，尽管通过大量的实验研究可以建立工艺参数与凝固组织之间的关系，却无法明确工艺参数变化对组织形成机理的影响。这是因为很难借助实验方法测量或跟踪凝固组织演化的整个冶金过程。而开展实验需要购买设备、原材料，不仅耗费大量的人力、物力和财力，同时对能源也是一种浪费。况且基于实验方法研究控制凝固组织必将大大延长工艺优化以及新工艺开发周期，不利于绿色环保和节能减排。与此同时为了提供关于凝固组织的详尽信息，还需要对实验所得晶粒组织进行三维表征。目前开展三维组织实验表征工作的难度很大，主要原因就是工作量巨大导致研究周期太长。相比于实验研究，基于对数学方程的求解，进行数值模拟可以精确描述凝固过程中所涉及的复杂物理现象，研究不同尺度物理现象之间的交互，直观且动态再现凝固过程中的组织形成，并捕捉到在实验中无法提取的物理量。因此，通过数值模拟手段建立工艺参数与凝固组织之间的关系，有助于缩短研制周期、提高工艺设计水平和产品出品率，更重要的是费用低且节能环保。

目前，对金属液单向凝固过程中晶粒组织的数值预测多采用元胞自动机模型进行二维模拟。二维模拟只能再现一个平面上的晶界分布特点，而实际凝固组织由于结晶取向不同，竞争生长在不同维度上特点也不同，因此很难提供全面的晶粒组织形貌信息。另外，有少数数值预测算法采用三维元胞自动机模型借助偏心算法模拟凝固组织，因为在计算过程中需要不断计算空间某一点与空间八面体之间的距离关系，来判断该点可否被空间八面体捕获；因此，偏心算法的计算量大；同时偏心算法还忽略了固液界面能各向异性对晶粒生长的影响。

发明内容

针对采用三维元胞自动机模型模拟金属液凝固过程并对晶粒组织进行三维表征，目前已有模拟计算方法存在三维数值预测计算量大，并且未考虑固液界面能对晶粒生长的影响的问题，本发明提供一种金属液单向凝固过程中晶粒组织三维数值预测方法。

本发明的一种金属液单向凝固过程中晶粒组织三维数值预测方法，包括：

步骤一：对假定铸件单向凝固过程的计算域进行设定，在常规三维坐标系下选择所述计算域在各个坐标轴方向的剖分步长Δx_CA，根据剖分步长Δx_CA对计算域划分网格并确定网格标识；对每个网格基于所述网格标识赋予状态变量；

步骤二：由小到大依次选择不同大小的剖分步长Δx_CA，在设定的温度梯度和冷却速度下，结合状态变量计算金属液单向凝固过程中直至凝固结束时每个网格对应的三个欧拉角；将每个时间步长内具有相同欧拉角的网格确定为同一晶粒；基于每个剖分步长Δx_CA在所有时间步长内计算获得的网格欧拉角计算结果形成一个剖分步长确认文件grain.txt；

步骤三：对由不同剖分步长Δx_CA计算获得的剖分步长确认文件grain.txt文件进行分析，提取晶界夹角值；若不同剖分步长Δx_CA对应的晶界夹角差值的绝对值大于或等于5°,则对已选择剖分步长Δx_CA中的最大值进行调整，再返回步骤一，直至不同剖分步长Δx_CA对应的晶界夹角差值的绝对值小于5°，选择最大的剖分步长Δx_CA；

步骤四：针对实际铸件进行凝固过程的预测：基于最大的剖分步长Δx_CA划分网格，采用网格标识对每个网格赋予状态变量；再在设定的温度梯度和冷却速度下，计算金属液单向凝固过程中直至凝固结束时每个时间步长内每个网格对应的三个欧拉角，将每个时间步长内具有相同欧拉角的网格作为同一晶粒；实现实际铸件凝固过程的晶粒组织三维数值预测。

本发明的有益效果：本发明方法以晶粒轮廓线(而非枝晶结构)为主要计算对象，在晶粒生长速度计算中引入固液界面能各向异性因子，从而模拟具有不同生长取向晶粒的生长。为了确保计算结果不依赖于剖分网格尺寸，针对每一合金体系都需要对剖分网格尺寸进行预先检验，确保计算结果在一定程度上不受网格尺寸影响。本发明提高了实际凝固过程数值预测的计算速度，可缩短工艺研发周期，加快组织以及最终性能调控进度，对工业生产具有重要指导意义。

附图说明

图1是本发明所述金属液单向凝固过程中晶粒组织三维数值预测方法的流程图；

图2是具体实施例中工艺参数一下Δx_CA＝2.5μm时计算结果对网格尺寸敏感性示意图；

图3是具体实施例中工艺参数一下Δx_CA＝5μm时计算结果对网格尺寸敏感性示意图；

图4是具体实施例中工艺参数二下柱状晶的欧拉角分别为(0°,10°,0°)和(45°,45°,45°)时Δx_CA＝2.5μm时计算结果对网格尺寸敏感性示意图；

图5是具体实施例中工艺参数二下柱状晶的欧拉角分别为(0°,10°,0°)和(45°,45°,45°)时Δx_CA＝5μm时计算结果对网格尺寸敏感性示意图；

图6是具体实施例中工艺参数二下柱状晶的欧拉角分别为(10°,10°,10°)和(45°,45°,45°)时Δx_CA＝2.5μm时计算结果对网格尺寸敏感性示意图；

图7是具体实施例中工艺参数二下柱状晶的欧拉角分别为(10°,10°,10°)和(45°,45°,45°)时Δx_CA＝5μm时计算结果对网格尺寸敏感性示意图；

图8是具体实施例中工艺参数二下柱状晶的欧拉角分别为(10°,10°,10°)和(25°,25°,25°)时Δx_CA＝2.5μm时计算结果对网格尺寸敏感性示意图；

图9是具体实施例中工艺参数二下柱状晶的欧拉角分别为(10°,10°,10°)和(25°,25°,25°)时Δx_CA＝5μm时计算结果对网格尺寸敏感性示意图；

图10是具体实施例中工艺参数二下柱状晶的欧拉角分别为(10°,10°,10°)和(5°,5°,5°)时Δx_CA＝2.5μm时计算结果对网格尺寸敏感性示意图；

图11是具体实施例中工艺参数二下柱状晶的欧拉角分别为(10°,10°,10°)和(5°,5°,5°)时Δx_CA＝5μm时计算结果对网格尺寸敏感性示意图；

图12是对应于工艺参数二时模拟所得单向凝固过程中晶粒组织；

图13是对应于工艺参数三时模拟所得单向凝固过程中晶粒组织；

图14是针对工艺参数二，模拟所得柱状晶粒组织的三维展示；

图15是图14铸件中垂直于X轴的中心截面示意图；

图16是图14铸件中垂直于Y轴的中心截面示意图；

图17是针对工艺参数三，模拟所得等轴晶粒组织的三维展示；

图18是图17铸件中垂直于X轴的中心截面示意图；

图19是图17铸件中垂直于Y轴的中心截面示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其它实施例，都属于本发明保护的范围。

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

下面结合附图和具体实施例对本发明作进一步说明，但不作为本发明的限定。

具体实施方式一、结合图1所示，本发明提供了一种金属液单向凝固过程中晶粒组织三维数值预测方法，包括：

步骤一：对假定铸件单向凝固过程的计算域进行设定，在常规三维坐标系下选择所述计算域在各个坐标轴方向的剖分步长Δx_CA，根据剖分步长Δx_CA对计算域划分网格并确定网格标识；对每个网格基于所述网格标识赋予状态变量；因为是单向凝固过程，所以只针对铸件进行计算。

步骤二：选定剖分步长Δx_CA后计算结果对网格尺寸敏感度检验：

由小到大依次选择不同大小的剖分步长Δx_CA，在设定的温度梯度和冷却速度下，结合状态变量计算金属液单向凝固过程中直至凝固结束时每个网格对应的三个欧拉角；将每个时间步长内具有相同欧拉角的网格确定为同一晶粒；基于每个剖分步长Δx_CA在所有时间步长内计算获得的网格欧拉角计算结果形成一个剖分步长确认文件grain.txt；

进一步，步骤一中，设定所述计算域尺寸为X米×Y米×Z米；计算域在X轴、Y轴、Z轴方向的最小值分别为X_min＝0、Y_min＝0、Z_min＝0(米)，最大值分别为X_max、Y_max、Z_max(米)；在X轴、Y轴、Z轴方向剖分步长相同均为Δx_CA，设定-Z轴为重力方向；

设定网格标识为(j,i,k)，其中j、i和k均为整数；

其中j表示X轴方向上网格的标号，j取值范围为0～(MTcell-1)：

i表示Y轴方向上网格的标号，i取值范围为0～(LTcell-1)：

k表示Z轴方向上网格的标号，k取值范围为0～(NTcell-1)：

且MTcell>100，LTcell＝5；

每一个网格均为一个正方体，边长为Δx_CA。

对每个网格基于所述网格标识赋予状态变量包括：

温度tem(j,i,k)，状态state(j,i,k)，固相分数fs(j,i,k)，一个时间步长内固相分数改变量dfs(j,i,k)，生长速度vtip(j,i,k)，界面曲率kcurve(j,i,k)，界面各向异性对应的X轴方向取向角gangx(j,i,k)，界面各向异性对应的Y轴方向取向角gangy(j,i,k)，界面各向异性对应的Z轴方向取向角gangz(j,i,k)，三个欧拉角ang1(j,i,k)，ang(j,i,k)，ang2(j,i,k)；

在凝固过程的初始时刻，金属液处于液态，设定状态state(j,i,k)＝0，对应固相分数fs(j,i,k)＝0；金属液处于生长态时，设定state(j,i,k)＝1，0<fs(j,i,k)<1；金属液处于凝固态时，金属液处于凝固态时，设定state(j,i,k)＝2，fs(j,i,k)＝1。

再进一步，步骤二中结合状态变量计算金属液单向凝固过程中直至凝固结束时每个网格对应的三个欧拉角ang1(j,i,k)，ang(j,i,k)，ang2(j,i,k)的方法包括：

步骤二一：计算每个网格(j,i,k)的温度tem(j,i,k)，设定温度只在Z轴方向发生变化：

式中T_l为液相线温度(K)，G为温度梯度(K/m)，

为冷却速度(K/s)，time为凝固时间(s)；

d_z为网格(j,i,k)在Z轴方向距离Z＝0平面的距离(m)：

d_z＝k×Δx_CA；

步骤二二：选定两个网格计算对应的温度，并选取对应的欧拉角：

选择两个网格，设定每个网格对应一个形核位置，记为形核网格(j_nu1,i_nu1,k_nu1)和(j_nu2,i_nu2,k_nu2)；

使j_nu1＝MTcell/4，取整数，i_nu1＝LTcell/2，取整数，k_nu1＝0；

使j_nu2＝3×MTcell/4，取整数，i_nu2＝LTcell/2，取整数，k_nu2＝0；

根据温度tem(j,i,k)的计算公式计算形核网格(j_nu1,i_nu1,k_nu1)和(j_nu2,i_nu2,k_nu2)的当前温度；

在所述当前温度小于液相线温度T_l时，使相应形核网格的状态state从0变为1；同时，对每个形核网格对应的欧拉角进行随机选取：欧拉角ang1为[0,360°]之间的一个随机值，欧拉角ang为[0,180°]之间的一个随机值，欧拉角ang2为[0,360°]之间一随机值。

再进一步，步骤二中，确定同一晶粒的过程包括：

步骤二三：对state＝1的形核网格进行界面曲率、界面各向异性对应的取向角、生长速度、固相分数和固相分数改变量的计算；

1)针对每个形核网格采用三个欧拉角构建旋转矩阵A：

其中，

pa1＝cos(ang1)cos(ang2)-sin(ang1)sin(ang2)cos(ang)，

pb1＝sin(ang1)cos(ang2)+cos(ang1)sin(ang2)cos(ang)，

pc1＝sin(ang2)sin(ang)，

pa2＝-cos(ang1)sin(ang2)-sin(ang1)cos(ang2)cos(ang)，

pb2＝-sin(ang1)sin(ang2)+cos(ang1)cos(ang2)cos(ang)，

pc2＝cos(ang2)sin(ang)，

pa3＝sin(ang1)sin(ang)，

pb3＝-cos(ang1)sin(ang)，

pc3＝cos(ang)；

2)同时计算常规三维坐标系下每个形核网格的的固相分数梯度fsx、fsy和fsz：

3)根据旋转矩阵A和三维坐标系下形核网格的的固相分数梯度计算旋转坐标系下的固相分数梯度(fsx)′、(fsy)′和(fsz)′：

(fsx)'＝pa1·fsx+pa2·fsy+pa3·fsz，

(fsy)'＝pb1·fsx+pb2·fsy+pb3·fsz，

(fsz)'＝pc1·fsx+pc2·fsy+pc3·fsz；

4)再计算旋转坐标系下的界面各向异性对应的X轴方向取向角gangx、Y轴方向取向角gangy和Z轴方向取向角gangz：

对于gangx：

若[(fsx)']²+[(fsz)']²＝0，则

若(fsx)'＝0且(fsz)'＞0，则

若(fsx)'＝0且(fsz)'＜0，则

若(fsz)'＝0且(fsx)'＞0，则gangx＝0；

若(fsz)'＝0且(fsx)'＜0，则gangx＝π；

若(fsz)'≠0且(fsx)'≠0，则

对于gangy：

若[(fsy)']²+[(fsz)']²＝0，则

若(fsy)'＝0且(fsz)'＞0，则

若(fsy)'＝0且(fsz)'＜0，则

若(fsz)'＝0且(fsy)'＞0，则gangy＝0；

若(fsz)'＝0且(fsy)'＜0，则gangy＝π；

若(fsz)'≠0且(fsx)'≠0，则

对于gangz：

若[(fsx)']²+[(fsy)']²＝0，则

若(fsx)'＝0且(fsy)'＞0，则

若(fsx)'＝0且(fsy)'＜0，则

若(fsy)'＝0且(fsx)'＞0，则gangz＝0；

若(fsy)'＝0且(fsx)'＜0，则gangz＝π；

若(fsy)'≠0且(fsx)'≠0，则

5)再计算考虑了界面各向异性影响的生长速度：

式中ΔT表示过冷度：ΔT＝T_l-tem；

6)固相分数改变量和固相分数的计算：

对计算域内非边缘区域的网格，每个网格存在26个相邻网格，分为第一层近邻，第二层近邻和第三层近邻；对于网格(j,i,k)，将第一层近邻网格记录为包括fs_left(j-1,i,k)，fs_right(j+1,i,k)，fs_up(j,i,k+1)，fs_down(j,i,k-1)，fs_front(j,i+1,k)，fs_back(j,i-1,k)；第二层近邻网格记录为包括fs_upleft(j-1,i,k+1)，fs_upright(j+1,i,k+1)，fs_downleft(j-1,i,k-1)，fs_downright(j+1,i,k-1)，fs_frontup(j,i+1,k+1)，fs_frontdown(j,i+1,k-1)，fs_frontleft(j-1,i+1,k)，fs_frontright(j+1,i+1,k)，fs_backup(j,i-1,k+1)，fs_backdown(j,i-1,k-1)，fs_backleft(j-1,i-1,k)，fs_backright(j+1,i-1,k)；第三层近邻网格记录为包括fs_frontupleft(j-1,i+1,k+1)，fs_frontupright(j+1,i+1,k+1)，fs_{frontdownleft}(j-1,i+1,k-1)，fs_{frontdownright}(j+1,i+1,k-1)，fs_backupleft(j-1,i-1,k+1)，fs_backupright(j+1,i-1,k+1)，fs_backdownleft(j-1,i-1,k-1)，fs_{backdownright}(j+1,i-1,k-1)；对26个相邻网格通过其状态state确定其形貌因子MFG；在每一个时间步长Δt内，MFG的初始值被赋予为0，然后进行累加赋值；

对于第一层近邻网格，以fs_left(j-1,i,k)为例，若fs_left(j-1,i,k)的state＝2，则其MFG＝MFG+1；若fs_left(j-1,i,k)的state≠2，则其MFG＝MFG+0；针对第一层近邻网格中的其余网格进行同样判断和计算。

对于第二层近邻网格，以fs_upleft(j-1,i,k+1)为例，若fs_upleft(j-1,i,k+1)的state＝2，则其

若fs_upleft(j-1,i,k+1)的state≠2，则MFG＝MFG+0；针对第二层近邻网格中其余网格进行同样判断和计算。

对于第三层近邻网格，以fs_frontupleft(j-1,i+1,k+1)为例，若fs_frontupleft(j-1,i+1,k+1)的state＝2，则其

若fs_frontupleft(j-1,i+1,k+1)的state≠2，则MFG＝MFG+0；针对第三层近邻网格中其余网格进行同样判断和计算。

获得MFG因子后，根据形貌因子MFG计算一个时间步长内相应网格的固相分数改变量和固相分数：

fs(time)＝fs(time-Δt)+dfs，

式中Δt为时间步长(s)，

D_l为溶质在液相中的扩散系数；dfs为一个时间步长内固相分数增量，fs(time)为当前时刻下固相分数，fs(time-Δt)为上一时刻下固相分数；min(MFG/3,1表示

和1比较取最小值；若fs(time)≥1，则表示网格(j,i,k)完成生长，对应的状态state值从1变为2。

再进一步，步骤二中，确定同一晶粒的过程还包括：

步骤二四：通过网格捕捉再现晶粒长大。将state＝2的网格(j,i,k)作为母网格，对母网格周围26个邻居网格的state值进行判断，若某一邻居网格的state＝0，则使该邻居网格的state从0变为1，表明被捕获的该邻居网格不再是液相，将进入生长状态，同时对被捕获的该邻居网格赋予母网格对应的三个欧拉角；

步骤二五：在设定的温度梯度和冷却速度下，逐次选择不同的剖分步长Δx_CA，使Δx_CA依次为step1，step2，step3，……且step1<step2<step3<……，对应执行步骤步骤一至步骤二四，每一次计算直至凝固过程结束，即所有(j,i,NTcell-1)网格的state值均等于2，(j,i,NTcell-1)代表计算域中k取最大值对应的最顶层网格；然后输出每一次计算结束时，每个网格(j,i,k)对应的三个欧拉角，所有网格对应的欧拉角形成剖分步长确认文件grain.txt。具有相同欧拉角的晶粒为同一晶粒。

再进一步，步骤三中，将对应Δx_CA依次为step1，step2，step3，……的剖分步长确认文件记为grain1.txt，grain2.txt，grain3.txt，……，分析步骤二所获得的grain.txt文件，对每个剖分步长确认文件分别提取晶界夹角值θ₁，θ₂，θ₃，……，若不同剖分步长Δx_CA对应的晶界夹角差值的绝对值大于或等于5°,表明欧拉角的计算结果受剖分步长取值的影响超过安全阈值，则减小step1，step2，step3，……中的最大值，重新执行步骤一；直至不同剖分步长Δx_CA对应的晶界夹角差值的绝对值小于5°，表明欧拉角的计算结果受剖分步长取值的影响在安全阈值以内，从step1，step2，step3，……中选择最大值，作为最终剖分步长。若晶界夹角差值的绝对值小于5°,则认为计算结果基本不受所选剖分网格尺寸的影响；所述安全阈值可为根据经验选定的，与晶界夹角差值等于5°时的对应值。

再进一步，步骤四中，采用网格标识对每个网格赋予状态变量包括：

步骤四一：根据实际铸件设定计算域尺寸为XR米×YR米×ZR米，计算域在X轴、Y轴、Z轴方向上的最小值分别为XR_min＝0、YR_min＝0、ZR_min＝0(米)，最大值分别为XR_max、YR_max、ZR_max(米)；对计算域进行网格剖分，在X轴、Y轴、Z轴方向的剖分步长均为最终剖分步长Δx_CA；最终剖分步长Δx_CA来自步骤三的优化结果；-Z轴为重力方向；

设定网格标识为(j,i,k)，其中j、i和k均为整数；

其中j表示X轴方向上网格的标号，j取值范围为0～(RMTcell-1)：

i表示Y轴方向上网格的标号，i取值范围为0～(RLTcell-1)：

k表示Z轴方向上网格的标号，k取值范围为0～(RNTcell-1)：

每一个网格均为一个正方体，边长为Δx_CA。

按照步骤一的方式对每个网格赋予状态变量，具体包括：

温度tem(j,i,k)，状态state(j,i,k)，固相分数fs(j,i,k)，一个时间步长内固相分数改变量dfs(j,i,k)，生长速度vtip(j,i,k)，界面曲率kcurve(j,i,k)，界面各向异性对应的X轴方向、Y轴方向和Z轴方向三个取向角gangx(j,i,k)，gangy(j,i,k)，gangz(j,i,k)，三个欧拉角ang1(j,i,k)，ang(j,i,k)和ang2(j,i,k)。状态state计算初始时刻为0表明液态，生长态时为1，凝固态为2。

再进一步，步骤四中，在设定的温度梯度和冷却速度下，计算金属液单向凝固过程中直至凝固结束时每个时间步长内每个网格对应的三个欧拉角包括：

步骤四二：针对计算域中每个网格(j,i,k)按照步骤二一的方法计算每个网格(j,i,k)的温度tem(j,i,k)；

步骤四三：柱状晶形核计算：采用基于高斯分布的形核公式计算随着过冷度ΔT变化，金属液中柱状晶的形核数：

其中ΔT为熔体过冷度＝T_l-tem(j,i,k),℃；

式中n为过冷度ΔT所对应的铸件合金三维形核密度,n_max为铸件合金三维最大形核密度，ΔT_N为铸件合金最大形核过冷度，ΔT_σ为铸件合金标准方差过冷度；

计算实际铸件在k＝(RNTcell-1)，i在[0,RLTcell-1]范围内变化，j在[0,RMTcell-1]范围内变化时对应的底部网格数N_bottom，当过冷度为ΔT时，采用底部网格数N_bottom计算获得铸件底部形核个数为：

这些形核粒子随机分布于底部网格中：在N_bottom个底部网格中，随机选定一个网格，若选定网格对应的state为0且温度小于液相线温度T_l，则所述选定网格为形核网格，将形核网格对应的state从0变为1；所述形核网格对应的欧拉角进行随机选取：欧拉角ang1为[0,360°]之间的一个随机值，欧拉角ang为[0,180°]之间的一个随机值，欧拉角ang2为[0,360°]之间一随机值。

再进一步，步骤四中，确定同一晶粒的方法包括：

步骤四四：对于state＝1的形核网格，按照步骤二三和步骤二四的方法进行晶粒生长计算以及捕捉周围液态网格的计算；

步骤四五：随着凝固时间的延长，重复步骤四二至步骤四四，直至计算域中所有网格所对应的state值均等于2，计算停止；输出每个网格对应的三个欧拉角，具有相同欧拉角网格中的晶粒作为同一晶粒；从而实现实际铸件凝固过程的晶粒组织三维数值预测。

本发明方法通过温度计算实现对金属液凝固过程的预测，包括步骤二一的具体过程，降温的过程即为金属液凝固过程。在给定冷却速度下，铸件中的温度从高温向低温进行转变，温度降低至合金液相线温度时，合金开始凝固(步骤二一)。凝固就意味着晶粒形成，晶粒的形成始于形核(步骤二二)，然后生长。本发明借助于元胞自动机技术再现晶粒的生长过程，通过一个元胞捕捉它周围的邻居元胞，并将元胞本身的属性(例如欧拉角)赋予被捕获的元胞。当凝固结束后具有相同欧拉角的元胞就组成一个晶粒。

作为常用飞机发动机涡轮叶片制备方法，单向凝固技术用来制备具有相似取向的晶粒，即柱状晶，从而使生产的叶片在特定方向上具有优异的高温抗蠕变性能。如果工艺参数选择不当引起等轴晶形成，则会导致产品报废。因此本发明方法模拟单向凝固过程中晶粒组织演化，首先评价选定剖分步长下计算结果对网格尺寸的敏感性、晶粒生长速度计算中考虑了界面能各向异性，更加准确地预测了单向凝固过程中晶粒组织形成，解决了现有方法计算量大、计算时间长的问题。

本发明方法适用于预测单向凝固过程中各种工艺条件下铸件内晶粒组织形成。利用本发明方法可以更为准确的预测晶粒组织，从多个方面为优化工艺提供帮助，市场应用潜力巨大，且绿色环保，一旦被广泛采用，将有百万元以上的产值。

下面通过具体实施例对本发明的效果进行验证：

具体实施例：

模型合金为Ni基高温合金(Ni-5wt％Nb)。合金热物性参数以及模拟参数示于表1。单向凝固过程工艺参数示于表2。

表1 Ni-5wt％Nb合金热物性参数和模拟所需参数

表2工艺参数

图2和图3对应的工艺参数一为G＝50℃/m，

对应的参数如表1所示。分别研究网格尺寸为2.5μm和5μm时柱状晶竞争生长。柱状晶的欧拉角分别为(0°,10°,0°)和(0°,30°,0°)。可以看出尽管网格尺寸大小对晶粒生长的最终高度有一定影响，但是对晶界夹角影响较小(晶界夹角分别为8.2°和7.2°)。晶界夹角是评价是否形成取向相似晶粒的重要指标。当选择网格尺寸不同时，模拟所得晶界夹角差值的绝对值小于5°，则计算结果不受网格尺寸影响。

图4至图11所示，对应的工艺参数二为G＝100℃/m，

对应的参数如表1所示。分别研究网格尺寸为2.5μm和5μm时柱状晶竞争生长。图4和图5中柱状晶的欧拉角分别为(0°,10°,0°)和(45°,45°,45°)；图6和图7中柱状晶的欧拉角分别为(10°,10°,10°)和(45°,45°,45°)；图8和图9中柱状晶的欧拉角分别为(10°,10°,10°)和(25°,25°,25°)；图10和图11中柱状晶的欧拉角分别为(10°,10°,10°)和(5°,5°,5°)。可以看出，具有不同欧拉角的柱状晶的竞争生长在图4至图11中均不受网格尺寸的影响，即晶界夹角差值的绝对值均小于5°。所以在后续模拟预测中可以选择网格尺寸为5μm。

图12对应的工艺为G＝100℃/m，

(工艺参数二)，展示了铸件中垂直于Y轴的中心截面。从实验结果中可以看出柱状晶在单向凝固过程中存在竞争生长。在所选择的工艺条件下没有等轴晶组织形成。图13对应的工艺为G＝10℃/m，

(工艺参数三)，为铸件中垂直于Y轴的中心截面。由实验结果可以看出降低温度梯度且增加冷却速度促进等轴晶形成，此时铸件内没有柱状晶。所以工艺参数选择对凝固组织形成有明显影响，若工艺参数选择不当将会引起等轴晶形成。

图14至图16是针对工艺参数二，模拟所得为柱状晶组织。其中图14为晶粒组织三维展示；图15为铸件中垂直于X轴的中心截面，有些柱状晶是由其位置延伸过来的，所以在截面上显示为短长形或近圆形。但这些短长形或近圆形晶粒并非等轴晶，而是柱状晶被截后的断面。图16为铸件中垂直于Y轴的中心截面，该截面可以较为明显的再现柱状晶形貌，可以看到柱状晶存在明显的竞争生长现象，凝固过程中没有等轴晶形成。图16中模拟结果与图12中的实验吻合较好。

图17至图19为针对工艺参数三，模拟所得等轴晶组织。其中图17为晶粒组织三维展示；图18为铸件中垂直于X轴的中心截面。降低温度梯度且增加冷却速度促进熔体内部局部形核及晶粒长大，所以形成等轴晶组织。铸件上方的等轴晶组织细小，下方等轴晶组织粗大。相比于图15，有些靠近空气的网格单元中也形成了等轴晶，这是冷却速度增大的结果。图19为铸件中垂直于Y轴的中心截面，该截面可以较为明显的再现等轴晶组织，且铸件上部和下部的晶粒尺寸存在差别。尽管图19中模拟所得等轴晶相比于图13中的晶粒略微粗大，但就晶粒组织形态而言与实验吻合较好。

虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其它所述实施例中。

Claims

1.一种金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于包括：

2.根据权利要求1所述的金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于，

步骤一中，设定所述计算域尺寸为X米×Y米×Z米；计算域在X轴、Y轴、Z轴方向的最小值分别为X_min＝0、Y_min＝0、Z_min＝0，最大值分别为X_max、Y_max、Z_max；设定-Z轴为重力方向；

设定网格标识为(j,i,k)，其中j、i和k均为整数；

其中j表示X轴方向上网格的标号，j取值范围为0～(MTcell-1)：

i表示Y轴方向上网格的标号，i取值范围为0～(LTcell-1)：

k表示Z轴方向上网格的标号，k取值范围为0～(NTcell-1)：

且MTcell>100，LTcell＝5；

对每个网格基于所述网格标识赋予状态变量包括：

3.根据权利要求2所述的金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于，

步骤二中结合状态变量计算金属液单向凝固过程中直至凝固结束时每个网格对应的三个欧拉角ang1(j,i,k)，ang(j,i,k)，ang2(j,i,k)的方法包括：

式中T_l为液相线温度，G为温度梯度，

为冷却速度，time为凝固时间；

d_z为网格(j,i,k)在Z轴方向距离Z＝0平面的距离：

d_z＝k×Δx_CA；

使j_nu1＝MTcell/4，取整数，i_nu1＝LTcell/2，取整数，k_nu1＝0；

4.根据权利要求3所述的金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于，步骤二中，确定同一晶粒的过程包括：

1)针对每个形核网格采用三个欧拉角构建旋转矩阵A：

其中，

pa1＝cos(ang1)cos(ang2)-sin(ang1)sin(ang2)cos(ang)，

pb1＝sin(ang1)cos(ang2)+cos(ang1)sin(ang2)cos(ang)，

pc1＝sin(ang2)sin(ang)，

pa2＝-cos(ang1)sin(ang2)-sin(ang1)cos(ang2)cos(ang)，

pb2＝-sin(ang1)sin(ang2)+cos(ang1)cos(ang2)cos(ang)，

pc2＝cos(ang2)sin(ang)，

pa3＝sin(ang1)sin(ang)，

pb3＝-cos(ang1)sin(ang)，

pc3＝cos(ang)；

(fsx)'＝pa1·fsx+pa2·fsy+pa3·fsz，

(fsy)'＝pb1·fsx+pb2·fsy+pb3·fsz，

(fsz)'＝pc1·fsx+pc2·fsy+pc3·fsz；

对于gangx：

若[(fsx)']²+[(fsz)']²＝0，则

若(fsx)'＝0且(fsz)'＞0，则

若(fsx)'＝0且(fsz)'＜0，则

若(fsz)'＝0且(fsx)'＞0，则gangx＝0；

若(fsz)'＝0且(fsx)'＜0，则gangx＝π；

若(fsz)'≠0且(fsx)'≠0，则

对于gangy：

若[(fsy)']²+[(fsz)']²＝0，则

若(fsy)'＝0且(fsz)'＞0，则

若(fsy)'＝0且(fsz)'＜0，则

若(fsz)'＝0且(fsy)'＞0，则gangy＝0；

若(fsz)'＝0且(fsy)'＜0，则gangy＝π；

若(fsz)'≠0且(fsx)'≠0，则

对于gangz：

若[(fsx)']²+[(fsy)']²＝0，则

若(fsx)'＝0且(fsy)'＞0，则

若(fsx)'＝0且(fsy)'＜0，则

若(fsy)'＝0且(fsx)'＞0，则gangz＝0；

若(fsy)'＝0且(fsx)'＜0，则gangz＝π；

若(fsy)'≠0且(fsx)'≠0，则

5)再计算考虑了界面各向异性影响的生长速度：

式中ΔT表示过冷度：ΔT＝T_l-tem；

6)固相分数改变量和固相分数的计算：

对计算域内非边缘区域的网格，每个网格存在26个相邻网格，分为第一层近邻，第二层近邻和第三层近邻；对于网格(j,i,k)，将第一层近邻网格记录为包括fs_left(j-1,i,k)，fs_right(j+1,i,k)，fs_up(j,i,k+1)，fs_down(j,i,k-1)，fs_front(j,i+1,k)，fs_back(j,i-1,k)；第二层近邻网格记录为包括fs_upleft(j-1,i,k+1)，fs_upright(j+1,i,k+1)，fs_downleft(j-1,i,k-1)，fs_downright(j+1,i,k-1)，fs_frontup(j,i+1,k+1)，fs_frontdown(j,i+1,k-1)，fs_frontleft(j-1,i+1,k)，fs_frontright(j+1,i+1,k)，fs_backup(j,i-1,k+1)，fs_backdown(j,i-1,k-1)，fs_backleft(j-1,i-1,k)，fs_backright(j+1,i-1,k)；第三层近邻网格记录为包括fs_frontupleft(j-1,i+1,k+1)，fs_frontupright(j+1,i+1,k+1)，fs_{frontdownleft}(j-1,i+1,k-1)，fs_{frontdownright}(j+1,i+1,k-1)，fs_backupleft(j-1,i-1,k+1)，fs_backupright(j+1,i-1,k+1)，fs_backdownleft(j-1,i-1,k-1)，fs_{backdownright}(j+1,i-1,k-1)；对26个相邻网格通过其状态state确定其形貌因子MFG；在每一个时间步长Δt内，MFG的初始值为0，然后进行累加赋值；

对于第一层近邻网格，以fs_left(j-1,i,k)为例，若fs_left(j-1,i,k)的state＝2，则其MFG＝MFG+1；若fs_left(j-1,i,k)的state≠2，则其MFG＝MFG+0；

若fs_upleft(j-1,i,k+1)的state≠2，则MFG＝MFG+0；

若fs_frontupleft(j-1,i+1,k+1)的state≠2，则MFG＝MFG+0；

根据形貌因子MFG计算一个时间步长内相应网格的固相分数改变量和固相分数：

fs(time)＝fs(time-Δt)+dfs，

式中Δt为时间步长，

D_l为溶质在液相中的扩散系数；fs(time)为当前时刻下固相分数，fs(time-Δt)为上一时刻下固相分数；min(MFG/3,1)表示

5.根据权利要求4所述的金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于，步骤二中，确定同一晶粒的过程还包括：

步骤二四：将state＝2的网格(j,i,k)作为母网格，对母网格周围26个邻居网格的state值进行判断，若某一邻居网格的state＝0，则使该邻居网格的state从0变为1，表明被捕获的该邻居网格不再是液相，将进入生长状态，同时对被捕获的该邻居网格赋予母网格对应的三个欧拉角；

步骤二五：逐次选择不同的剖分步长Δx_CA，使Δx_CA依次为step1，step2，step3，……且step1<step2<step3<……，对应执行步骤步骤一至步骤二四，每一次计算直至凝固过程结束，即所有(j,i,NTcell-1)网格的state值均等于2，(j,i,NTcell-1)代表计算域中k取最大值对应的最顶层网格；然后输出每一次计算结束时，每个网格(j,i,k)对应的三个欧拉角，所有网格对应的欧拉角形成剖分步长确认文件grain.txt。

6.根据权利要求5所述的金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于，

步骤三中，将对应Δx_CA依次为step1，step2，step3，……的剖分步长确认文件记为grain1.txt，grain2.txt，grain3.txt，……，对每个剖分步长确认文件分别提取晶界夹角值θ₁，θ₂，θ₃，……，若不同剖分步长Δx_CA对应的晶界夹角差值的绝对值大于或等于5°,表明欧拉角的计算结果受剖分步长取值的影响超过安全阈值，则减小step1，step2，step3，……中的最大值，重新执行步骤一；直至不同剖分步长Δx_CA对应的晶界夹角差值的绝对值小于5°，表明欧拉角的计算结果受剖分步长取值的影响在安全阈值以内，从step1，step2，step3，……中选择最大值，作为最终剖分步长。

7.根据权利要求6所述的金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于，步骤四中，采用网格标识对每个网格赋予状态变量包括：

步骤四一：根据实际铸件设定计算域尺寸为XR米×YR米×ZR米，计算域在X轴、Y轴、Z轴方向上的最小值分别为XR_min＝0、YR_min＝0、ZR_min＝0，最大值分别为XR_max、YR_max、ZR_max；对计算域进行网格剖分，在X轴、Y轴、Z轴方向的剖分步长均为最终剖分步长Δx_CA；

设定网格标识为(j,i,k)，其中j、i和k均为整数；

其中j表示X轴方向上网格的标号，j取值范围为0～(RMTcell-1)：

i表示Y轴方向上网格的标号，i取值范围为0～(RLTcell-1)：

k表示Z轴方向上网格的标号，k取值范围为0～(RNTcell-1)：

按照步骤一的方式对每个网格赋予状态变量。

8.根据权利要求7所述的金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于，步骤四中，在设定的温度梯度和冷却速度下，计算金属液单向凝固过程中直至凝固结束时每个时间步长内每个网格对应的三个欧拉角包括：

计算实际铸件在k＝(RNTcell-1)，i在[0,RLTcell-1]范围内变化，j在[0,RMTcell-1]范围内变化时对应的底部网格数N_bottom，采用底部网格数N_bottom计算获得铸件底部形核个数为：

在N_bottom个底部网格中，随机选定一个网格，若选定网格对应的state为0且温度小于液相线温度T_l，则所述选定网格为形核网格，将形核网格对应的state从0变为1；所述形核网格对应的欧拉角进行随机选取：欧拉角ang1为[0,360°]之间的一个随机值，欧拉角ang为[0,180°]之间的一个随机值，欧拉角ang2为[0,360°]之间一随机值。

9.根据权利要求8所述的金属液单向凝固过程中晶粒组织三维数值预测方法，其特征在于，步骤四中，确定同一晶粒的方法包括：