CN105057642A - 铸件晶粒组织形成相关数值的模拟方法 - Google Patents

Abstract

铸件晶粒组织形成相关数值的模拟方法，属于铸件晶粒组织形成模拟领域。现有的目前晶粒组织模拟中计算晶粒移动所导致的计算量大以及计算时间长的问题。一种铸件晶粒组织形成相关数值的模拟方法，进行宏观尺度网格剖分；获得t时刻下铸件内等轴晶晶核密度的分布；获得温度T分布，直到t^*时刻所有z＝0的计算网格(j，k)_char＝0中最高温度小于液相线温度T_L，并记录对应的等轴晶晶核密度采用CA法则进行晶粒组织模拟；进行晶粒生长速度的计算；设定t^*时刻后，固相均停止运动，即重复步骤二、步骤三和步骤六，直至凝固结束，输出晶粒组织。本发明准确地预测了凝固晶粒组织形成，并缩降低预测过程计算量、缩短计算时间。

Description

铸件晶粒组织形成相关数值的模拟方法

技术领域

本发明涉及一种铸件晶粒组织形成相关数值的模拟方法。

背景技术

铸件凝固过程中所形成的晶粒组织直接决定着铸件的力学性能和使用性能。例如，螺旋桨的叶片希望获得单方向生长的柱状晶或是单晶，因为该种铸件在使用过程中强调单方向性能。柱状晶组织或单晶组织具有突出的抗热冲击性能、高的高温抗蠕变力以及中温塑性，能够延长叶片的疲劳寿命，从而延长叶片的使用寿命。而汽车发动机则希望在凝固过程中获得大量细小等轴晶。由于等轴晶的晶界面积较大，杂质和缺陷分布比较分散，且等轴晶粒各向同性，这会提高铸件力学性能的均匀性。细小等轴晶还有助于提高抗疲劳性能。因此控制铸件晶粒组织是提高铸件产品力学性能的有效手段。采用实验方法探寻晶粒组织形成特点具有盲目性，而且还会增加产品成本，浪费了人力、物力和财力，也污染了环境。因此通过计算机模拟的方法预测铸件晶粒组织，根据模拟结果不断优化铸造工艺，直至得到具有所需晶粒组织的铸件。利用计算机模拟技术预测晶粒组织形成，是提高铸件产品质量的一个重要手段，也是高科技推动传统铸造产业发展的关键。

铸件晶粒组织主要呈现三种形貌：紧贴铸型内表面的细小等轴晶(晶粒轮廓线为球形)，也称激冷区；邻近激冷区的柱状晶(晶粒轮廓线近似长形)，柱状晶组成的区域称为柱状晶区；邻近柱状晶区的较粗大等轴晶(晶粒轮廓线近似圆形)，较粗大等轴晶组成的区域为等轴晶区。晶粒组织的形成是凝固过程中冶金传输的结果。凝固过程中存在四种基本传输现象：热量传输、溶质传输、液相动量传输和固相动量传输。目前广泛使用的晶粒组织模拟方法为有限元-元胞自动机法(FE-CA)和有限差分-元胞自动机法(FDM-CA)，这两种方法在预测凝固组织方面存在的问题是：当耦合固相动量传输(等轴晶粒运动)时计算非常复杂，计算时间长，这难以符合生产实际中“提高效率、缩短时间”的要求。这就强调所开发的铸件晶粒组织数值模拟方法不但包含了等轴晶粒运动对凝固组织形成的影响，还要尽量简化模型并且缩短计算时间和计算量。这在深入探寻晶粒移动对组织形成的影响规律和开发铸造工艺控制晶粒组织方面均具有重要的意义。

发明内容

本发明的目的是为了解决现有的目前晶粒组织模拟中计算晶粒移动所导致的计算量大以及计算时间长的问题，而提出一种铸件晶粒组织形成相关数值的模拟方法。

一种铸件晶粒组织形成相关数值的模拟方法，铸件晶粒组织形成相关数值的模拟方法通过以下步骤实现：

步骤一、在X×Y米铸锭系统的X方向和Y方向上，分别以X方向网格剖分步长△x和Y方向网格剖分步长△y进行宏观尺度网格剖分，形成一系列计算网格，且由(j，k)_char＝z表示计算网格的类型，z取0、4、5、6、7和8时分别铸锭网格、内冷铁网格、外冷铁网格、冒口套网格、保温材料网格和绝热材料网格；

其中，△x＝△y，△x和△y的取值范围都为1×10^-3～4×10^-3米；j、k和z均为整数，j表示二维数组中X轴方向上网格的标号，j的取值范围为1～M，k表示二维数组中Y轴方向上网格的标号，k的取值范围是1～N，

步骤二、在t时刻针对所有z＝0类型的计算网格(j，k)_char＝0，获得t时刻下铸件内等轴晶晶核密度的分布；

步骤三、针对所有z≠0类型的计算网格(j，k)_char≠0，计算能量守恒方程，获得温度T分布：

[H]＝c_pT，

其中，c_p表示比热，ρ表示密度，λ表示导热系数；

步骤四、重复步骤二和步骤三，直到在t^*时刻下所有z＝0的计算网格(j，k)_char＝0中最高温度小于液相线温度T_L，记录t^*时刻下所有计算网格(j，k)_char＝0对应的等轴晶晶核密度

