CN106944607A - 一种孕育合金晶粒组织数值预测方法 - Google Patents

Abstract

本发明涉及一种孕育合金晶粒组织数值预测的方法，属于晶粒组织的仿真预测方法。本发明为了解决现有技术中孕育合金晶粒组织数值预测中三维宏观场的计算效率不高、无法准确预测晶粒组织的缺点，而提出一种孕育合金晶粒组织数值预测的方法，包括：对铸造系统进行宏观尺度网格剖分；对于所有非铸件网格，计算能量守恒方程，获得温度场分布；对于所有铸件的网格，计算能量守恒方程和成分守恒方程；对于铸件网格，计算动量守恒方程；采用元胞自动机法进行晶粒组织模拟，得到当前时刻的铸件内晶粒组织分布；重复上述步骤，直至所有铸件网格所对应的固相分数为1；最终输出铸件内晶粒组织分布。本发明适用于孕育合金晶粒组织的仿真及数值预测。

Description

一种孕育合金晶粒组织数值预测方法

技术领域

本发明涉及一种孕育合金晶粒组织数值预测方法，属于晶粒组织的仿真预测方法。

背景技术

晶粒组织是评价铸造产品性能的重要指标。铸造产品的使用性能根据晶粒组织的差别会呈现显著变化。表征晶粒组织的参量如晶粒大小、晶粒形态、晶粒分布均匀度对铸造产品的机械性能和物理性能有着强烈影响。呈圆柱状的晶粒称为柱状晶，其特点为晶界面积小、位向一致。具有大量柱状晶组织的铸造产品性能具有明显的方向性，沿着柱状晶生长方向的性能好且垂直于柱状晶生长方向的性能差，同时柱状晶生长前沿为气体和第二相杂质富集区，该区域极易产生热裂。呈近圆形状的晶粒称为等轴晶，相比于柱状晶，等轴晶晶粒之间位向随机分布，因此铸造产品的性能更均匀稳定。等轴晶尺寸小、个数多、晶界面积大，晶界面积大促使杂质和缩松缺陷分布更加离散，避免形成杂质和缩松缺陷聚集区。采用钢铁材料和塑性较差的有色金属材料制备铸造产品时希望获得全部细小等轴晶组织，从而提高产品塑性和抗腐蚀性。目前细化晶粒组织的主要手段为孕育处理，即向金属液中添加孕育剂达到细化晶粒的目的。细化剂加入后对合金的细化程度受很多因素影响，如冷却条件、合金浇注时的过热度、细化剂含量等。不同因素之间的交互作用非常复杂，同时为保证铸造产品制备后的结构完整性，很多凝固过程中的变量无法进行实时监测，因此采用实验手段分析、研究和控制孕育合金晶粒组织具有盲目性且浪费大量物力和财力，不利于环保。

孕育合金晶粒组织数值预测结合了基础凝固理论-计算机学-实验三方面研究，多学科交叉研究为孕育合金材料的铸造产品成型提供参数选择和理论指导。通过晶粒组织演化的数值预测可以明晰工艺参数改变如何影响凝固组织形成，获得关键工艺参数组合，有效缩短铸造工艺研发时间，提高研发效率。随着凝固模型的逐渐完善，晶粒组织数值预测将更精确，成为提高铸造产品质量的新途径。

目前孕育合金晶粒组织数值预测方法的局限性：第一仅考虑二维宏观温度场、流场和成分场变化，而实际铸件都为三维，但三维宏观场计算量大且计算时间长；第二仅通过增加形核密度体现孕育处理的作用，因此晶粒细化机理为机械碰撞导致晶粒停止生长，而大量研究表明溶质富集降低已形核晶粒生长速度从而给继续形核提供时间为主要晶粒细化机理，因此考虑形核密度的同时还要计算局部溶质富集对晶粒生长速度的影响。这就要求所开发的孕育合金晶粒组织数值预测方法既可以在三维方向上计算宏观场又可以减少计算量且缩短计算时间，同时晶粒生长速度计算中需要考虑局部溶质富集的作用。

发明内容

本发明的目的是为了解决现有技术中孕育合金晶粒组织数值预测中三维宏观场的计算效率不高、无法准确预测晶粒组织的缺点，而提出一种孕育合金晶粒组织数值预测方法，包括：

步骤一：对铸造系统进行宏观尺度网格剖分，将铸造系统划分为标号为(i,j,k)_chan的若干网格，i、j、k分别表示沿X轴、Y轴、Z轴的坐标分量，其中X轴、Y轴、Z轴为相互正交的任意坐标轴；chan表示网格的类型；

步骤二：对于所有类型不为“铸件”的网格，计算能量守恒方程，获得温度场分布；

步骤三：对于所有类型为“铸件”的网格，计算能量守恒方程和成分守恒方程，得到温度场以及成分场；

步骤四：对于所有类型为“铸件”的网格，计算动量守恒方程，得到网格中的金属液流动速度；

