CN102289542A - 一种铸件宏观偏析数值模拟的方法 - Google Patents

Abstract

一种铸件宏观偏析数值模拟的方法，它涉及一种铸件中宏观偏析模拟的方法。针对目前宏观偏析计算时间长、精度低且不适合预测大尺寸或形状复杂的铸件中宏观偏析形成的问题。所述方法包括如下步骤：选取枝晶生长微观尺度计算域，进行网格剖分；在微观尺度上计算不同最大形核密度下合金凝固时的枝晶生长形貌，输出固相分数-平均枝晶固相特征直径曲线；对铸件进行宏观尺度网格剖分；在宏观尺度上计算质量、动量、能量和成分守恒方程；采用牛顿-下山法计算固相分数；采用线性插值技术，宏观尺度上糊状区渗透率模型的计算基于给定形核密度下的固相分数-平均枝晶固相特征直径曲线。本发明用于各类尺寸和复杂形状的砂型铸造及金属型铸造中宏观偏析的数值预测。

Description

一种铸件宏观偏析数值模拟的方法

技术领域

本发明涉及一种铸件中宏观偏析模拟的方法。

背景技术

铸件的宏观偏析主要是指铸件内部的化学成分不均匀，并由此导致的铸件内部在组织结构和性能方面的不均匀性。这种不均匀性降低了铸件的冲击韧性和塑性，增加了热裂倾向，有时还导致铸件难于加工或铸件废品率高。对于在空气中或在腐蚀介质中工作的铸件，宏观偏析部位更容易遭受腐蚀破坏。因此利用计算机模拟的方法准确预测铸件中宏观偏析形成，为科学确定铸造工艺并寻找减小宏观偏析程度的工艺措施提供参数选择和理论指导。

合金的凝固是一个包含着宏观范畴(10^-2-1m)热、溶质、动量传输与微观范畴(10^-4-10^-2m)晶粒形核长大的复杂物理过程。对宏观尺度上宏观偏析形成行为的准确描述需要微观组织的信息(二次枝晶臂间距SDAS)。目前广泛使用的经典混合理论的缺点是因忽略凝固组织变化而过高估计枝晶间液体流动速度；固-液两相流模型考虑了凝固组织变化对宏观偏析的影响，但存在的问题是没有考虑组织形成过程中的随机现象(随机形核、随机结晶取向)，多个变量之间需要进行高强度的迭代耦合求解，因此计算时间长，不适合预测大尺寸或形状复杂的铸件中宏观偏析形成。这就要求所开发的宏观偏析数值模拟方法不但要考虑微观组织形成对宏观偏析的影响，还要具有较高的计算效率，这在实际应用方面均具有重要意义。

发明内容

本发明的目的是提供一种铸件宏观偏析数值模拟的方法，以解决目前宏观偏析计算时间长、精度低且不适合预测大尺寸或形状复杂的铸件中宏观偏析形成的问题。

本发明为解决上述技术问题采取的技术方案是：所述方法包括如下步骤：

步骤一、将枝晶生长计算域进行微观尺度网格剖分，剖分为n×m个网格；计算域为长方形，计算域面积(X米×Y米)和网格剖分个数的选择需满足

网格的面积Δy·Δx≤25×10^-12(米)²；

步骤二、确定最大形核密度NU_max和最大形核过冷度ΔT_max；

步骤二(一)、一个凝固时间t下，基于高斯连续形核分布原理，采用双曲余弦函数计算微观尺度计算域内形核质点个数MU(t)：

$\mathrm{\ΔT}\left(t\right)=T_l-\stackrel{\·}{T}\×t$

$\mathrm{MU}\left(t\right)=\frac{{\mathrm{NU}}_{\max }}{2}[\tanh \left(\frac{\mathrm{\ΔT}\left(t\right)-{\mathrm{\ΔT}}_{\max }}{1.25}\right)+\tanh \left(\frac{{\mathrm{\ΔT}}_{\max }}{1.25}\right)]\·X\·Y$

其中ΔT(t)为过冷度，T_l为液相线温度，

为冷却速率；

在n×m个网格内随机选取MU(t)个过冷度ΔT达到最大形核过冷度ΔT_max的网格(ΔT≥ΔT_max)；对这MU(t)个网格进行筛选，去除已经成为核心和已被生长枝晶捕获的网格MUsub(t)，余下(MU(t)-MUsub(t))个网格为形核质点存在网格，该类网格被赋予随机结晶取向θ_o(0°≤θ_o≤90°)；

步骤二(二)、一个凝固时间t下，基于尖锐界面模型，针对形核质点存在网格，计算枝晶尖端生长速率V_tip(t)和固相分数f_s(t)：

$\mathrm{\θ}=\arctan (\frac{{\∂f}_s(t-\mathrm{\Δt})}{\∂y}/\frac{{\∂f}_s(t-\mathrm{\Δt})}{\∂x})$

${\mathrm{\μ}}_k\left(\mathrm{\θ}\right)=\stackrel{\‾}{{\mathrm{\μ}}_k}\{1+{\mathrm{\ξ}}_k\cos [4(\mathrm{\θ}-{\mathrm{\θ}}_0)\left]\right\}$

V_tip(t)＝μ_k(θ)·ΔT(t)

$\frac{{\∂f}_s\left(t\right)}{\∂t}=\frac{V_{\mathrm{tip}}\left(t\right)}{\sqrt{\mathrm{\Δy}\·\mathrm{\Δx}}}$

其中t-Δt表示上一时刻，固相分数的初始值

为平均界面动力学系数，ξ_k为动力学各向异性系数，θ固/液界面法线与水平方向的夹角；