步骤五、采用CA法则进行晶粒组织模拟：

将宏观尺度计算网格(j，k)_char＝0进行CA网格剖分形成M_cell·N_cell个CA网格，CA网格的标号为(j，k)_char＝0-(j_cell,k_cell)；之后将步骤四中所有计算网格(j，k)_char＝0对应的等轴晶晶核密度转化为晶核个数,并在CA网格中进行分配；其中，j_cell的取值范围是1～M_cell，k_cell的取值范围是1～N_cell， $M_{cell}=\frac{\mathrm{\Δ}x}{{\mathrm{\Δx}}_{cell}},N_{cell}=\frac{\mathrm{\Δ}y}{{\mathrm{\Δy}}_{cell}};$

步骤六、针对标号为(j_cell,k_cell)_nuclei的CA网格通过生长速度计算公式：

进行晶粒生长速度的计算，其中：D₁表示液相溶质扩散率，单位为m²s^-1，m_l表示液相线斜率，Kwt％^-1，q表示平衡分配系数，C_o表示合金初始成分；Γ表示吉布斯-汤姆森系数，T_L表示液相线温度；

步骤七、设定t^*时刻后，固相均停止运动，即重复步骤二、步骤三和步骤六，直至凝固结束，输出晶粒组织。

本发明的有益效果为：

本发明是设计了铸件中晶粒组织模拟的方法。由于固相运动是凝固过程中固有的一种物理现象，在预测凝固晶粒组织形成过程中，考虑固相移动对晶粒组织的影响，分别针对铸锭系统中铸锭网格、内冷铁网格、外冷铁网格、冒口套网格、保温材料网格和绝热材料网格等不同类型网格进行能量守恒方程、成分守恒方程、等轴晶晶核密度、固相分数守恒方程、液相动量守恒方程、固相动量守恒方程和质量守恒方程的计算，获得铸件内等轴晶晶核密度的分布、获得温度T分布，然后采用CA法则进行晶粒组织模拟，计算标号为(j_cell,k_cell)_nuclei的CA网格计算生长速度，直至凝固结束，输出晶粒组织。从而更加准确地预测了凝固晶粒组织形成，相比目前晶粒组织预测中在考虑固相移动所需的计算时间而言，缩短了60-70％，且将计算量降低80％左右，解决目前晶粒组织预测计算时间长和计算量大的问题。

另外，本发明适用于各类尺寸的砂型和金属型中晶粒组织的预测。利用本发明可以更为准确的预测晶粒组织形成，从多个方面为工艺改进提供帮助，市场应用潜力巨大，一旦被广泛采用，将有百万元以上的产值。

附图说明

图1为本发明涉及的铸造系统；

图2为与实施例1形成对比的凝固实验所得凝固晶粒组织；

图3为本发明实施例1中采用本方法，考虑固相移动模拟所得晶粒组织；

图4为本发明与施例1形成对比的不考虑固相移动模拟所得晶粒组织；

图5为本发明流程图。

具体实施方式

具体实施方式一：

本实施方式的铸件晶粒组织形成相关数值的模拟方法，铸件晶粒组织形成相关数值的模拟方法通过以下步骤实现：

步骤一、铸锭系统的大小为X×Y米，在此铸锭系统的X方向和Y方向上，分别以X方向网格剖分步长△x和Y方向网格剖分步长△y进行宏观尺度网格剖分，形成一系列计算网格，且由(j，k)_char＝z表示计算网格的类型，z取0、4、5、6、7和8时分别铸锭网格、内冷铁网格、外冷铁网格、冒口套网格、保温材料网格和绝热材料网格；

其中，△x＝△y，△x和△y的取值范围都为1×10^-3～4×10^-3米；j、k和z均为整数，j表示二维数组中X轴方向上网格的标号，j的取值范围为1～M，k表示二维数组中Y轴方向上网格的标号，k的取值范围是1～N，

步骤二、在t时刻针对所有z＝0类型的计算网格(j，k)_char＝0，获得t时刻下铸件内等轴晶晶核密度的分布；

步骤三、针对所有z≠0类型的计算网格(j，k)_char≠0，计算能量守恒方程，获得温度T分布：

[H]＝c_pT，

其中，c_p表示比热，ρ表示密度，λ表示导热系数；

步骤四、重复步骤二和步骤三，直到在t^*时刻下所有z＝0的计算网格(j，k)_char＝0中最高温度小于液相线温度T_L，记录t^*时刻下所有计算网格(j，k)_char＝0对应的等轴晶晶核密度

步骤五、采用CA法则进行晶粒组织模拟：

将宏观尺度计算网格(j，k)_char＝0进行CA网格剖分形成M_cell·N_cell个CA网格，CA网格的标号为(j，k)_char＝0-(j_cell,k_cell)；之后将步骤四中所有计算网格(j，k)_char＝0对应的等轴晶晶核密度转化为晶核个数,并在CA网格中进行分配；其中，j_cell的取值范围是1～M_cell，k_cell的取值范围是1～N_cell， $M_{cell}=\frac{\mathrm{\Δ}x}{{\mathrm{\Δx}}_{cell}},N_{cell}=\frac{\mathrm{\Δ}y}{{\mathrm{\Δy}}_{cell}};$

步骤六、针对标号为(j_cell,k_cell)_nuclei的CA网格通过生长速度计算公式：