步骤五：采用元胞自动机法进行晶粒组织模拟，得到当前时刻的铸件内晶粒组织分布；

步骤六：重复步骤二至步骤五，直至所有类型为“铸件”的网格所对应的固相分数为1；最终输出铸件内晶粒组织分布。

本发明的有益效果为：本发明模拟铸造系统三维方向的传热、传质和对流传输与实际更为接近且采用不同尺寸的网格对铸造系统进行剖分减少计算时间、提高计算效率。采用自由生长形核模型和基于高斯分布的形核模型描述孕育合金凝固过程中的形核现象、随机选择晶粒以球形方式或枝形方式完成生长、计算枝晶生长速度时考虑了局部溶质富集对晶粒生长速度的影响，解决了目前晶粒组织数值预测更多基于二维计算、三维计算量大、对物理现象捕捉不全面的问题，为铸造产品制备过程中铸造工艺改进以及产品力学性能预测提供了理论指导和数据参考。本发明适用于各类尺寸的砂型和金属型铸造过程中孕育合金晶粒组织预测。利用本发明可以更为准确的预测晶粒组织形貌分布，为铸造工艺设计和优化提供帮助，市场应用潜力巨大，一旦被广泛采用，将有千万元以上的产值。

附图说明

图1为本发明的孕育合金晶粒组织数值预测方法的流程图；

图2为本发明选用的铸造型腔的一个实施例的实物图；

图3为熔体过热度为50℃时实验所得Al-15wt％Cu二元合金铸件中间截面；

图4为对图3中的“PP”区域电子背散射衍射测试所得等轴晶组织形貌；

图5为一个实施例中对铸件温度场、流场和成分场进行三维计算所得铸件中间截面处凝固晶粒组织效果图；

图6为一个实施例中对铸件温度场、流场和成分场进行二维计算所得铸件中间截面处凝固晶粒组织；

图7为对铸件温度场、流场和成分场进行三维计算所得S1截面50s时流场分布；

图8为对铸件温度场、流场和成分场进行二维计算所得S1截面50s时流场分布。

具体实施方式

具体实施方式一：本实施方式的孕育合金晶粒组织数值预测方法，包括如下步骤：

步骤一：对铸造系统进行宏观尺度网格剖分，将铸造系统划分为标号为(i,j,k)_chan的若干网格，i、j、k分别表示沿X轴、Y轴、Z轴的坐标分量，其中X轴、Y轴、Z轴为相互正交的任意坐标轴；chan表示网格的类型。X轴、Y轴、Z轴的选取可以依照实际情况而定。图2示出了铸造型腔的实物图，并且是俯视图，选取三轴的方法可以为：沿图像方向垂直向下为X轴正方向、沿图像方向水平向右为Y轴正方向、沿图像纸面向外为Z轴正方向。为了便于计算，选择坐标轴时可以使得任意两个坐标轴形成的平面与铸造系统平面保持平行。

步骤二：对于所有类型不为“铸件”的网格，计算能量守恒方程，获得温度场分布。

温度场分布可以按照如下公式计算：

[H]＝c_pT

其中c_p为比热(J/kg K)，ρ_m为密度(kg/m³)，λ_m为导热系数(W/m K)，t为时间(s)，T为温度(℃),[H]为热焓(J/kg)。

步骤三：对于所有类型为“铸件”的网格，计算能量守恒方程和成分守恒方程，得到

温度场以及成分场。

能量守恒方程可以按照如下公式计算：

h_s＝c_pT

h_l＝c_pT+ΔH

其中h_s和h_l分别为固相和液相热焓(J/kg)，λ_l为液体导热系数(W/m·K)，ΔH为结晶潜热(J/kg)，T_s固相线温度(℃)，T_l液相线温度(℃)，U_l液体流动速度(m/s),f_l为液相分数。

成分守恒方程可以按照如下公式计算：

C_mix＝f_sC_s+f_lC_l

C_l＝(C_o+C_l ^*)/2

C_s＝k_p·C_l

其中C_mix为混合成分(wt％)，C_l为平均液相成分(wt％)，C_l ^*为固液界面处液相成分(wt％)，C_o为合金初始成分(wt％)，C_s为固相成分(wt％)，k_p为平衡分配系数(无量纲)，m_l为液相线斜率(℃/wt％),D_l为液相溶质扩散系数(m²/s),f_s＝1-f_l为固相分数(无量纲)。

步骤四：对于所有类型为“铸件”的网格，计算动量守恒方程，得到网格中的金属液流动速度。可以按照如下公式计算：

λ_c＝7.63(D_lΓ)^1/3(Cool_rate)^-1/3C_o ^-1/4