基于等式 $f_s\left(t\right)=[\mathrm{MU}\left(t\right)-\mathrm{MUsub}\left(t\right)]\·\mathrm{\π}\·{\left(\frac{\mathrm{SDAS}\left(t\right)}{2}\right)}^2/(X\·Y),$ 计算平均枝晶固相特征直径SDAS(t)；将一个凝固时间下的f_s(t)和SDAS(t)值输出到Nuclei-i.dat文件中，其中固相分数写入第一列，平均枝晶固相特征直径写入第二列，i代表不同的最大形核密度；

步骤二(三)、重复步骤二(一)和步骤二(二)，直到凝固结束；所得数据文件Nuclei-i.dat记录了某一最大形核密度NU_max下，凝固过程中固相分数f_s-平均枝晶固相特征直径SDAS曲线；

步骤三、改变最大形核密度NU_max，重复重复步骤二(一)、步骤二(二)和步骤二(三)，得到不同最大形核密度NU_max下的数据文件Nuclei-i.dat；

步骤四、对铸件(AX米×AY米)进行宏观尺度网格剖分，网格剖分个数(An×Am)的选择需满足

步骤五、确定铸件的最大形核密度NU_max；读取数据文件Nuclei-i.dat，将读取文件得到的固相分数存放在fs_normal数组，将平均枝晶固相特征直径存放在sdas_normal数组；fs_normal[j]对应sdas_normal[j]，fs_normal[j+1]对应sdas_normal[j+1]，j表示某一凝固时间t，j+1表示经过了一个时间步长后的凝固时间t+Δt；

步骤六、计算能量、成分、动量和质量守恒方程，获得铸件内速度、温度T和平均成分的分布；

计算能量守恒方程，求得铸件内温度T分布：

h_s＝c_pT

h_l＝c_pT+ΔH

[H]＝F_sh_s+(1-F_s)h_l

$\frac{\∂\left[\mathrm{\ρH}\right]}{\∂t}+\▿\·({\mathrm{\ρc}}_p\stackrel{\→}{U}T+\mathrm{\ρ}\stackrel{\→}{U}\mathrm{\ΔH})=\▿\·(\mathrm{\λ}\▿T)$

其中h_s和h_l分别为固相和液相热焓，c_p为比热，[H]为混合热焓，ρ为密度，λ为导热系数，ΔH为潜热，T为温度，

为液体流动速度在X，Y方向上的矢量和，

的初始值为0m/s；

计算成分守恒方程，求得获得铸件内平均成分[C]分布：

$\frac{\∂\left[C\right]}{\∂t}+\▿\·\left(\stackrel{\→}{U}{C}_l\right)=0;$

[C]＝F_skC_l+F_lC_l；

T＝T_M+m_lC_l；

其中[C]为平均成分，C_l为液相成分，k为溶质平衡分配系数，T_M为熔点，m_l为液相线m斜率，

的初始值为0m/s；

计算动量守恒方程，求得铸件内流场分布： $K\left(t\right)=\frac{\mathrm{FSDA}{S}^2\×{(1-F_s)}^3}{180\×F_s^2}$

X方向： $\frac{\∂}{\∂t}\left(\mathrm{\ρV}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}V\right)=-\frac{\∂P}{\∂x}-\frac{{\mathrm{\μ}}_l}{K}V+\▿\·({\mathrm{\μ}}_l\▿V)$

Y方向： $\frac{\∂}{\∂t}\left(\mathrm{\ρW}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}W\right)=-\frac{\∂P}{\∂y}-\frac{{\mathrm{\μ}}_l}{K}W+\▿\·({\mathrm{\μ}}_l\▿W)+\mathrm{\ρg}[{\mathrm{\β}}_T(T-T_{\mathrm{ref}})+{\mathrm{\β}}_C(C_l-C_{\mathrm{ref}})]$

其中V、W分别为

在X，Y坐标轴上的分量，P为压力，μ_l为液相粘度，β_T、β_C分别为温度和溶质膨胀系数，T_ref、C_ref分别为金属材料参考温度和参考成分，K为渗透率，FSDAS为平均枝晶固相特征直径，g为重力加速度，K的初始值K(0s)为∞；

所以动量方程的初始时刻计算表达式为：

$\frac{\∂}{\∂t}\left(\mathrm{\ρV}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}V\right)=-\frac{\∂P}{\∂x}+\▿\·({\mathrm{\μ}}_l\▿V)$

$\frac{\∂}{\∂t}\left(\mathrm{\ρW}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}W\right)=-\frac{\∂P}{\∂y}+\▿\·({\mathrm{\μ}}_l\▿W)+\mathrm{\ρg}[{\mathrm{\β}}_T(T-T_{\mathrm{ref}})+{\mathrm{\β}}_C(C_l-C_{\mathrm{ref}})]$

计算质量守恒方程，检验动量方程求解是否准确：

$\▿\·\left(\stackrel{\→}{U}\right)=0$

其中

液体流动速度在X，Y方向上的矢量和；

步骤七、用牛顿下山法计算固相分数F_s，获得铸件内固相分数的分布；

$F_s\left(t\right)=F_s(t-\mathrm{\Δt})-\frac{G_1}{G_{11}}$

$G_1={\mathrm{\ρc}}_p{T}_M-\left[H\left(t\right)\right]+\mathrm{\ρL}(1-F_s(t-\mathrm{\Δt}))-\frac{m_l(1-F_s(t-\mathrm{\Δt})){\mathrm{\ρc}}_p{C}_l\left(t\right)}{(k-1)\·F_s(t-\mathrm{\Δt})+(1-k\·F_s(t-\mathrm{\Δt}))}$

