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

Abstract

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

Description

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

技术领域

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

背景技术

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

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

目前合金晶粒组织数值预测方法的局限性：第一仅考虑二维宏观温度场、流场和成分场变化，而实际铸件都为三维，但三维宏观场计算量大且计算时间长；第二仅考虑流动对温度场和成分场分布的影响，通过温度场和成分场变化间接影响晶体生长，没有直接考虑流动对晶体生长速度的影响。这就要求所开发的合金晶粒组织数值预测方法既可以在三维方向上计算宏观场又可以减少计算量且缩短计算时间，同时晶粒生长速度计算中需要考虑金属液流动的作用。

发明内容

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

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

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

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

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

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

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

本发明的有益效果为：本发明的合金晶粒组织数值预测方法，模拟铸造系统三维方向的传热、传质和对流传输与实际更为接近且采用不同尺寸的网格对铸造系统进行剖分减少计算时间、提高计算效率。采用基于高斯分布的形核模型描述合金凝固过程中的形核现象、计算枝晶生长速度时考虑了液体流动速度的影响，解决了目前晶粒组织数值预测更多基于二维计算、三维计算量大、对物理现象捕捉不全面的问题，为铸造产品制备过程中铸造工艺改进以及产品力学性能预测提供了理论指导和数据参考。

本发明适用于各类尺寸的砂型和金属型铸造过程中合金晶粒组织预测。利用本发明可以更为准确的预测晶粒组织形貌分布，为铸造工艺设计和优化提供帮助，市场应用潜力巨大，一旦被广泛采用，将有千万元以上的产值。

附图说明

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

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

图3为实验所得Sn-6wt％Pb合金铸件中截面上晶粒组织；

图4为模拟铸件三维凝固所得Sn-6wt％Pb合金S1截面上晶粒组织；

图5模拟铸件二维凝固所得Sn-6wt％Pb合金S1截面上晶粒组织；

图6(a)对铸件温度场、流场和成分场进行三维计算所得S1截面90s时流场分布；

图6(b)对铸件温度场、流场和成分场进行二维计算所得S1截面90s时流场分布。

具体实施方式

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

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

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

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

[H]＝c_mpT

其中c_mp为比热(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轴则

若截面S1垂直于X轴，则

若截面S1垂直于Y轴，则

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

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

步骤一二中，δ＝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₁。

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

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

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

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

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

步骤五包括：

步骤五一：将截面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)；

计算形核核心的计算方式采用基于高斯分布的形核模型，具体计算公式为：

Δ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)

其中，形核核心个数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之间；形核核心个数N_num若大于0且A_ssn(i,j,k,t)小于1，则执行步骤五四；

步骤五四：设在t＝0s时，枝晶尖端生长半径的初始值为10^-8m，枝晶尖端生长速度的初始值为0m/s；且枝晶尖端溶质过饱和度Ω_AG(O)(i,j,k,t)与生长贝克利特数P_V、流动贝克利特数P_U、施密特数Sc的关系为：

其中

为t时刻下计算网格(i,j,k,t)_chan＝0所对应的平均枝晶尖端生长半径，

为t时刻下计算网格(i,j,k,t)_chan＝0所对应的平均枝晶尖端生长速度；

则进行如下判断：

若

则使用5点Gauss-Legendre积分方法计算枝晶尖端溶质过饱和度Ω_AG(O)(i,j,k,t)，且积分上限取60；其中U_l(i,j,k,t)(max)为t时刻下计算网格(i,j,k,t)_chan＝0所对应的液体流动速度最大值；U_l(i,j,k,t-Δt)(max)为t-Δt时刻下计算网格(i,j,k,t)_chan＝0所对应的液体流动速度最大值；同时采用伊万卓夫函数计算相同生长贝克利特数P_V所对应的溶质过饱和度Ω_AB(M)：

为枝晶生长方向与液体流动方向的夹角，为0°～45°之间的一个随机值；调整常数a₁、a₂和a₃直至|Ω_AG(O)(i,j,k,t)-Ω_AB(M)(i,j,k,t)|≤10^-5,存储调整后的a₁、a₂和a₃；并执行步骤五五；

若

则执行步骤五五；

步骤五五：

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

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

σ^*＝1/(4π²)

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

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

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

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

<实施例1>

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

表1

表2

图2为本实施例中所使用的铸造型腔实物图。铸件型腔三维几何尺寸:0.078m(X轴)×0.252m(Y轴)×0.078m(Z轴)。铸件型腔右侧放置石墨冷铁，其三维几何尺寸：0.078m(X轴)×0.064m(Y轴)×0.078m(Z轴)。铸件型腔左侧为直浇道和横浇道。重力方向(Z轴)垂直于纸面。