其中U_l为液体流动速度且0s时的值为0m/s，ρ_l为液相密度(kg/m³)，P为压强(Pa)，μ_l为液体粘度(Pa·s)，β_T为温度膨胀系数(1/℃)，β_C为成分膨胀系数(1/wt％)，T_ref为参考温度(℃)等于液相线温度T_l，C_ref为参考成分(wt％)等于合金初始成分C_o，为重力加速度m/s²，K_per为渗透率(m²)，Γ为吉布斯汤姆森系数(℃m)，λ_c为枝晶臂间距(m)，Cool_rate为平均冷却速率(℃/s),Δt为时间步长(s)，T^t为t时刻下温度，T^t-Δt为(t-Δt)时刻下温度。

步骤五：采用元胞自动机法(即CA法)进行晶粒组织模拟，得到当前时刻的铸件内晶粒组织分布。

步骤六：重复步骤二至步骤五，直至所有类型为“铸件”的网格所对应的固相分数为1；最终输出铸件内晶粒组织分布。

具体实施方式二：本实施方式与具体实施方式一不同的是：

铸造系统在X轴、Y轴、Z轴方向上的最小值分别为X_min、Y_min、Z_min，在X轴、Y轴、Z轴方向上的最大值分别为X_max、Y_max、Z_max，步骤一具体为：

步骤一一：选取需要模拟的晶粒组织的截面S1、第一辅助平面S2、第二辅助平面S3以及坐标轴；选取的方法为：

若选取垂直于Z轴且所处位置为Z_s1的平面，则截面S1在X轴方向和Y轴方向采用的网格剖分步长分别为Δx₁和Δy₁；以S1为基准面，平行于S1截面且沿着Z轴正方向、距S1截面的距离为δ米的截面为S2，平行于S1截面且沿着Z轴负方向、距S1截面为δ米的界面为S3；垂直于Z轴且分别平行于S2和S3截面且位于S2和S3截面之间，每隔Δz₁米选取一个截面，所述选取的截面以及截面S2和截面S3的网格剖分步长均为Δx₁米和Δy₁米；垂直于Z轴、平行于S1截面且与S1截面之间的距离大于δ米的其他截面的剖分步长为Δx₂米和Δy₂米，沿着Z轴的剖分步长为Δz₂米。

若选取垂直于X轴且所处位置为X_s1的平面，则截面S1在Y轴方向和Z轴方向采用的网格剖分步长分别为Δy₁和Δz₁；以S1为基准面，平行于S1截面且沿着X轴正方向、距S1截面的距离为δ米的截面为S2，平行于S1截面且沿着X轴负方向、距S1截面为δ米的界面为S3；垂直于X轴且分别平行于S2和S3截面且位于S2和S3截面之间，每隔Δx₁米选取一个截面，所述选取的截面以及截面S2和截面S3的网格剖分步长均为Δy₁米和Δz₁米；垂直于X轴、平行于S1截面且与S1截面之间的距离大于δ米的其他截面的剖分步长为Δy₂米和Δz₂米，沿着X轴的剖分步长为Δz₂米。

若选取垂直于Y轴且所处位置为Y_s1的平面，则截面S1在Y轴方向和Z轴方向采用的网格剖分步长分别为Δz₁和Δx₁；以S1为基准面，平行于S1截面且沿着Y轴正方向、距S1截面的距离为δ米的截面为S2，平行于S1截面且沿着Y轴负方向、距S1截面为δ米的界面为S3；垂直于Y轴且分别平行于S2和S3截面且位于S2和S3截面之间，每隔Δy₁米选取一个截面，所述选取的截面以及截面S2和截面S3的网格剖分步长均为Δz₁米和Δx₁米；垂直于Y轴、平行于S1截面且与S1截面之间的距离大于δ米的其他截面的剖分步长为Δz₂米和Δx₂米，沿着X轴的剖分步长为Δy₂米。

步骤一二：按照步骤一一划分完成的每个网格的标号记为(i,j,k)_chan，其中chan＝0表示铸件网格，i、j、k均为整数；i的取值范围是1～M,j的取值范围是1～N,k的取值范围是1～L。

若截面S1垂直于Z轴则(0.006m都用δ代替，0.012m为2δ)

若截面S1垂直于X轴，则

若截面S1垂直于Y轴，则

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：

步骤一二中，δ＝0.006m。

若选取垂直于Z轴且所处位置为Z_s1的平面，则Δx₁＝Δy₁，Δz₁＝0.002m，Δx₂＝3Δx₁，Δy₂＝3Δy₁,Δz₂＝3Δz₁。

若选取垂直于X轴且所处位置为X_s1的平面，则Δy₁＝Δz₁，Δx₁＝0.002m，Δz₂＝3Δz₁，Δy₂＝3Δy₁,Δx₂＝3Δx₁。

若选取垂直于Y轴且所处位置为Y_s1的平面，则Δz₁＝Δx₁，Δy₁＝0.002m，Δz₂＝3Δz₁，Δx₂＝3Δx₁,，Δy₂＝3Δy₁。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：