$G_{11}=-\mathrm{\ρL}+\frac{(k-1)m_l(1-F_s(t-\mathrm{\Δt})){\mathrm{\ρc}}_p{C}_l\left(t\right)}{{[(k-1)\·F_s(t-\mathrm{\Δt})+(1-k\·F_s(t-\mathrm{\Δt}))]}^2},$ 其中F_s(t-Δt)为上一时刻的固相分数，固相分数的初始值F_s(0s)＝0；C_l(t)为当前时刻下的液相成分；[H(t)]为当前时刻下的热焓；

步骤八、一个凝固时间t下，采用线性插值技术

即 $\mathrm{FSDAS}\left(t\right)-\mathrm{sdas}\_\mathrm{normal}\left[j\right]=\left(\frac{\mathrm{sdas}\_\mathrm{normal}[j+1]-\mathrm{sdas}\_\mathrm{normal}\left[j\right]}{\mathrm{fs}\_\mathrm{normal}[j+1]-\mathrm{fs}\_\mathrm{normal}\left[j\right]}\right)\×(F_s\left(t\right)-\mathrm{fs}\_\mathrm{normal}[j\left]\right)$ 求解F_s(t)对应的FSDAS(t)，再根据糊状区渗透率模型

获得K(t)；

步骤九、重复步骤六、步骤七和步骤八，直到凝固结束，输出温度、平均成分、液体流动速度和固相分数场。

本发明具有以下有益效果：1.本发明设计了铸件中宏观偏析模拟的方法，对宏观偏析形成进行快速预测，解决了目前宏观偏析计算时间长、精度低的问题。本发明将微观组织模拟与宏观偏析计算进行弱耦合连接，不需要两相流模型那样进行强耦合计算，该方法的速度明显快于两相流模型求解方法的速度，适合大尺寸和形状复杂铸件宏观偏析形成，该方法的精度明显高于经典混合理论。

2.本发明计算过程中微观尺度模拟所输出的不同最大形核密度下的固相分数-平均枝晶固相特征直径曲线，作为输入参数，在各类尺寸和复杂形状的砂型和金属型宏观偏析模拟中直接使用，大大节约了计算时间。

3.本发明适用于各类尺寸和复杂形状的砂型和金属型中宏观偏析的预测。利用本发明可以缩短模拟时间，铸件废品率低，很容易得到生产者的认可，由于市场潜力巨大，一旦被广泛采用，将有几十亿元以上的产值。

附图说明

图1是具体实施例中铸件不同冷却速度下枝晶生长形貌；图2是具体实施例中铸件不同最大形核密度下枝晶生长形貌；图3是具体实施例中铸件不同冷却速度下固相分数-平均枝晶固相特征直径曲线；图4是具体实施例中铸件不同最大形核密度下固相分数-平均枝晶固相特征直径曲线；图5具体实施例中铸件选定的模拟铸件；图6是具体实施例中铸件固相分数-平均枝晶固相特征直径曲线，计算所得糊状区渗透率；图7是具体实施例中铸件基于经典混合理论，计算所得糊状区渗透率；图8是具体实施例中铸件固相分数-平均枝晶固相特征直径曲线模拟所得铸件宏观偏析分布；图9是具体实施例中铸件基于经典混合理论模拟所得铸件宏观偏析分布；图10是具体实施例中铸件宏观偏析分布图，其中——表示铸件基于最大形核密度10⁹/m³对应的固相分数-平均枝晶固相特征直径曲线，模拟所得铸件轴线上宏观偏析分布(提取图8中轴线上宏观偏析分布)，

表示铸件基于经典混合理论，模拟所得铸件轴线上宏观偏析分布(提取图9中轴线上宏观偏析分布)；图11是具体实施例中所得铸件轴线上宏观偏析分布。

具体实施方式

具体实施方式一：本实施方式的铸件宏观偏析数值模拟的方法，所述方法包括如下步骤：

步骤一、将枝晶生长计算域进行微观尺度网格剖分，剖分为n×m个网格；计算域为长方形，计算域面积(X米×Y米)和网格剖分个数的选择需满足

网格的面积Δy·Δx≤25×10^-12(米)²；

步骤二、确定最大形核密度NU_max和最大形核过冷度ΔT_max；

步骤二(一)、一个凝固时间t下，基于高斯连续形核分布原理，采用双曲余弦函数计算微观尺度计算域内形核质点个数MU(t)：

$\mathrm{\ΔT}\left(t\right)=T_l-\stackrel{\·}{T}\×t$

$\mathrm{MU}\left(t\right)=\frac{{\mathrm{NU}}_{\max }}{2}[\tanh \left(\frac{\mathrm{\ΔT}\left(t\right)-{\mathrm{\ΔT}}_{\max }}{1.25}\right)+\tanh \left(\frac{{\mathrm{\ΔT}}_{\max }}{1.25}\right)]\·X\·Y$

其中ΔT(t)为过冷度，T_l为液相线温度，

为冷却速率；

在n×m个网格内随机选取MU(t)个过冷度ΔT达到最大形核过冷度ΔT_max的网格(ΔT≥ΔT_max)；对这MU(t)个网格进行筛选，去除已经成为核心和已被生长枝晶捕获的网格MUsub(t)，余下(MU(t)-MUsub(t))个网格为形核质点存在网格，该类网格被赋予随机结晶取向θ_o(0°≤θ_o≤90°)；

步骤二(二)、一个凝固时间t下，基于尖锐界面模型，针对形核质点存在网格，计算枝晶尖端生长速率V_tip(t)和固相分数f_s(t)：