进行晶粒生长速度的计算，其中：D₁表示液相溶质扩散率，单位为m²s^-1，m_l表示液相线斜率，Kwt％^-1，q表示平衡分配系数，C_o表示合金初始成分；Γ表示吉布斯-汤姆森系数，T_L表示液相线温度；CA法则技术是指采用四邻居法则进行网格捕捉，体现晶粒的生长；

步骤七、设定t^*时刻后，固相均停止运动，即重复步骤二、步骤三和步骤六，直至凝固结束，输出晶粒组织。

具体实施方式二：

与具体实施方式一不同的是，本实施方式的铸件晶粒组织形成相关数值的模拟方法，步骤二所述获得t时刻下铸件内等轴晶晶核密度的分布的过程为，

步骤二一、计算能量守恒方程，

h_s＝c_PT，

h_l＝c_PT+△H，

[H]＝f_sh_s+(1-f_s)h₁，求得铸锭内温度T的分布：

其中，h_s表示固相热焓，h_l表示液相热焓，c_p表示合金比热，[H]表示混合热焓，ρ表示合金密度，λ表示合金导热系数，△H表示结晶潜热，T表示铸锭网格内部温度，f_s表示固相分数，f₁表示液相分数f₁＝1-f_s，t表示时间，表示固相移动速度，表示液相流动速度，和在0s时的值均为0ms^-1；

步骤二二、计算成分守恒方程：

$\frac{\∂\left(f_s{C}_s\right)}{\∂t}+\▿\·\left(f_s{C}_s\stackrel{\→}{v_s}\right)=0,$

求得铸件内平均成分[C_mix]的分布：

[C_mix]＝f_sC_s+f_lC_l

其中：C_s表示固相成分，C₁表示液相成分，[C_mix]表示平均成分；

步骤二三、计算固相分数守恒方程：

$\frac{\∂\left(f_s\right)}{\∂t}+\▿\·\left(f_s\stackrel{\→}{v_s}\right)=0;$

步骤二四、计算等轴晶晶核密度守恒方程：

首先，判断在t时刻下，z＝0类型的计算网格(j，k)_char＝0是否存在等轴晶形核，并针对z＝0类型的计算网格(j，k)_char＝0计算通过过冷度根据公式ΔT(j,k)^t＝T_l-T(j,k)^t，若ΔT(j,k)^t>ΔT(j,k)^t-Δt，则表明计算网格(j，k)_char＝0中存在等轴晶形核现象，形核密度计算基于高斯形核分布公式：

$N_{nuclei}^t=N_{nuclei}^{t-\mathrm{\Δ}t}+\[\mathrm{\Δ}T{(j,k)}^t-\mathrm{\Δ}T{(j,k)}^{t-\mathrm{\Δ}t}\]\×\left(\frac{N_{max}}{{\mathrm{\ΔT}}_{\mathrm{\σ}}\sqrt{2\mathrm{\π}}}\right)e^{\[-\frac{1}{2}{\left(\frac{\mathrm{\Δ}T{(j,k)}^t-{\mathrm{\ΔT}}_{max}}{{\mathrm{\ΔT}}_{\mathrm{\σ}}}\right)}^2\]};$

若ΔT(j,k)^t≤ΔT(j,k)^t-Δt，则表明计算网格(j，k)_char＝0中没有发生等轴晶形核现象；

其次，计算轴晶晶核密度守恒方程：

其中：ΔT(j,k)^t表示计算网格(j，k)_char＝0在t时刻下的过冷度，单位为℃；ΔT(j,k)^t-Δt表示计算网格(j，k)_char＝0在t-△t时刻下的过冷度，单位为℃；△t表示时间步长；表示计算网格(j，k)_char＝0在t时刻下的等轴晶晶核密度，单位为m^-3；表示计算网格(j，k)_char＝0在t-△t时刻下的等轴晶晶核密度；N_max为最大等轴晶形核密度，△T_σ表示标准方差形核过冷度，△T_max表示最大等轴晶形核过冷度；

若一个计算网格(j，k)_char＝0，在t时刻下则进行如下处理：

$N_{nuclei}^t=0,$

$f_s^t=0,$

$C_s^t=0;$

步骤二五、求解固相分数f_s的一元二次方程：aa·f_s ²+bb·f_s+cc＝0，计算固相分数f_s；

其中：aa、bb和cc为一元二次方程的系数：

aa＝ρ·△H·(1-k_p)；

bb＝ρ△H(k_p-1)(1+f_s)+(k_p-1)(ρc_PT_M-[H])-(2△tD₁)(N_nucleiπd_P)ρ△H；

cc＝ρ·△H·(1-k_p)f_s+(ρc_PT_M-[H])(1-k_p)f_s

+(2△tD_l)(N_nucleiπd_p)(ρc_pT_M-[H]-C_lm_lρc_p)^；

其中：D_l表示液相中溶质扩散系数；d_p表示等轴晶晶核直径，初始时刻0s时d_p的初始值为10^-6m；T_M表示合金熔点；m_l表示液相线斜率；k_p表示溶质平衡分配系数；

步骤二六、分别计算凝固引起的温度、液相成分和固相成分的变化：

(1)凝固引起的温度变化：

$T^t=\frac{\[H\]-(1-f_s)\mathrm{\Δ}H}{c_p}$

(2)凝固引起的液相成分变化：