步骤一二中，chan＝1表示砂型网格，chan＝2表示金属铸型网格，chan＝4表示内冷铁网格，chan＝5表示外冷铁网格，chan＝6表示冒口套网格，chan＝7表示保温材料网格，chan＝8表示绝热材料网格。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤五包括：

步骤五一：将截面S1作为晶粒组织模拟计算域，在t时刻下，判断S1截面上的网格(i,j,k)_chan的温度T＝(i,j,k,t)与局部熔体所对应的温度T_local＝(i,j,k,t)＝[T₁-m₁C_o+m₁C_l(i,j,k,t)]的大小关系。若T＝(i,j,k,t)≥T_local＝(i,j,k,t)，则该网格无形核现象发生。

若T＝(i,j,k,t)＜T_local＝(i,j,k,t)且t时刻下的过冷度大于(t-Δt)时刻下的过冷度，则该网格发生形核现象；并执行步骤五二。

步骤五二：针对网格(i,j,k)_chan＝0计算t时刻下的形核密度n＝(i,j,k,t)以及形核核心个数N_num＝(i,j,k,t)；设置第一数值NR，所述第一数值NR随机为1或2。

若NR＝1，则形核核心来源于熔体中的杂质，且

ΔT_nucl(i,j,k,t)＝T_local(i,j,k,t)-T(i,j,k,t)

ΔT_nucl(i,j,k,t-Δt)＝T_local(i,j,k,t-Δt)-T(i,j,k,t-Δt)

若NR＝2，则形核核心来源于孕育剂粒子，在孕育剂粒子的直径范围内随机选取数值dd(i,j,k,t)，计算网格(i,j,k)_chan＝0对应的形核过冷度其中σ为固液界面表面能，ΔS_v为熵；

其中，形核核心个数N_num＝(i,j,k,t)的确定方法为：

若S1截面为垂直于Z轴的平面，则形核核心个数为：

若S1截面为垂直于X轴的平面，则形核核心个数为：

若S1截面为垂直于Y轴的平面，则形核核心个数为：

步骤五三：在t时刻下，S1截面的网格(i,j,k)_chan＝0所对应的形核核心个数N_num＝(i,j,k,t)，每个形核核心的编号记为ni，其中ni介于1至N_num之间；设置第二数值GR，第二数值GR随机为1或2；形核核心个数N_num若大于0且A_ssn(i,j,k,t)小于1，则根据第二数值GR的具体值判断晶粒的生长方式并进行计算；

若GR＝1则形核核心以球形方式生长，球形方式的生长速度计算公式为：

Vs(i,j,k,t)_ni＝para_sphere·D_l/[2R_g(i,j,k,t)_ni]

R_g(i,k,j,t)_ni＝R_g(i,j,k,t-Δt)_ni+Δt·Vs(i,j,k,t)_ni；

若GR＝2则形核核心以枝晶方式生长，枝晶方式的生长速度计算公式为：

Vs(i,j,k,t)_ni＝K₁·(ΔT_c)²+K₂·(ΔT_c)³

K₁＝1.16×10^-4×C_o ^-1.24319

K₂＝5.4×10^-4×C_o ^-2.13518

R_g(i,j,k,t)_ni＝R_g(i,j,k,t-Δt)_ni+Δt·Vs(i,j,k,t)_ni

其中，固液界面前沿熔体中的平均液相成分为：

C_lnear(i,j,k,t)＝C_l(i,j,k,t)；

由于晶粒形核或长大所引起的液相体积减少为：

当A_ssn(i,j,k,t)＞1时，取A_ssn(i,j,k,t)＝1。

<实施例1>

本试验选择Al-15wt％Cu二元合金，采用砂型铸造，单侧放置石墨冷铁。Al-15wt％Cu二元合金的热物性参数和相图数据列于表1。砂型和石墨冷铁的热物性参数列于表2。

表1

表2

图2为本实施例选择的一个铸造型腔的实物图，其中铸件型腔三维几何尺寸为：76mm(X轴)×254mm(X轴)×76mm(Z轴)。铸件型腔右侧放置石墨冷铁，其三维几何尺寸为：76mm(X轴)×64mm(Y轴)×76mm(Z轴)。铸件型腔左侧为直浇道和横浇道。重力方向(Z轴)垂直于纸面。图3为熔体过热度为50℃时实验所得Al-15wt％Cu二元合金铸件中间截面。

对图3中的“PP”区域电子背散射衍射测试所得等轴晶组织形貌，得到如图4所示的等轴晶组织形貌图。

对铸件温度场、流场和成分场进行三维计算，得到如图5所示的铸件中间截面处凝固晶粒的组织效果图。晶粒组织模拟截面为S1，其垂直于X轴且所处位置X_s1＝0.038m。沿着X轴截面S2和S3的位置分别为0.044m和0.032。网格剖分尺寸:Δy₁＝Δz₁＝0.001m，Δx₁＝0.002m，Δz₂＝0.003m，Δy₂＝0.003m。针对S1截面，模拟所得凝固组织均为细小等轴晶，与图3实验所得凝固组织对比较好吻合。在S1截面上选取“PPS”区域，计算该区域平均晶粒尺寸为0.49mm，与图4中实验测得平均晶粒尺寸0.55mm进行对比较好吻合。