$\mathrm{\θ}=\arctan (\frac{{\∂f}_s(t-\mathrm{\Δt})}{\∂y}/\frac{{\∂f}_s(t-\mathrm{\Δt})}{\∂x})$

${\mathrm{\μ}}_k\left(\mathrm{\θ}\right)=\stackrel{\‾}{{\mathrm{\μ}}_k}\{1+{\mathrm{\ξ}}_k\cos [4(\mathrm{\θ}-{\mathrm{\θ}}_0)\left]\right\}$

V_tip(t)＝μ_k(θ)·ΔT(t)

$\frac{{\∂f}_s\left(t\right)}{\∂t}=\frac{V_{\mathrm{tip}}\left(t\right)}{\sqrt{\mathrm{\Δy}\·\mathrm{\Δx}}}$

其中t-Δt表示上一时刻，固相分数的初始值f_s(0s)＝0，为平均界面动力学系数，ξ_k为动力学各向异性系数，θ固/液界面法线与水平方向的夹角；

基于等式 $f_s\left(t\right)=[\mathrm{MU}\left(t\right)-\mathrm{MUsub}\left(t\right)]\·\mathrm{\π}\·{\left(\frac{\mathrm{SDAS}\left(t\right)}{2}\right)}^2/(X\·Y),$ 计算平均枝晶固相特征直径SDAS(t)；将一个凝固时间下的f_s(t)和SDAS(t)值输出到Nuclei-i.dat文件中，其中固相分数写入第一列，平均枝晶固相特征直径写入第二列，i代表不同的最大形核密度；

步骤二(三)、重复步骤二(一)和步骤二(二)，直到凝固结束；所得数据文件Nuclei-i.dat记录了某一最大形核密度NU_max下，凝固过程中固相分数f_s-平均枝晶固相特征直径SDAS曲线；

步骤三、改变最大形核密度NU_max，重复重复步骤二(一)、步骤二(二)和步骤二(三)，得到不同最大形核密度NU_max下的数据文件Nuclei-i.dat；

步骤四、对铸件(AX米×AY米)进行宏观尺度网格剖分，网格剖分个数(An×Am)的选择需满足

步骤五、确定铸件的最大形核密度NU_max；读取数据文件Nuclei-i.dat，将读取文件得到的固相分数存放在fs_normal数组，将平均枝晶固相特征直径存放在sdas_normal数组；fs_normal[j]对应sdas_normal[j]，fs_normal[j+1]对应sdas_normal[j+1]，j表示某一凝固时间t，j+1表示经过了一个时间步长后的凝固时间t+Δt；

步骤六、计算能量、成分、动量和质量守恒方程，获得铸件内速度、温度T和平均成分的分布；

计算能量守恒方程，求得铸件内温度T分布：

h_s＝c_pT

h_l＝c_pT+ΔH

[H]＝F_sh_s+(1-F_s)h_l

$\frac{\∂\left[\mathrm{\ρH}\right]}{\∂t}+\▿\·({\mathrm{\ρc}}_p\stackrel{\→}{U}T+\mathrm{\ρ}\stackrel{\→}{U}\mathrm{\ΔH})=\▿\·(\mathrm{\λ}\▿T)$

其中h_s和h_l分别为固相和液相热焓，c_p为比热，[H]为混合热焓，ρ为密度，λ为导热系数，ΔH为潜热，T为温度，

为液体流动速度在X，Y方向上的矢量和，

的初始值为0m/s；

计算成分守恒方程，求得获得铸件内平均成分[C]分布：

$\frac{\∂\left[C\right]}{\∂t}+\▿\·\left(\stackrel{\→}{U}{C}_l\right)=0;$

[C]＝F_skC_l+F_lC_l；

T＝T_M+m_lC_l；

其中[C]为平均成分，C_l为液相成分，k为溶质平衡分配系数，T_M为熔点，m_l为液相线斜率，

的初始值为0m/s；

计算动量守恒方程，求得铸件内流场分布： $K\left(t\right)=\frac{\mathrm{FSDA}{S}^2\×{(1-F_s)}^3}{180\×{F_s}^2}$

X方向： $\frac{\∂}{\∂t}\left(\mathrm{\ρV}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}V\right)=-\frac{\∂P}{\∂x}-\frac{{\mathrm{\μ}}_l}{K}V+\▿\·({\mathrm{\μ}}_l\▿V)$

Y方向： $\frac{\∂}{\∂t}\left(\mathrm{\ρW}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}W\right)=-\frac{\∂P}{\∂y}-\frac{{\mathrm{\μ}}_l}{K}W+\▿\·({\mathrm{\μ}}_l\▿W)+\mathrm{\ρg}[{\mathrm{\β}}_T(T-T_{\mathrm{ref}})+{\mathrm{\β}}_C(C_l-C_{\mathrm{ref}})]$

其中V、W分别为

在X，Y坐标轴上的分量，P为压力，μ_l为液相粘度，β_T、β_C分别为温度和溶质膨胀系数，T_ref、C_ref分别为金属材料参考温度和参考成分，K为渗透率，FSDAS为平均枝晶固相特征直径，g为重力加速度，K的初始值K(0s)为∞；

所以动量方程的初始时刻计算表达式为：

$\frac{\∂}{\∂t}\left(\mathrm{\ρV}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}V\right)=-\frac{\∂P}{\∂x}+\▿\·({\mathrm{\μ}}_l\▿V)$