$C_1^t=\frac{C_1^{t-\mathrm{\Δ}t}(1-f_s^{t-\mathrm{\Δ}t})-\left(\frac{T_M-T^t}{m_1}{k}_p\right)(f_s^t-f_s^{t-\mathrm{\Δ}t})}{(1-f_s^t)}$

(3)凝固引起的固相成分变化：

$C_s^t=\frac{C_s^{t-\mathrm{\Δ}t}{f}_s^{t-\mathrm{\Δ}t}+\left(\frac{T_M-T^t}{m_1}{k}_p\right)(f_s^t-f_s^{t-\mathrm{\Δ}t})}{f_s^t};$

步骤二七、计算等轴晶晶核直径：

针对所有z＝0的计算网格(j，k)_char＝0，判断在t时刻是否存在等轴晶晶核密度：

若则通过计算等轴晶晶核直径；

否则，说明不存在晶核，不需计算该网格对应的等轴晶晶核直径；

步骤二八、分别计算固相动量守恒方程和液相动量守恒方程：

(1)计算固相动量守恒方程：

针对所有z＝0的计算网格(j，k)_char＝0，若计算网格(j，k)_char＝0的固相分数小于临界固相分数，即：则计算如下的固相动量守恒方程：

$\frac{\∂f_s{u}_s}{\∂t}+\▿\·\left(f_s{u}_s{u}_s\right)=-\frac{f_s}{\mathrm{\ρ}}\▿P+f_s{\mathrm{\μ}}_s\▿\·(\▿u_s)+\mathrm{\ρ}(1+{\mathrm{\β}}_{s1})\stackrel{\→}{G}+K_{1s}(u_1-u_s);$

否则，计算网格(j，k)_char＝0对应的固相速度为0ms^-1；

(2)液相动量守恒方程：

通过如下的液相动量守恒方程进行计算，

$\frac{\∂f_1{u}_1}{\∂t}+\▿\·\left(f_1{u}_1{u}_1\right)=-\frac{f_1}{\mathrm{\ρ}}\▿P+f_1{\mathrm{\μ}}_1\▿\·(\▿u_1)+\mathrm{\ρ}\[1+{\mathrm{\β}}_T(T_{ref}-T)+{\mathrm{\β}}_C(C_{ref}-C_1)\]\stackrel{\→}{G}+K_{1s}(u_s-1)$

其中： ${\mathrm{\μ}}_s=\frac{{\[1-f_s/(f_s^c+{10}^{-5})\]}^{-2.5f_s^c}-f_s}{f_s}{\mathrm{\μ}}_1,$

$K_{1s}=\frac{3f_s{\mathrm{\ρW}}_d}{4d_p}\left|v_1-v_s\right|,$

$W_d=\frac{48W_{ke}{f}_s}{\mathrm{Re}}+W_{ie},$

$\mathrm{Re}=\frac{(1-f_s){\mathrm{\ρd}}_p}{{\mathrm{\μ}}_1}\left|v_1-v_s\right|,$

$W_{ke}=\frac{1}{2}\[\frac{{(1-f_s)}^3}{f_s}\]\frac{1+4.7f_s}{1+1.83f_s},$

$W_{ie}=\frac{24({10}^{EE}-1)}{\mathrm{Re}{\[1-0.9{(0.75-f_s)}^{1/3}{f}_s^{2/3}\]}^3},$

$EE=0.261{\mathrm{Re}}^{0.369}-0.105{\mathrm{Re}}^{0.431}-\frac{0.124}{1+{\left({\log }_{10}\mathrm{Re}\right)}^2};$ 式中，

μ_s表示固相粘度，μ_l表示液相粘度，P表示压力，β_sl表示合金的凝固收缩系数，表示重力加速度，K_ls表示固液拖拽系数，β_T和β_C分别表示热膨胀系数和溶质膨胀系数，T_ref和C_ref分别表示参考温度和参考成分，表示临界固相分数，Re表示雷诺数，W_d、W_ke、W_ie和EE均表示中间变量，没有实际物理含义，只具有数学含义；

步骤二九、计算质量守恒方程，检验动量方程求解是否准确。

具体实施方式三：

与具体实施方式一或二不同的是，本实施方式的铸件晶粒组织形成相关数值的模拟方法，步骤二九所述检验动量方程求解是否准确的过程为，通过判断是否成立来检验动量方程求解是否准确，若成立，则动量守恒方程求解准确，若不成立，则动量方程求解不准确，重新求解动量守恒方程。

具体实施方式四：

与具体实施方式三不同的是，本实施方式的铸件晶粒组织形成相关数值的模拟方法，步骤五所述将宏观尺度计算网格(j，k)_char＝0进行CA网格剖分形成M_cell·N_cell个CA网格的过程为，选取CA网格剖分尺寸，X方向和Y方向的剖分步长分别为△x_cell和△y_cell，且△x_cell＝△y_cell，在X方向上，网格剖分步长△x是CA网格剖分步长△x_cell的整数倍，△x≥5△x_cell，且在Y方向上，网格剖分步长△y是CA网格剖分步长△y_cell的整数倍，且△y≥5△y_cell；CA网格的标号为(j，k)_char＝0-(j_cell,k_cell)；