图3为本实施例中熔体过热度为30℃时实验所得Sn-6wt％Pb二元合金铸件中间截面(该截面垂直于Y轴，在Y轴的截距位置0.126m)处凝固晶粒组织分布。截面上的晶粒组织形貌为柱状晶和等轴晶混合组织。

图4为本实施例对铸件温度场、流场和成分场进行三维计算所得铸件中间截面处凝固晶粒组织。晶粒组织模拟截面为S1，其垂直于Y轴且所处位置Y_s1＝0.126m。沿着Y轴截面S2和S3的位置分别为0.132m和0.12m。网格剖分尺寸:Δx₁＝Δz₁＝0.001m，Δy₁＝0.002m，Δx₂＝0.003m，Δz₂＝0.003m，Δy₂＝0.006m。针对S1截面，模拟所得凝固组织为柱状晶和等轴晶混合组织且等轴晶区存在于铸件右下方，与图3实验所得凝固组织进行对比，较好吻合。

图5为本实施例对铸件温度场、流场和成分场进行二维计算所得铸件中间截面处凝固晶粒组织。晶粒组织模拟截面为S1，其垂直于Y轴且所处位置Y_s1＝0.126m。网格剖分尺寸:Δx₁＝Δz₁＝0.001m。针对S1截面，模拟所得凝固组织为粗大柱状晶和粗大等轴晶，且柱状晶区占据面积较大，与图3实验所得凝固组织进行对比，存在较大差别。

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

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

Claims (6)

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

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

若截面S1垂直于X轴，则

若截面S1垂直于Y轴，则

步骤一二中，δ＝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₁；

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

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

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

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

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

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

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

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

[H]＝c_mpT

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

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

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

能量守恒方程为：

h_s＝c_pT

h_l＝c_pT+ΔH

其中h_s和h_l分别为固相和液相热焓，λ_l为液体导热系数，ΔH为结晶潜热，T_s固相线温度，T_l液相线温度，U_l液体流动速度,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为混合成分，C_l为平均液相成分，C_l ^*为固液界面处液相成分，C_o为合金初始成分，C_s为固相成分，k_p为平衡分配系数，m_l为液相线斜率,D_l为液相溶质扩散系数,f_s＝1-f_l为固相分数。

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

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

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

其中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时刻下温度。

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

步骤五一：将截面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)；

计算形核核心的计算方式采用基于高斯分布的形核模型，具体计算公式为：

Δ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_iocal(i,j,k,t-Δt)-T(i,j,k,t-Δt)

其中，形核核心个数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之间；形核核心个数N_num若大于0且A_ssn(i,j,k,t)小于1，则执行步骤五四；

步骤五四：设在t＝0s时，枝晶尖端生长半径的初始值为10^-8m，枝晶尖端生长速度的初始值为0m/s；且枝晶尖端溶质过饱和度Ω_AG(O)(i,j,k,t)与生长贝克利特数P_V、流动贝克利特数P_U、施密特数Sc的关系为：

Ω_AG(O)(i,j,k,t)

其中

为t时刻下计算网格(i,j,k,t)_chan＝0所对应的平均枝晶尖端生长半径，

为t时刻下计算网格(i,j,k,t)_chan＝0所对应的平均枝晶尖端生长速度；

则进行如下判断：

若

则使用5点Gauss-Legendre积分方法计算枝晶尖端溶质过饱和度Ω_AG(O)(i,j,k,t)，且积分上限取60；其中U_l(i,j,k,t)(max)为t时刻下计算网格(i,j,k,t)_chan＝0所对应的液体流动速度最大值；U_l(i,j,k,t-Δt）（max)为t-Δt时刻下计算网格(i,j,k,t)_chan＝0所对应的液体流动速度最大值；同时采用伊万卓夫函数计算相同生长贝克利特数P_V所对应的溶质过饱和度Ω_AB(M）：

为枝晶生长方向与液体流动方向的夹角，为0°～45°之间的一个随机值；调整常数a₁、a₂和a₃直至|Ω_AG(O)(i,j,k,t)-Ω_AB(M)(i,j,k,t)|≤10^-5,存储调整后的a₁、a₂和a₃；并执行步骤五五；

若

则执行步骤五五；

步骤五五：

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

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

σ^*＝1/(4π²)

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

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

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

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

Images