$\frac{\∂}{\∂t}\left(\mathrm{\ρW}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}W\right)=-\frac{\∂P}{\∂y}+\▿\·({\mathrm{\μ}}_l\▿W)+\mathrm{\ρg}[{\mathrm{\β}}_T(T-T_{\mathrm{ref}})+{\mathrm{\β}}_C(C_l-C_{\mathrm{ref}})]$

计算质量守恒方程，检验动量方程求解是否准确：

$\▿\·\left(\stackrel{\→}{U}\right)=0$

其中

液体流动速度在X，Y方向上的矢量和；

步骤七、用牛顿下山法计算固相分数F_s，获得铸件内固相分数的分布；

$F_s\left(t\right)=F_s(t-\mathrm{\Δt})-\frac{G_1}{G_{11}}$

$G_1={\mathrm{\ρc}}_p{T}_M-\left[H\left(t\right)\right]+\mathrm{\ρL}(1-F_s(t-\mathrm{\Δt}))-\frac{m_l(1-F_s(t-\mathrm{\Δt})){\mathrm{\ρc}}_p{C}_l\left(t\right)}{(k-1)\·F_s(t-\mathrm{\Δt})+(1-k\·F_s(t-\mathrm{\Δt}))}$

$G_{11}=-\mathrm{\ρL}+\frac{(k-1)m_l(1-F_s(t-\mathrm{\Δt})){\mathrm{\ρc}}_p{C}_l\left(t\right)}{{[(k-1)\·F_s(t-\mathrm{\Δt})+(1-k\·F_s(t-\mathrm{\Δt}))]}^2},$ 其中F_s(t-Δt)为上一时刻的固相分数，固相分数的初始值F_s(0s)＝0；C_l(t)为当前时刻下的液相成分；[H(t)]为当前时刻下的热焓；

步骤八、一个凝固时间t下，采用线性插值技术

即

$\mathrm{FSDAS}\left(t\right)-\mathrm{sdas}\_\mathrm{normal}\left[j\right]=\left(\frac{\mathrm{sdas}\_\mathrm{normal}[j+1]-\mathrm{sdas}\_\mathrm{normal}\left[j\right]}{\mathrm{fs}\_\mathrm{normal}[j+1]-\mathrm{fs}\_\mathrm{normal}\left[j\right]}\right)\×(F_s\left(t\right)-\mathrm{fs}\_\mathrm{normal}[j\left]\right)$

求解F_s(t)对应的FSDAS(t)，再根据糊状区渗透率模型

获得K(t)；

步骤九、重复步骤六、步骤七和步骤八，直到凝固结束，输出温度、平均成分、液体流动速度和固相分数场。

具体实施方式二：本实施方式的步骤八中所述糊状区渗透率模型的计算采用两种方法，一种计算方法为基于最大形核密度10⁹/m³对应的固相分数-平均枝晶固相特征直径曲线，一种计算方法为基于经典混合理论，计算得出的结果选取其中一组与实验最为接近的数据。其他组成及连接关系与具体实施方式一相同。

具体实施方式三：本实施方式的步骤四中设定(Bn×Bm)个网格为柱状晶生长网格，(An×Am-Bn×Bm)个网格为等轴晶生长网格；步骤五中确定铸件的读取所对应的数据文件Nuclei.dat；步骤九中柱状晶最大形核密度LNU_max和等轴晶最大形核密度ENU_max，状晶生长区糊状区渗透率计算基于最大形核密度LNU_max对应的固相分数-平均枝晶固相特征直径曲线，等轴晶生长区糊状区渗透率计算基于最大形核密度ENU_max对应的固相分数-平均枝晶固相特征直径曲线，此步骤为了模拟凝固过程中当微观组织发生演变时，即组织由柱状晶变为等轴晶时，宏观偏析的形成。其他组成及连接关系与具体实施方式一相同。

具体实施例：

步骤一：以Fe-0.22wt％C合金为研究对象，计算域的面积选择0.25cm×0.25cm，进行微观尺度网格剖分，剖分500×500个网格，保证

步骤二：采用双曲余弦函数计算形核个数；进行晶核的随机分布；赋予随机结晶取向θ_o；采用元胞自动机机模型、尖锐界面模型，针对形核质点存在网格，计算枝晶尖端生长速率V_tip和固相分数f_s(t)；计算所需输入参数列于表1；

表1

步骤三：基于等式 $f_s\left(t\right)=[\mathrm{MU}\left(t\right)-\mathrm{MUsub}\left(t\right)]\·\mathrm{\π}\·{\left(\frac{\mathrm{SDAS}\left(t\right)}{2}\right)}^2/(X\·Y),$ 计算平均枝晶固相特征直径SDAS(t)；将一个凝固时间下的f_s(t)和SDAS(t)值输出到Nuclei-i.dat文件中，其中固相分数写入第一列，平均枝晶固相特征直径写入第二列，i代表不同的最大形核密度；

步骤四：改变冷却速率和最大形核密度，重复步骤(2)、(3)，模拟合金凝固过程中枝晶生长形貌，结果示于图1和图2；

步骤五：将输出的固相分数-平均枝晶固相特征直径曲线(Nuclei-i.dat)，用绘图软件进行显示，结果示于图3和图4；

步骤六：模拟铸件如图5所示，铸件被铸型包围，铸件上部有覆盖剂和保温砖，铸件底部有钢板和保温板。宏观尺度网格剖分尺寸为2cm，满足

铸件材质为Fe-0.22wt％C合金，计算输入参数列于表2；

表2

步骤七：计算能量、成分、动量和质量守恒万程，获得铸件内速度、温度T和平均成分的分布；

计算能量守恒方程，求得铸件内温度T分布：