其中，(j，k)_char＝0表示z＝0的宏观尺度剖分所得的计算网格；j_cell和k_cell均为整数，j_cell的取值范围是1～M_cell，k_cell的取值范围是1～N_cell，每一个计算网格(j，k)_char＝0中都包含M_cell·N_cell个CA网格。

具体实施方式五：

与具体实施方式一、二或四不同的是，本实施方式的铸件晶粒组织形成相关数值的模拟方法，步骤五所述将步骤四中所有计算网格(j，k)_char＝0对应的等轴晶晶核密度转化为晶核个数,并在CA网格中进行分配的过程为，

第一，针对计算网格(j，k)_char＝0，将晶核密度N_nuclei转化为晶核个数Num_nuclei；

第二，将上一步得到的计算网格(j，k)_char＝0对应的Num_nuclei个晶核，分配到当前计算网格(j，k)_char＝0包含的标号为(j,k)_char＝0-(j_cell,k_cell)的CA网格中。

具体实施方式六：

与具体实施方式五不同的是，本实施方式的铸件晶粒组织形成相关数值的模拟方法，将晶核密度N_nuclei转化为晶核个数Num_nuclei的过程为，通过如下转化公式：

将计算网格(j，k)_char＝0的晶核密度N_nuclei转化为晶核个数Num_nuclei。

具体实施方式七：

与具体实施方式六不同的是，本实施方式的铸件晶粒组织形成相关数值的模拟方法，计算网格(j，k)_char＝0对应的Num_nuclei个晶核，分配到当前计算网格(j，k)_char＝0包含的标号为(j,k)_char＝0-(j_cell,k_cell)的CA网格中的过程为：

若M_cell·N_cell≤Num_nuclei，则计算网格(j，k)_char＝0包含的M_cell·N_cell个CA网格均为形核核心：

若M_cell·N_cell＞Num_nuclei，则随机选取一个CA网格放置一个晶核，且该CA网格的标号变为(j_cell,k_cell)_nuclei，此时剩余晶核个数为Num_nuclei-1；并重复此操作直至剩余晶核个数为0。

实施例1：

如图1所示的铸造系统结构示意图，并通过表1列出Al-4.0wt％Cu合金热物性参数以及高斯形核分布，和表2列出的砂型和石墨冷铁涉及的物理特性及其对应的特性值：

表1Al-4.0wt％Cu合金热物性参数以及高斯形核分布公式

表2砂型和石墨冷铁涉及的物理特性及其对应的特性值

	砂型	石墨冷铁
			ρ(kg/m3)	1520.0	2250.00
λ(W/mK)	0.73	140.0
			cp(J/kgK)	1070.0	710.0

利用本发明方法模拟所得晶粒组织的形成：

步骤一、在X×Y米铸锭系统的X方向和Y方向上，分别以X方向网格剖分步长△x和Y方向网格剖分步长△y进行宏观尺度网格剖分，形成一系列计算网格，且由(j，k)_char＝z表示计算网格的类型，z取0、4、5、6、7和8时分别铸锭网格、内冷铁网格、外冷铁网格、冒口套网格、保温材料网格和绝热材料网格；

其中，△x＝△y，△x和△y的取值范围都为1×10^-3～4×10^-3米；j、k和z均为整数，j表示二维数组中X轴方向上网格的标号，j的取值范围为1～M，k表示二维数组中Y轴方向上网格的标号，k的取值范围是1～N，

步骤二、在t时刻针对所有z＝0类型的计算网格(j，k)_char＝0，获得t时刻下铸件内等轴晶晶核密度的分布；

步骤三、针对所有z≠0类型的计算网格(j，k)_char≠0，计算能量守恒方程，获得温度T分布：

[H]＝c_pT，

其中，c_p表示比热，ρ表示密度，λ表示导热系数；

步骤四、重复步骤二和步骤三，直到在t^*时刻下所有z＝0的计算网格(j，k)_char＝0中最高温度小于液相线温度T_L，记录t^*时刻下所有计算网格(j，k)_char＝0对应的等轴晶晶核密度

步骤五、采用CA法则进行晶粒组织模拟：

将宏观尺度计算网格(j，k)_char＝0进行CA网格剖分形成M_cell·N_cell个CA网格，CA网格的标号为(j，k)_char＝0-(j_cell,k_cell)；之后将步骤四中所有计算网格(j，k)_char＝0对应的等轴晶晶核密度转化为晶核个数,并在CA网格中进行分配；其中，j_cell的取值范围是1～M_cell，k_cell的取值范围是1～N_cell， $M_{cell}=\frac{\mathrm{\Δ}x}{{\mathrm{\Δx}}_{cell}},N_{cell}=\frac{\mathrm{\Δ}y}{{\mathrm{\Δy}}_{cell}};$

步骤六、针对标号为(j_cell,k_cell)_nuclei的CA网格通过生长速度计算公式：

进行晶粒生长速度的计算，其中：D₁表示液相溶质扩散率，单位为m²s^-1，m_l表示液相线斜率，Kwt％^-1，q表示平衡分配系数，C_o表示合金初始成分；Γ表示吉布斯-汤姆森系数，T_L表示液相线温度；

步骤七、设定t^*时刻后，固相均停止运动，即重复步骤二、步骤三和步骤六，直至凝固结束，输出晶粒组织。