对铸件温度场、流场和成分场进行二维计算，得到如图6所示的铸件中间截面处凝固晶粒组织效果图。图6中，晶粒组织模拟截面为S1，其垂直于X轴且所处位置X_s1＝0.038m。网格剖分尺寸:Δy₁＝Δz₁＝0.001m。针对S1截面，模拟所得凝固组织为粗大等轴晶和细小等轴晶混合型组织，与图3实验所得凝固组织对比存在较大差别。在S1截面上选取“PPSS”区域，计算该区域平均晶粒尺寸为1.43mm，与图4中实验测得平均晶粒尺寸0.55mm进行对比差别较大。说明二维模拟无法更准确再现实际凝固过程。

图7和图8为铸件温度场、流场和成分场进行三维计算和二维计算所得S1截面凝固过程中流场分布特征的对比图。根据连续性方程二维模拟所得流动强度大于三维模拟所得流动强度符合流体力学理论。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (8)

1.一种孕育合金晶粒组织数值预测方法，其特征在于，包括如下步骤：

步骤一：对铸造系统进行宏观尺度网格剖分，将铸造系统划分为标号为(i,j,k)_chan的若干网格，i、j、k分别表示沿X轴、Y轴、Z轴的坐标分量，其中X轴、Y轴、Z轴为相互正交的任意坐标轴；chan表示网格的类型；

步骤二：对于所有类型不为“铸件”的网格，计算能量守恒方程，获得温度场分布；

步骤三：对于所有类型为“铸件”的网格，计算能量守恒方程和成分守恒方程，得到温度场以及成分场；

步骤四：对于所有类型为“铸件”的网格，计算动量守恒方程，得到网格中的金属液流动速度；

步骤五：采用元胞自动机法进行晶粒组织模拟，得到当前时刻的铸件内晶粒组织分布；

步骤六：重复步骤二至步骤五，直至所有类型为“铸件”的网格所对应的固相分数为1；最终输出铸件内晶粒组织分布。

2.根据权利要求1所述的孕育合金晶粒组织数值预测方法，其中，所述铸造系统在X轴、Y轴、Z轴方向上的最小值分别为X_min、Y_min、Z_min，在X轴、Y轴、Z轴方向上的最大值分别为X_max、Y_max、Z_max，其特征在于，所述步骤一具体为：

步骤一一：选取需要模拟的晶粒组织的截面S1、第一辅助平面S2、第二辅助平面S3以及坐标轴；选取的方法为：

若选取垂直于Z轴且所处位置为Z_s1的平面，则截面S1在X轴方向和Y轴方向采用的网格剖分步长分别为Δx₁和Δy₁；以S1为基准面，平行于S1截面且沿着Z轴正方向、距S1截面的距离为δ米的截面为S2，平行于S1截面且沿着Z轴负方向、距S1截面为δ米的界面为S3；垂直于Z轴且分别平行于S2和S3截面且位于S2和S3截面之间，每隔Δz₁米选取一个截面，所述选取的截面以及截面S2和截面S3的网格剖分步长均为Δx₁米和Δy₁米；垂直于Z轴、平行于S1截面且与S1截面之间的距离大于δ米的其他截面的剖分步长为Δx₂米和Δy₂米，沿着Z轴的剖分步长为Δz₂米；

若选取垂直于X轴且所处位置为X_s1的平面，则截面S1在Y轴方向和Z轴方向采用的网格剖分步长分别为Δy₁和Δz₁；以S1为基准面，平行于S1截面且沿着X轴正方向、距S1截面的距离为δ米的截面为S2，平行于S1截面且沿着X轴负方向、距S1截面为δ米的界面为S3；垂直于X轴且分别平行于S2和S3截面且位于S2和S3截面之间，每隔Δx₁米选取一个截面，所述选取的截面以及截面S2和截面S3的网格剖分步长均为Δy₁米和Δz₁米；垂直于X轴、平行于S1截面且与S1截面之间的距离大于δ米的其他截面的剖分步长为Δy₂米和Δz₂米，沿着X轴的剖分步长为Δz₂米；

若选取垂直于Y轴且所处位置为Y_s1的平面，则截面S1在Y轴方向和Z轴方向采用的网格剖分步长分别为Δz₁和Δx₁；以S1为基准面，平行于S1截面且沿着Y轴正方向、距S1截面的距离为δ米的截面为S2，平行于S1截面且沿着Y轴负方向、距S1截面为δ米的界面为S3；垂直于Y轴且分别平行于S2和S3截面且位于S2和S3截面之间，每隔Δy₁米选取一个截面，所述选取的截面以及截面S2和截面S3的网格剖分步长均为Δz₁米和Δx₁米；垂直于Y轴、平行于S1截面且与S1截面之间的距离大于δ米的其他截面的剖分步长为Δz₂米和Δx₂米，沿着X轴的剖分步长为Δy₂米；