h_s＝c_pT

h_l＝c_pT+ΔH

[H]＝F_sh_s+(1-F_s)h_l

$\frac{\∂\left[\mathrm{\ρH}\right]}{\∂t}+\▿\·({\mathrm{\ρc}}_p\stackrel{\→}{U}T+\mathrm{\ρ}\stackrel{\→}{U}\mathrm{\ΔH})=\▿\·\left(\mathrm{\λ\ΔT}\right)$

其中h_s和h_l分别为固相和液相热焓，c_p为比热，[H]为混合热焓，ρ为密度，λ为导热系数，ΔH为潜热，T为温度，U为液体流动速度在X，Y方向上的矢量和，

的初始值为0m/s；

计算成分守恒方程，求得获得铸件内平均成分[C]分布：

$\frac{\∂\left[C\right]}{\∂t}+\▿\·\left(\stackrel{\→}{U}{C}_l\right)=0;$

[C]＝F_skC_l+F_lC_l；

T＝T_M+m_lC_l；

其中[C]为平均成分，C_l为液相成分，k为溶质平衡分配系数，T_M为熔点，m_l为液相线斜率，的初始值为0m/s；

计算动量守恒方程，求得铸件内流场分布： $K\left(t\right)=\frac{\mathrm{FSDA}{S}^2\×{(1-F_s)}^3}{180\×F_s^2}$

X方向： $\frac{\∂}{\∂t}\left(\mathrm{\ρV}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}V\right)=-\frac{\∂P}{\∂x}-\frac{{\mathrm{\μ}}_l}{K}V+\▿\·({\mathrm{\μ}}_l\▿V)$

Y方向： $\frac{\∂}{\∂t}\left(\mathrm{\ρW}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}W\right)=-\frac{\∂P}{\∂y}-\frac{{\mathrm{\μ}}_l}{K}W+\▿\·({\mathrm{\μ}}_l\▿W)+\mathrm{\ρg}[{\mathrm{\β}}_T(T-T_{\mathrm{ref}})+{\mathrm{\β}}_C(C_l-C_{\mathrm{ref}})]$

其中V、W分别为

在X，Y坐标轴上的分量，P为压力，μ_l为液相粘度，β_T、β_C分别为温度和溶质膨胀系数，T_ref、C_ref分别为金属材料参考温度和参考成分，K为渗透率，FSDAS为平均枝晶固相特征直径，g为重力加速度，K的初始值K(0s)为∞；

所以动量方程的初始时刻计算表达式为：

$\frac{\∂}{\∂t}\left(\mathrm{\ρV}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}V\right)=-\frac{\∂P}{\∂x}+\▿\·({\mathrm{\μ}}_l\▿V)$

$\frac{\∂}{\∂t}\left(\mathrm{\ρW}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}W\right)=-\frac{\∂P}{\∂y}+\▿\·({\mathrm{\μ}}_l\▿W)+\mathrm{\ρg}[{\mathrm{\β}}_T(T-T_{\mathrm{ref}})+{\mathrm{\β}}_C(C_l-C_{\mathrm{ref}})]$

计算质量守恒方程，检验动量方程求解是否准确：

$\▿\·\left(\stackrel{\→}{U}\right)=0$

其中

液体流动速度在X，Y方向上的矢量和；计算所得的平均成分[C]示于图8和图9；

步骤八、用牛顿下山法计算固相分数F_s，获得铸件内固相分数的分布；

$F_s\left(t\right)=F_s(t-\mathrm{\Δt})-\frac{G_1}{G_{11}}$

$G_1={\mathrm{\ρc}}_p{T}_M-\left[H\left(t\right)\right]+\mathrm{\ρL}(1-F_s(t-\mathrm{\Δt}))-\frac{m_l(1-F_s(t-\mathrm{\Δt})){\mathrm{\ρc}}_p{C}_l\left(t\right)}{(k-1)\·F_s(t-\mathrm{\Δt})+(1-k\·F_s(t-\mathrm{\Δt}))}$

$G_{11}=-\mathrm{\ρL}+\frac{(k-1)m_l(1-F_s(t-\mathrm{\Δt})){\mathrm{\ρc}}_p{C}_l\left(t\right)}{{[(k-1)\·F_s(t-\mathrm{\Δt})+(1-k\·F_s(t-\mathrm{\Δt}))]}^2},$ 其中F_s(t-Δt)为上一时刻的固相分数，固相分数的初始值F_s(0s)＝0；C_l(t)为当前时刻下的液相成分；[H(t)]为当前时刻下的热焓；

步骤九：宏观尺度糊状区渗透率计算采用两种方法，一种为基于最大形核密度10⁹/m³对应的固相分数-平均枝晶固相特征直径曲线，一种为基于经典混合理论，结果示于图6和图7，对比可知图6与实际更为接近，图7中对凝固初期糊状区渗透率的计算偏大；

步骤十：模拟所得铸件上宏观偏析分布示于图8和图9。基于固相分数-平均枝晶固相特征直径曲线计算所得宏观偏析分布程度较小，因为凝固初期糊状区渗透率的计算较小，液体流动速度减弱；

步骤十一：模拟所得铸件轴线上宏观偏析分布示于图10，实验所得铸件轴线上宏观偏析分布示于图11。基于固相分数-平均枝晶固相特征直径曲线计算所得宏观偏析分布与实验吻合良好；