(1)自浇道浇筑并凝固后，截去浇道位置的凝固部分保留铸件部分，为方便操作，将铸件截成两部分进行砂纸打磨、化学腐蚀剂浸泡等操作，完毕后再将截后的两部分拼在一起形成完整的铸件，用于实验对照，如图2。

(2)对比附图2和附图3，得出：采用本方法考虑固相移动模拟所得晶粒组织与实验观察吻合较好。本模拟可以再现靠近石墨冷铁附近的柱状晶组织和铸件顶部形成的柱状晶组织，即虚线内的区域，以及铸件中部和底部的等轴晶组织。这是因为Al-4wt％Cu合金凝固过程中，固相密度大于液相密度(凝固收缩系数β_sl>0)，因此在重力场作用下，固相向铸件底部沉淀。形成的等轴晶晶核向铸件下部运动，导致铸件上部等轴晶晶核数量减少，柱状晶生长过程中没有等轴晶对其进行阻碍，因此铸件顶部表现为柱状晶组织，而铸件中部和下部为等轴晶组织。

对比附图2和附图4，得出：不考虑固相移动模拟所得晶粒组织无法再现实验结果。除了靠近石墨冷铁附近的柱状晶组织，铸件其它位置均为等轴晶组织。

Claims (7)

1.一种铸件晶粒组织形成相关数值的模拟方法，其特征在于：铸件晶粒组织形成相关数值的模拟方法通过以下步骤实现：

步骤一、在X×Y米铸锭系统的X方向和Y方向上，分别以X方向网格剖分步长△x和Y方向网格剖分步长△y进行宏观尺度网格剖分，形成一系列计算网格，且由(j，k)_char＝z表示计算网格的类型，z取0、4、5、6、7和8时分别铸锭网格、内冷铁网格、外冷铁网格、冒口套网格、保温材料网格和绝热材料网格；

其中，△x＝△y，△x和△y的取值范围都为1×10^-3～4×10^-3米；j、k和z均为整数，j表示二维数组中X轴方向上网格的标号，j的取值范围为1～M，k表示二维数组中Y轴方向上网格的标号，k的取值范围是1～N，

步骤二、在t时刻针对所有z＝0类型的计算网格(j，k)_char＝0，获得t时刻下铸件内等轴晶晶核密度的分布；

步骤三、针对所有z≠0类型的计算网格(j，k)_char≠0，计算能量守恒方程，获得温度T分布：

[H]＝c_pT，

其中，c_p表示比热，ρ表示密度，λ表示导热系数；

步骤四、重复步骤二和步骤三，直到在t^*时刻下所有z＝0的计算网格(j，k)_char＝0中最高温度小于液相线温度T_L，记录t^*时刻下所有计算网格(j，k)_char＝0对应的等轴晶晶核密度

步骤五、采用CA法则进行晶粒组织模拟：

将宏观尺度计算网格(j，k)_char＝0进行CA网格剖分形成M_cell·N_cell个CA网格，CA网格的标号为(j，k)_char＝0-(j_cell,k_cell)；之后将步骤四中所有计算网格(j，k)_char＝0对应的等轴晶晶核密度转化为晶核个数,并在CA网格中进行分配；其中，j_cell的取值范围是1～M_cell，k_cell的取值范围是1～N_cell， $M_{cell}=\frac{\mathrm{\Δ}x}{{\mathrm{\Δx}}_{cell}},$ $N_{cell}=\frac{\mathrm{\Δ}y}{{\mathrm{\Δy}}_{celll}};$

步骤六、针对标号为(j_cell,k_cell)_nuclei的CA网格通过生长速度计算公式：

进行晶粒生长速度的计算，其中：D₁表示液相溶质扩散率，单位为m²s^-1，m_l表示液相线斜率，Kwt％^-1，q表示平衡分配系数，C_o表示合金初始成分；Γ表示吉布斯-汤姆森系数，T_L表示液相线温度；

步骤七、设定t^*时刻后，固相均停止运动，即重复步骤二、步骤三和步骤六，直至凝固结束，输出晶粒组织。

2.根据权利要求1所述铸件晶粒组织形成相关数值的模拟方法，其特征在于：步骤二所述获得t时刻下铸件内等轴晶晶核密度的分布的过程为，

步骤二一、计算能量守恒方程，

h_s＝c_PT，

h_l＝c_PT+△H，

[H]＝f_sh_s+(1-f_s)h₁，求得铸锭内温度T的分布：

其中，h_s表示固相热焓，h_l表示液相热焓，c_p表示合金比热，[H]表示混合热焓，ρ表示合金密度，λ表示合金导热系数，△H表示结晶潜热，T表示铸锭网格内部温度，f_s表示固相分数，f₁表示液相分数f₁＝1-f_s，t表示时间，表示固相移动速度，表示液相流动速度，和在0s时的值均为0ms^-1；

步骤二二、计算成分守恒方程：

$\frac{\∂\left(f_s{C}_s\right)}{\∂t}+\▿\·\left(f_s{C}_s\stackrel{\→}{V_s}\right)=0,$

求得铸件内平均成分[C_mix]的分布：

[C_mix]＝f_sC_s+f_lC_l

其中：C_s表示固相成分，C₁表示液相成分，[C_mix]表示平均成分；

步骤二三、计算固相分数守恒方程：

$\frac{\∂\left(f_s\right)}{\∂t}+\▿\·\left(f_s\stackrel{\→}{v_s}\right)=0;$