步骤一二：按照步骤一一划分完成的每个网格的标号记为(i,j,k)_chan，其中chan＝0表示铸件网格，i、j、k均为整数；i的取值范围是1～M,j的取值范围是1～N,k的取值范围是1～L；

若截面S1垂直于Z轴则

$M=\frac{X_{max}-X_{min}}{{\mathrm{\&Delta;x}}_1}N=\frac{Y_{max}-Y_{min}}{{\mathrm{\&Delta;y}}_1},L=\frac{2\mathrm{\&delta;}}{{\mathrm{\&Delta;z}}_1}+\frac{Z_{max}-Z_{min}-2\mathrm{\&delta;}}{{\mathrm{\&Delta;z}}_2};$

若截面S1垂直于X轴，则

$N=\frac{Y_{max}-Y_{min}}{{\mathrm{\&Delta;y}}_1},L=\frac{Z_{max}-Z_{min}}{{\mathrm{\&Delta;z}}_1},M=\frac{2\mathrm{\&delta;}}{{\mathrm{\&Delta;x}}_1}+\frac{X_{max}-X_{min}-2\mathrm{\&delta;}}{{\mathrm{\&Delta;x}}_2};$

若截面S1垂直于Y轴，则

$L=\frac{Z_{max}-Z_{min}}{{\mathrm{\&Delta;z}}_1},M=\frac{X_{max}-X_{min}}{{\mathrm{\&Delta;x}}_1},N=\frac{2\mathrm{\&delta;}}{{\mathrm{\&Delta;y}}_1}+\frac{Y_{max}-Y_{min}-2\mathrm{\&delta;}}{{\mathrm{\&Delta;y}}_2}.$

3.根据权利要求2所述的孕育合金晶粒组织数值预测方法，其特征在于，步骤一二中，δ＝0.006；

若选取垂直于Z轴且所处位置为Z_s1的平面，则Δx₁＝Δy₁，Δz₁＝0.002m，Δx₂＝3Δx₁，Δy₂＝3Δy₁；Δz₂＝3Δz₁；

若选取垂直于X轴且所处位置为X_s1的平面，则Δy₁＝Δz₁，Δx₁＝0.002m，Δz₂＝3Δz₁，Δy₂＝3Δy₁；Δx₂＝3Δx₁；

若选取垂直于Y轴且所处位置为Y_s1的平面，则Δz₁＝Δx₁，Δy₁＝0.002m，Δz₂＝3Δz₁，Δx₂＝3Δx₁；Δy₂＝3Δy₁。

4.根据权利要求3所述的孕育合金晶粒组织数值预测方法，其特征在于，步骤一二中，chan＝1表示砂型网格，chan＝2表示金属铸型网格，chan＝4表示内冷铁网格，chan＝5表示外冷铁网格，chan＝6表示冒口套网格，chan＝7表示保温材料网格，chan＝8表示绝热材料网格。

5.根据权利要求4所述的孕育合金晶粒组织数值预测方法，其特征在于，步骤二具体为：

对于所有类型不为“铸件”的网格，根据如下公式计算能量守恒方程，获得温度场分布：

[H]＝c_pT

$\frac{{\mathrm{\&rho;}}_m\&part;\&lsqb;H\&rsqb;}{\&part;t}={\mathrm{\&lambda;}}_m\&dtri;\&CenterDot;(\&dtri;T)$

其中c_p为比热，ρ_m为密度，λ_m为导热系数，t为时间，T为温度,[H]为热焓。

6.根据权利要求5所述的孕育合金晶粒组织数值预测方法，其特征在于，步骤三具体为：

对于所有类型为“铸件”的网格，根据如下公式计算能量守恒方程和成分守恒方程，得到温度场以及成分场；

能量守恒方程为：

h_s＝c_pT

h_l＝c_pT+ΔH

$\&lsqb;H\&rsqb;=(c_p+\frac{\mathrm{\&Delta;}H}{T_l-T_s})T-\frac{T_s\&CenterDot;\mathrm{\&Delta;}H}{T_l-T_s}$

$\frac{\&part;\&lsqb;H\&rsqb;}{\&part;t}+\&dtri;\&CenterDot;\left(f_l{U}_l{h}_l\right)={\mathrm{\&lambda;}}_l\&dtri;\&CenterDot;(\&dtri;T)$

其中h_s和h_l分别为固相和液相热焓，λ_l为液体导热系数，ΔH为结晶潜热，T_s固相线温度，T_l液相线温度，U_l液体流动速度,f_l为液相分数；

成分守恒方程为：

$\frac{\&part;\left(C_{mix}\right)}{\&part;t}+U_l\&CenterDot;\&dtri;\left(C_l\right)=\&dtri;\&CenterDot;(D_l\&dtri;(C_l\left)\right)$