其中当需要在宏观偏析预测过程中考虑柱状晶向等轴晶转变时，在具体实施例所述的一种铸件宏观偏析数值模拟的方法基础上，调整步骤九：将铸件划分为两个区域：柱状晶生长区和等轴晶生长区。柱状晶生长区糊状区渗透率计算基于最大形核密度10⁹/m³对应的固相分数-平均枝晶固相特征直径曲线，等轴晶生长区糊状区渗透率计算基于最大形核密度10⁸/m³对应的固相分数-平均枝晶固相特征直径曲线。

Claims (3)

1.一种铸件宏观偏析数值模拟的方法，其特征在于：所述方法包括如下步骤：

步骤一、将枝晶生长计算域进行微观尺度网格剖分，剖分为n×m个网格；计算域为长方形，计算域面积(X米×Y米)和网格剖分个数的选择需满足

网格的面积Δy·Δx≤25×10^-12(米)²；

步骤二、确定最大形核密度NU_max和最大形核过冷度ΔT_max；

步骤二(一)、一个凝固时间t下，基于高斯连续形核分布原理，采用双曲余弦函数计算微观尺度计算域内形核质点个数MU(t)：

$\mathrm{\ΔT}\left(t\right)=T_l-\stackrel{\·}{T}\×t$

$\mathrm{MU}\left(t\right)=\frac{{\mathrm{NU}}_{\max }}{2}[\tanh \left(\frac{\mathrm{\ΔT}\left(t\right)-{\mathrm{\ΔT}}_{\max }}{1.25}\right)+\tanh \left(\frac{{\mathrm{\ΔT}}_{\max }}{1.25}\right)]\·X\·Y$

其中ΔT(t)为过冷度，T_l为液相线温度，

为冷却速率；

在n×m个网格内随机选取MU(t)个过冷度ΔT达到最大形核过冷度ΔT_max的网格(ΔT≥ΔT_max)；对这MU(t)个网格进行筛选，去除已经成为核心和已被生长枝晶捕获的网格MUsub(t)，余下(MU(t)-MUsub(t))个网格为形核质点存在网格，该类网格被赋予随机结晶取向θ_o(0°≤θ_o≤90°)；

步骤二(二)、一个凝固时间t下，基于尖锐界面模型，针对形核质点存在网格，计算枝晶尖端生长速率V_tip(t)和固相分数f_s(t)：

$\mathrm{\θ}=\arctan (\frac{{\∂f}_s(t-\mathrm{\Δt})}{\∂y}/\frac{{\∂f}_s(t-\mathrm{\Δt})}{\∂x})$

${\mathrm{\μ}}_k\left(\mathrm{\θ}\right)=\stackrel{\‾}{{\mathrm{\μ}}_k}\{1+{\mathrm{\ξ}}_k\cos [4(\mathrm{\θ}-{\mathrm{\θ}}_0)\left]\right\}$

V_tip(t)＝μ_k(θ)·ΔT(t)

$\frac{{\∂f}_s\left(t\right)}{\∂t}=\frac{V_{\mathrm{tip}}\left(t\right)}{\sqrt{\mathrm{\Δy}\·\mathrm{\Δx}}}$

其中t-Δt表示上一时刻，固相分数的初始值f_s(0s)＝0，

为平均界面动力学系数，ξ_k为动力学各向异性系数，θ固/液界面法线与水平方向的夹角；

基于等式 $f_s\left(t\right)=[\mathrm{MU}\left(t\right)-\mathrm{MUsub}\left(t\right)]\·\mathrm{\π}\·{\left(\frac{\mathrm{SDAS}\left(t\right)}{2}\right)}^2/(X\·Y),$ 计算平均枝晶固相特征直径SDAS(t)；将一个凝固时间下的f_s(t)和SDAS(t)值输出到Nuclei-i.dat文件中，其中固相分数写入第一列，平均枝晶固相特征直径写入第二列，i代表不同的最大形核密度；

步骤二(三)、重复步骤二(一)和步骤二(二)，直到凝固结束；所得数据文件Nuclei-i.dat记录了某一最大形核密度NU_max下，凝固过程中固相分数f_s-平均枝晶固相特征直径SDAS曲线；

步骤三、改变最大形核密度NU_max，重复重复步骤二(一)、步骤二(二)和步骤二(三)，得到不同最大形核密度NU_max下的数据文件Nuclei-i.dat；

步骤四、对铸件(AX米×AY米)进行宏观尺度网格剖分，网格剖分个数(An×Am)的选择需满足

步骤五、确定铸件的最大形核密度NU_max；读取数据文件Nuclei-i.dat，将读取文件得到的固相分数存放在fs_normal数组，将平均枝晶固相特征直径存放在sdas_normal数组；fs_normal[j]对应sdas_normal[j]，fs_normal[j+1]对应sdas_normal[j+1]，j表示某一凝固时间t，j+1表示经过了一个时间步长后的凝固时间t+Δt；

步骤六、计算能量、成分、动量和质量守恒方程，获得铸件内速度、温度T和平均成分的分布；

计算能量守恒方程，求得铸件内温度T分布：

h_s＝c_pT

h_l＝c_pT+ΔH

[H]＝F_sh_s+(1-F_s)h_l

$\frac{\∂\left[\mathrm{\ρH}\right]}{\∂t}+\▿\·({\mathrm{\ρc}}_p\stackrel{\→}{U}T+\mathrm{\ρ}\stackrel{\→}{U}\mathrm{\ΔH})=\▿\·(\mathrm{\λ}\▿T)$

其中h_s和h_l分别为固相和液相热焓，c_p为比热，[H]为混合热焓，ρ为密度，λ为导热系数，ΔH为潜热，T为温度，

为液体流动速度在X，Y方向上的矢量和，

的初始值为0m/s；

计算成分守恒方程，求得获得铸件内平均成分[C]分布：