步骤二四、计算等轴晶晶核密度守恒方程：

首先，判断在t时刻下，z＝0类型的计算网格(j，k)_char＝0是否存在等轴晶形核，并针对z＝0类型的计算网格(j，k)_char＝0计算通过过冷度根据公式ΔT(j,k)^t＝T_l-T(j,k)^t，若ΔT(j,k)^t>ΔT(j,k)^t-Δt，则表明计算网格(j，k)_char＝0中存在等轴晶形核现象，形核密度计算基于高斯形核分布公式：

$N_{nuclei}^t=N_{nuclei}^{t-\mathrm{\Δ}t}+\[\mathrm{\Δ}T{(j,k)}^t-\mathrm{\Δ}T{(j,k)}^{t-\mathrm{\Δ}t}\]\×\left(\frac{N_{max}}{{\mathrm{\ΔT}}_{\mathrm{\σ}}\sqrt{2\mathrm{\π}}}\right)e^{\[-\frac{1}{2}{\left(\frac{\mathrm{\Δ}T{(j,k)}^t-{\mathrm{\ΔT}}_{max}}{{\mathrm{\ΔT}}_{\mathrm{\σ}}}\right)}^2\]};$

若则表明计算网格(j，k)_char＝0中没有发生等轴晶形核现象；

其次，计算轴晶晶核密度守恒方程：

其中：ΔT(j,k)^t表示计算网格(j，k)_char＝0在t时刻下的过冷度，单位为℃；ΔT(j,k)^t-Δt表示计算网格(j，k)_char＝0在t-△t时刻下的过冷度，单位为℃△t表示时间步长；表示计算网格(j，k)_char＝0在t时刻下的等轴晶晶核密度，单位为m^-3；表示计算网格(j，k)_char＝0在t-△t时刻下的等轴晶晶核密度；N_max为最大等轴晶形核密度，△T_σ表示标准方差形核过冷度，△T_max表示最大等轴晶形核过冷度；

若一个计算网格(j，k)_char＝0，在t时刻下则进行如下处理：

$N_{nuclei}^t=0,$

$f_s^t=0,$

$C_s^t=0;$

步骤二五、求解固相分数f_s的一元二次方程：aa·f_s ²+bb·f_s+cc＝0，计算固相分数f_s；

其中：aa、bb和cc为一元二次方程的系数：

aa＝ρ·△H·(1-k_p)；

bb＝ρ△H(k_p-1)(1+f_s)+(k_p-1)(ρc_PT_M-[H])-(2△tD₁)(N_nucleiπd_P)ρ△H；

cc＝ρ·△H·(1-k_p)f_s+(ρc_PT_M-[H])(1-k_p)f_s

；

+(2△tD_l)(N_nucleiπd_p)(ρc_pT_M-[H]-C_lm_lρ_cp)

其中：D_l表示液相中溶质扩散系数；d_p表示等轴晶晶核直径，初始时刻0s时d_p的初始值为10^-6m；T_M表示合金熔点；m_l表示液相线斜率；k_p表示溶质平衡分配系数；

步骤二六、分别计算凝固引起的温度、液相成分和固相成分的变化：

(1)凝固引起的温度变化：

$T^t=\frac{\[H\]-(1-f_s)\mathrm{\Δ}H}{c_p}$

(2)凝固引起的液相成分变化：

$C_1^t=\frac{C_1^{t-\mathrm{\Δ}t}(1-f_s^{t-\mathrm{\Δ}t})-\left(\frac{T_M-T^t}{m_1}{k}_p\right)(f_s^t-f_s^{t-\mathrm{\Δ}t})}{(1-f_s^t)}$

(3)凝固引起的固相成分变化：

$C_s^t=\frac{C_s^{t-\mathrm{\Δ}t}{f}_s^{t-\mathrm{\Δ}t}+\left(\frac{T_M-T^t}{m_1}{k}_p\right)(f_s^t-f_s^{t-\mathrm{\Δ}t})}{f_s^t};$

步骤二七、计算等轴晶晶核直径：

针对所有z＝0的计算网格(j，k)_char＝0，判断在t时刻是否存在等轴晶晶核密度：

若 $N_{nuclei}^t>0,$ 则通过 $d_p^t={\left(\frac{3f_s^t}{4{\mathrm{\πN}}_{nuclei}^t}\right)}^{\frac{1}{3}}$ 计算等轴晶晶核直径；

否则，说明不存在晶核，不需计算该网格对应的等轴晶晶核直径；

步骤二八、分别计算固相动量守恒方程和液相动量守恒方程：

(1)计算固相动量守恒方程：

针对所有z＝0的计算网格(j，k)_char＝0，若计算网格(j，k)_char＝0的固相分数小于临界固相分数，即：则计算如下的固相动量守恒方程：

$\frac{\∂f_s{u}_s}{\∂t}+\▿\·\left(f_s{u}_s{u}_s\right)=-\frac{f_s}{\mathrm{\ρ}}\▿P+f_s{\mathrm{\μ}}_s\▿\·(\▿u_s)+\mathrm{\ρ}(1+{\mathrm{\β}}_{s1})\stackrel{\→}{G}+K_{1s}(u_1-u_s);$

否则，计算网格(j，k)_char＝0对应的固相速度为0ms^-1；

(2)液相动量守恒方程：