$f_s=\frac{T_l-T}{T_l-T_s}$

C_mix＝f_sC_s+f_lC_l

${C_l}^{*}=\frac{T_l-m_l{C}_o-T}{-m_l}$

C_l＝(C_o+C_l ^*)/2

C_s＝k_p·C_l

其中C_mix为混合成分，C_l为平均液相成分，C_l ^*为固液界面处液相成分，C_o为合金初始成分，C_s为固相成分，k_p为平衡分配系数，m_l为液相线斜率,D_l为液相溶质扩散系数,f_s＝1-f_l为固相分数。

7.根据权利要求6所述的孕育合金晶粒组织数值预测方法，其特征在于，步骤四具体为：

对于所有类型为“铸件”的网格，根据如下公式计算动量守恒方程，得到网格中的金属液流动速度：

$\begin{array}{c}\frac{\&part;\left(f_l{U}_l\right)}{\&part;t}+\&dtri;\&CenterDot;\left(f_l{U}_l{U}_l\right)=\\ -\frac{1}{{\mathrm{\&rho;}}_l}\&dtri;P+\frac{1}{{\mathrm{\&rho;}}_l}\&dtri;\&CenterDot;({\mathrm{\&mu;}}_l\&dtri;(f_l{U}_l\left)\right)+f_l\stackrel{\&RightArrow;}{g}\left({\mathrm{\&beta;}}_T\right(T-T_{ref})+{\mathrm{\&beta;}}_C(C_l-C_{ref}\left)\right)-\frac{f_l^2{\mathrm{\&mu;}}_l}{{\mathrm{\&rho;}}_l{K}_{per}}{U}_l\end{array}$

$K_{per}=\frac{f_l^3}{{(1-f_l)}^2}\frac{{\mathrm{\&lambda;}}_c^2}{180}$

λ_c＝7.63(D_lΓ)^1/3(Cool_rate)^-1/3C_o ^-1/4

${\mathrm{Cool}}_{rate}=\frac{T^t-T^{t-\mathrm{\&Delta;}t}}{\mathrm{\&Delta;}t}$

其中U_l为液体流动速度且0s时的值为0m/s；ρ_l为金属液密度；P为压强；μ_l为液体粘度；β_T为温度膨胀系数；β_C为成分膨胀系数，T_ref为参考温度等于液相线温度T_l；C_ref为参考成分，等于合金初始成分C_o；为重力加速度m/s²；K_per为渗透率；Γ为吉布斯汤姆森系数；λ_c为枝晶臂间距；Cool_rate为平均冷却速率；Δt为时间步长；T^t为t时刻下温度；T^t-Δt为t-Δt时刻下温度。

8.根据权利要求7所述的孕育合金晶粒组织数值预测方法，其特征在于，步骤五包括：

步骤五一：将截面S1作为晶粒组织模拟计算域，在t时刻下，判断S1截面上的网格(i,j,k)_chan＝0的温度T＝(i,j,k,t)与局部熔体所对应的温度T_local＝(i,j,k,t)＝[T₁-m₁C_o+m₁C₁(i,j,k,t)]的大小关系；

若T＝(i,j,k,t)≥T_local＝(i,j,k,t)，则该网格无形核现象发生；

若T＝(i,j,k,t)＜T_local＝(i,j,k,t)且t时刻下的过冷度大于(t-Δt)时刻下的过冷度，则该网格发生形核现象；并执行步骤五二；

步骤五二：针对网格(i,j,k)_chan＝0计算t时刻下的形核密度n＝(i,j,k,t)以及形核核心个数N_num＝(i,j,k,t)；设置第一数值NR，所述第一数值NR随机为1或2；

若NR＝1，则形核核心来源于熔体中的杂质，且

$\begin{array}{c}n(i,j,k,t)=n(i,j,k,t-\mathrm{\&Delta;}t)+\&lsqb;A_{ssn}(i,j,k,t)\&rsqb;\&CenterDot;\frac{N_{\max 1}}{{\mathrm{\&Delta;T}}_{\mathrm{\&sigma;}}\sqrt{2\mathrm{\&pi;}}}\\ \&CenterDot;\exp \&lsqb;-\frac{1}{2}{\left(\frac{{\mathrm{\&Delta;T}}_{nucl}(i,j,k,t)-{\mathrm{\&Delta;T}}_m}{{\mathrm{\&Delta;T}}_{\mathrm{\&sigma;}}}\right)}^2\&rsqb;\&CenterDot;\&lsqb;{\mathrm{\&Delta;T}}_{nucl}(i,j,k,t)-{\mathrm{\&Delta;T}}_{nucl}(i,j,k,t-\mathrm{\&Delta;}t)\&rsqb;\end{array}$

ΔT_nucl(i,j,k,t)＝T_local(i,j,k,t)-T(i,j,k,t)