$\frac{\∂\left[C\right]}{\∂t}+\▿\·\left(\stackrel{\→}{U}{C}_l\right)=0;$

[C]＝F_skC_l+F_lC_l；

T＝T_M+m_lC_l；

其中[C]为平均成分，C_l为液相成分，k为溶质平衡分配系数，T_M为熔点，m_l为液相线斜率，

的初始值为0m/s；

计算动量守恒方程，求得铸件内流场分布： $K\left(t\right)=\frac{{\mathrm{FSDAS}}^2\×{(1-F_s)}^3}{180\×F_s^2}$

X方向： $\frac{\∂}{\∂t}\left(\mathrm{\ρV}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}V\right)=-\frac{\∂P}{\∂x}-\frac{{\mathrm{\μ}}_l}{K}V+\▿\·({\mathrm{\μ}}_l\▿V)$

Y方向： $\frac{\∂}{\∂t}\left(\mathrm{\ρW}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}W\right)=-\frac{\∂P}{\∂y}-\frac{{\mathrm{\μ}}_l}{K}W+\▿\·({\mathrm{\μ}}_l\▿W)+\mathrm{\ρg}[{\mathrm{\β}}_T(T-T_{\mathrm{ref}})+{\mathrm{\β}}_C(C_l-C_{\mathrm{ref}})]$

其中V、W分别为

在X，Y坐标轴上的分量，P为压力，μ_l为液相粘度，β_T、β_C分别为温度和溶质膨胀系数，T_ref、C_ref分别为金属材料参考温度和参考成分，K为渗透率，FSDAS为平均枝晶固相特征直径，g为重力加速度，K的初始值K(0s)为∞；

所以动量方程的初始时刻计算表达式为：

$\frac{\∂}{\∂t}\left(\mathrm{\ρV}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}V\right)=-\frac{\∂P}{\∂x}+\▿\·({\mathrm{\μ}}_l\▿V)$

$\frac{\∂}{\∂t}\left(\mathrm{\ρW}\right)+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U}W\right)=-\frac{\∂P}{\∂y}+\▿\·({\mathrm{\μ}}_l\▿W)+\mathrm{\ρg}[{\mathrm{\β}}_T(T-T_{\mathrm{ref}})+{\mathrm{\β}}_C(C_l-C_{\mathrm{ref}})]$

计算质量守恒方程，检验动量方程求解是否准确：

$\▿\·\left(\stackrel{\→}{U}\right)=0$

其中

液体流动速度在X，Y方向上的矢量和；

步骤七、用牛顿下山法计算固相分数F_s，获得铸件内固相分数的分布；

$F_s\left(t\right)=F_s(t-\mathrm{\Δt})-\frac{G_1}{G_{11}}$

$G_1={\mathrm{\ρc}}_p{T}_M-\left[H\left(t\right)\right]+\mathrm{\ρL}(1-F_s(t-\mathrm{\Δt}))-\frac{m_l(1-F_s(t-\mathrm{\Δt})){\mathrm{\ρc}}_p{C}_l\left(t\right)}{(k-1)\·F_s(t-\mathrm{\Δt})+(1-k\·F_s(t-\mathrm{\Δt}))}$

$G_{11}=-\mathrm{\ρL}+\frac{(k-1)m_l(1-F_s(t-\mathrm{\Δt})){\mathrm{\ρc}}_p{C}_l\left(t\right)}{{[(k-1)\·F_s(t-\mathrm{\Δt})+(1-k\·F_s(t-\mathrm{\Δt}))]}^2},$ 其中F_s(t-Δt)为上一时刻的固相分数，固相分数的初始值F_s(0s)＝0；C_l(t)为当前时刻下的液相成分；[H(t)]为当前时刻下的热焓；

步骤八、一个凝固时间t下，采用线性插值技术

即

$\mathrm{FSDAS}\left(t\right)-\mathrm{sdas}\_\mathrm{normal}\left[j\right]=\left(\frac{\mathrm{sdas}\_\mathrm{normal}[j+1]-\mathrm{sdas}\_\mathrm{normal}\left[j\right]}{\mathrm{fs}\_\mathrm{normal}[j+1]-\mathrm{fs}\_\mathrm{normal}\left[j\right]}\right)\×(F_s\left(t\right)-\mathrm{fs}\_\mathrm{normal}[j\left]\right)$

求解F_s(t)对应的FSDAS(t)，再根据糊状区渗透率模型获得K(t)；

步骤九、重复步骤六、步骤七和步骤八，直到凝固结束，输出温度、平均成分、液体流动速度和固相分数场。

2.根据权利要求1所述一种铸件宏观偏析数值模拟的方法，其特征在于步骤八中所述糊状区渗透率模型的计算采用两种方法中的一种，一种计算方法为基于最大形核密度10⁹/m³对应的固相分数-平均枝晶固相特征直径曲线，一种计算方法为基于经典混合理论。

3.根据权利要求1或2所述一种铸件宏观偏析数值模拟的方法，其特征在于步骤四中设定(Bn×Bm)个网格为柱状晶生长网格，(An×Am-Bn×Bm)个网格为等轴晶生长网格，步骤五中确定铸件的读取所对应的数据文件Nuclei.dat，步骤九中柱状晶最大形核密度LNU_max和等轴晶最大形核密度ENU_max，状晶生长区糊状区渗透率计算基于最大形核密度LNU_max对应的固相分数-平均枝晶固相特征直径曲线，等轴晶生长区糊状区渗透率计算基于最大形核密度ENU_max对应的固相分数-平均枝晶固相特征直径曲线。

Images