通过如下的液相动量守恒方程进行计算，

$\frac{\∂f_1{u}_1}{\∂t}+\▿\·\left(f_1{u}_1{u}_1\right)=-\frac{f_1}{\mathrm{\ρ}}\▿P+f_1{\mathrm{\μ}}_1\▿\·(\▿u_1)+\mathrm{\ρ}\[1+{\mathrm{\β}}_T(T_{ref}-T)+{\mathrm{\β}}_C(C_{ref}-C_1)\]\stackrel{\→}{G}+K_{1s}(u_s-1);$

其中： ${\mathrm{\μ}}_s=\frac{{\[1-f_s/(f_s^c+{10}^{-5})\]}^{-2.5f_s^c}-f_s}{f_s}{\mathrm{\μ}}_1,$

$K_{1s}=\frac{3f_s{\mathrm{\ρW}}_d}{4d_p}\left|v_1-v_s\right|,$

$W_d=\frac{48W_{ke}{f}_s}{\mathrm{Re}}+W_{ie},$

$\mathrm{Re}=\frac{(1-f_s){\mathrm{\ρd}}_p}{{\mathrm{\μ}}_1}\left|v_1-v_s\right|,$

$W_{ke}=\frac{1}{2}\[\frac{{(1-f_s)}^3}{f_s}\]\frac{1+4.7f_s}{1+1.83f_s},$

$W_{ie}=\frac{24({10}^{EE}-1)}{\mathrm{Re}{\[1-0.9{(0.75-f_s)}^{1/3}{f}_s^{2/3}\]}^3},$

$EE=0.261{\mathrm{Re}}^{0.369}-0.105{\mathrm{Re}}^{0.431}-\frac{0.124}{1+{\left({\log }_{10}\mathrm{Re}\right)}^2};$ 式中，

μ_s表示固相粘度，μ_l表示液相粘度，P表示压力，β_sl表示合金的凝固收缩系数，表示重力加速度，K_ls表示固液拖拽系数，β_T和β_C分别表示热膨胀系数和溶质膨胀系数，T_ref和C_ref分别表示参考温度和参考成分，表示临界固相分数，Re表示雷诺数，W_d、W_ke、W_ie和EE均表示中间变量；

步骤二九、计算质量守恒方程，检验动量方程求解是否准确。

3.根据权利要求1或2所述铸件晶粒组织形成相关数值的模拟方法，其特征在于：步骤二九所述检验动量方程求解是否准确的过程为，通过判断是否成立来检验动量方程求解是否准确，若成立，则动量守恒方程求解准确，若不成立，则动量方程求解不准确，重新求解动量守恒方程。

4.根据权利要求3所述铸件晶粒组织形成相关数值的模拟方法，其特征在于：步骤五所述将宏观尺度计算网格(j，k)_char＝0进行CA网格剖分形成M_cell·N_cell个CA网格的过程为，选取CA网格剖分尺寸，X方向和Y方向的剖分步长分别为△x_cell和△y_cell，且△x_cell＝△y_cell，在X方向上，网格剖分步长△x是CA网格剖分步长△x_cell的整数倍，△x≥5△x_cell，且在Y方向上，网格剖分步长△y是CA网格剖分步长△y_cell的整数倍，且△y≥5△y_cell；CA网格的标号为(j，k)_char＝0-(j_cell,k_cell)；

其中，(j，k)_char＝0表示z＝0的宏观尺度剖分所得的计算网格；j_cell和k_cell均为整数，j_cell的取值范围是1～M_cell，k_cell的取值范围是1～N_cell， 每一个计算网格(j，k)_char＝0中都包含M_cell·N_cell个CA网格。

5.根据权利要求1、2或4所述铸件晶粒组织形成相关数值的模拟方法，其特征在于：步骤五所述将步骤四中所有计算网格(j，k)_char＝0对应的等轴晶晶核密度转化为晶核个数,并在CA网格中进行分配的过程为，

第一，针对计算网格(j，k)_char＝0，将晶核密度N_nuclei转化为晶核个数Num_nuclei；

第二，将上一步得到的计算网格(j，k)_char＝0对应的Num_nuclei个晶核，分配到当前计算网格(j，k)_char＝0包含的标号为(j,k)_char＝0-(j_cell,k_cell)的CA网格中。

6.根据权利要求5所述铸件晶粒组织形成相关数值的模拟方法，其特征在于：将晶核密度N_nuclei转化为晶核个数Num_nuclei的过程为，通过如下转化公式：

将计算网格(j，k)_char＝0的晶核密度N_nuclei转化为晶核个数Num_nuclei。

7.根据权利要求6所述铸件晶粒组织形成相关数值的模拟方法，其特征在于：计算网格(j，k)_char＝0对应的Num_nuclei个晶核，分配到当前计算网格(j，k)_char＝0包含的标号为(j,k)_char＝0-(j_cell,k_cell)的CA网格中的过程为：

若M_cell·N_cell≤Num_nuclei，则计算网格(j，k)_char＝0包含的M_cell·N_cell个CA网格均为形核核心：

若M_cell·N_cell＞Num_nuclei，则随机选取一个CA网格放置一个晶核，且该CA网格的标号变为(j_cell,k_cell)_nuclei，此时剩余晶核个数为Num_nuclei-1；并重复此操作直至剩余晶核个数为0。