ΔT_nucl(i,j,k,t-Δt)＝T_local(i,j,k,t-Δt)-T(i,j,k,t-Δt)；

若NR＝2，则形核核心来源于孕育剂粒子，在孕育剂粒子的直径范围内随机选取数值dd＝(i,j,k,t)，计算网格(i,j,k)_chan＝0对应的形核过冷度其中σ为固液界面表面能，ΔS_v为熵；

其中，形核核心个数N_num＝(i,j,k,t)的确定方法为：

若S1截面为垂直于Z轴的平面，则形核核心个数为：

$N_{num}(i,j,k,t)={\&lsqb;\sqrt{\frac{6}{\mathrm{\&pi;}}}\&CenterDot;n(i,j,k,t)\&rsqb;}^{2/3}\&CenterDot;{\mathrm{\&Delta;x}}_1\&CenterDot;{\mathrm{\&Delta;y}}_1.$

若S1截面为垂直于X轴的平面，则形核核心个数为：

$N_{num}(i,j,k,t)={\&lsqb;\sqrt{\frac{6}{\mathrm{\&pi;}}}\&CenterDot;n(i,j,k,t)\&rsqb;}^{2/3}\&CenterDot;{\mathrm{\&Delta;y}}_1\&CenterDot;{\mathrm{\&Delta;z}}_1.$

若S1截面为垂直于Y轴的平面，则形核核心个数为：

$N_{num}(i,j,k,t)={\&lsqb;\sqrt{\frac{6}{\mathrm{\&pi;}}}\&CenterDot;n(i,j,k,t)\&rsqb;}^{2/3}\&CenterDot;{\mathrm{\&Delta;z}}_1\&CenterDot;{\mathrm{\&Delta;x}}_1.$

步骤五三：在t时刻下，S1截面的网格(i,j,k)_chan＝0所对应的形核核心个数N_num＝(i,j,k,t)，每个形核核心的编号记为ni，其中ni介于1至N_num之间；设置第二数值GR，第二数值GR随机为1或2；形核核心个数N_num若大于0且A_ssn(i,j,k,t)小于1，则根据第二数值GR的具体值判断晶粒的生长方式并进行计算；

若GR＝1则形核核心以球形方式生长，球形方式的生长速度计算公式为：

Vs(i,j,k,t)_ni＝para_sphere·D_l/[2R_g(i,j,k,t)_ni]

${\mathrm{para}}_{sphere}=\left(\frac{-{\mathrm{ppara}}_{sphere}}{2{\mathrm{\&pi;}}^{1/2}}\right)+{(\frac{{{\mathrm{ppara}}_{sphere}}^2}{4\mathrm{\&pi;}}-{\mathrm{ppara}}_{sphere})}^{1/2}$

${\mathrm{ppara}}_{sphere}=\frac{2\&lsqb;{C_l}^{*}(i,j,k,t)-C_{lnear}(i,j,k,t)\&rsqb;}{(1-k_p)\&CenterDot;{C_l}^{*}(i,j,k,t)}$

R_g(i,k,j,t)_ni＝R_g(i,j,k,t-Δt)_ni+Δt·Vs(i,j,k,t)_ni；

若GR＝2则形核核心以枝晶方式生长，枝晶方式的生长速度计算公式为：

Vs(i,j,k,t)_ni＝K₁·(ΔT_c)²+K₂·(ΔT_c)³

K₁＝1.16×10^-4×C_o ^-1.24319

K₂＝5.4×10^-4×C_o ^-2.13518

${\mathrm{\&Delta;T}}_c=m_l{C}_o\&lsqb;1-\frac{1}{1-{\mathrm{\&Omega;}}_c(1-k_p)}\&rsqb;$

${\mathrm{\&Omega;}}_c=\frac{\&lsqb;{C_l}^{*}(i,j,k,t)-C_{lnear}(i,j,k,t)\&rsqb;}{(1-k_p)\&CenterDot;{C_l}^{*}(i,j,k,t)}$

R_g(i,j,k,t)_ni＝R_g(i,j,k,t-Δt)_ni+Δt·Vs(i,j,k,t)_ni

其中，固液界面前沿熔体中的平均液相成分为：

C_lnear(i,j,k,t)＝C_l(i,j,k,t)；

由于晶粒形核或长大所引起的液相体积减少为：

$A_{ssn}(i,j,k,t)=\left\{{\mathrm{\&Sigma;}}_{ni=1}^{N_{num}(i,j,k,t)}\frac{4\mathrm{\&pi;}}{3}{\&lsqb;R_g{(i,j,k,t)}_{ni}+\frac{2D_l}{Vs{(i,j,k,t)}_{ni}}\&rsqb;}^3\right\}/({\mathrm{\&Delta;x}}_1\&CenterDot;{\mathrm{\&Delta;y}}_1\&CenterDot;{\mathrm{\&Delta;z}}_1);$

当A_ssn(i,j,k,t)＞1时，取A_ssn(i,j,k,t)＝1。