CN109885984A - 一种球墨铸铁铸锭石墨球尺寸数值预测的方法 - Google Patents

Abstract

本发明涉及一种球墨铸铁铸锭石墨球尺寸数值预测的方法，是为了解决现有技术中仅针对中心截面的二维计算由于没有考虑壁面效应，从而影响过冷度计算；在三维方向计算温度场和速度场会增加计算量、大大延长计算时间的缺点而提出的，包括：对铸造系统进行网格划分；对所有网格计算能量守恒方程，获得三维方向铸造系统温度场分布；对于铸锭中平行于重力方向的中心截面网格计算动量守恒方程，获得该计算网格中的金属液流动速度；对于铸锭网格和铸锭中平行于重力方向的中心截面网格，计算石墨球尺寸；重复上述步骤直至所有铸锭网格和铸锭中平行于重力方向的中心截面网格的温度均小于共晶线温度。本发明适用于各类尺寸的砂型和金属型中石墨球尺寸的预测。

Description

一种球墨铸铁铸锭石墨球尺寸数值预测的方法

本发明申请是申请日为2018年06月22日、申请号为201810654222.3、发明名称为一种球墨铸铁铸锭石墨球尺寸数值预测的方法的发明申请的优先权申请。

技术领域

本发明涉及球墨铸铁铸锭仿真领域，具体涉及一种球墨铸铁铸锭石墨球尺寸数值预测的方法。

背景技术

铸铁中球状石墨的形成可以有效提高铸铁的塑性、韧性和强度。含有大量石墨球的铸铁被称为球墨铸铁，基于其优异的力学性能，该材料被成功应用于一些受力复杂或对强度、韧性、耐磨性有较高要求的零件。

球墨铸铁铸锭是一种母材，需经塑性加工得到各种用途。工业上可以用于球墨铸铁管、球墨铸铁井盖、船用气缸套的加工。为了提高球墨铸铁铸锭的综合力学性能，需要对石墨球尺寸进行控制。球化率高(石墨球个数多)、石墨球细小是提高球墨铸铁铸锭力学性能的关键。石墨球产生于铸锭凝固过程，凝固过程中的过冷度是石墨球形核和生长的驱动力。铸锭凝固是一个复杂的物理过程，涉及多个尺度物理现象的相互叠交：宏观尺度动量、热量、质量传输，介观尺度晶粒组织形成，微观尺度溶质扩散，不同尺度的物理现象之间相互作用、相互影响。采用实验方法研究铸锭凝固过程中石墨球尺寸不仅浪费了人力和财力，而且无法获得不同凝固条件下石墨球尺寸的变化特点，无法获得凝固过程中石墨球尺寸随时间变化曲线。更重要的是由于铸锭尺寸较大，每一次浇注成型都会耗费大量能源，因此大量的实验研究必然造成能源浪费，使环境受到破坏和污染。所以利用计算机模拟技术预测球墨铸铁铸锭凝固过程中石墨球尺寸变化，是控制和提升球墨铸铁铸锭产品质量的一个行之有效的手段。

石墨球的形成由形核和长大两阶段组成。过冷度是形核驱动力，因为它反映了液相和固相自由能的相差程度。过冷度越大，自由能相差程度就越大，球状石墨更易形核。过冷度也是生长驱动力，过冷度增大表明固液界面前沿的碳元素成分梯度增大，大的成分梯度下碳元素从液相中向石墨球中的扩散通量增加，从而提高石墨球的生长速率。准确预测铸锭凝固过程温度场从而获得温度随时间变化曲线是预测石墨球尺寸的关键。目前针对铸锭的研究多采用二维计算或三维计算。

其中仅针对中心截面的二维计算由于没有考虑壁面效应，会低估冷速，从而影响过冷度计算；在三维方向计算温度场和速度场会增加计算量、大大延长计算时间。

发明内容

本发明的目的是为了解决现有技术中仅针对中心截面的二维计算由于没有考虑壁面效应，会低估冷速，从而影响过冷度计算；在三维方向计算温度场和速度场会增加计算量、大大延长计算时间缺点，而提出一种球墨铸铁铸锭石墨球尺寸数值预测的方法，包括：

步骤一、对铸造系统进行网格划分；

步骤二、对所有网格计算能量守恒方程，获得三维方向铸造系统温度场分布；

步骤三、对于铸锭中平行于重力方向的中心截面网格计算动量守恒方程，获得该计算网格中的金属液流动速度；

步骤四、对于铸锭网格和铸锭中平行于重力方向的中心截面网格，计算石墨球尺寸；

步骤五、重复步骤二至四，直至所有铸锭网格和铸锭中平行于重力方向的中心截面网格的温度T_in均小于共晶线温度T_E。

本发明具有以下有益效果：本发明的方法采用了连续形核二维石墨生长模型，考虑了三维温度场计算，更加准确地预测了凝固过程中石墨球尺寸随时间变化特点，解决了目前石墨球尺寸数值预测中不考虑三维能量传输、不考虑石墨连续生长特点的问题，为分析球墨铸铁中石墨球形成供了数据参考。

附图说明

图1为本发明一个实施例的铸造系统三维示意图；

图2(a)为本发明一个实施例的铸锭三维示意图；图2(b)为本发明一个实施例的铸锭中心截面图，图中P1和P2是在中心界面上选取的两个相同高度的点；

图3(a)为本发明一个实施例中铸锭温度场二维模拟、液体流动仅针对中心截面二维模拟时的P1点和P2点的冷却曲线；图3(b)为本发明一个实施例中铸锭温度场三维模拟、液体流动仅针对中心截面二维模拟时P1点和P2点的冷却曲线；图3(a)和图3(b)的横轴均为时间，纵轴均为温度；

图4(a)为本发明一个实施例中铸锭温度场二维模拟、液体流动仅针对中心截面二维模拟时4805s时刻铸锭中心截面温度场分布图；图4(b)为本发明一个实施例中铸锭温度场二维模拟、液体流动仅针对中心截面二维模拟时和12005s时刻铸锭中心截面温度场分布图；

图5(a)为本发明一个实施例中铸锭温度场三维模拟、液体流动仅针对中心截面二维模拟时4805s时刻铸锭中心截面温度场分布图；图5(b)为本发明一个实施例中铸锭温度场三维模拟、液体流动仅针对中心截面二维模拟时12005s时刻铸锭中心截面温度场分布图；

图6(a)为铸锭温度场二维模拟、液体流动仅针对中心截面二维模拟时P1点和P2点处石墨球半径随时间变化曲线图；图6(b)为铸锭温度场三维模拟、液体流动仅针对中心截面二维模拟时P1点和P2点处石墨球半径随时间变化曲线图；图6(a)和图6(b)的横轴均为时间，纵轴均为石墨球半径；

图7为本发明一个实施例的球墨铸铁铸锭石墨球尺寸数值预测的方法的流程图。

具体实施方式

具体实施方式一：本实施方式的球墨铸铁铸锭石墨球尺寸数值预测的方法，如图7所示，包括：

步骤一、对铸造系统进行网格划分。其中X轴、Y轴、Z轴可以为相互正交的任意坐标轴。X轴、Y轴、Z轴的选取可以依照实际情况而定。坐标系的选择不同不会影响预测结果。

步骤二、对所有网格计算能量守恒方程，获得三维方向铸造系统温度场分布。

步骤三、对于铸锭中平行于重力方向的中心截面网格计算动量守恒方程，获得该计算网格中的金属液流动速度。

步骤四、对于铸锭网格和铸锭中平行于重力方向的中心截面网格，计算石墨球尺寸。

步骤五、重复步骤二至四，直至所有铸锭网格和铸锭中平行于重力方向的中心截面网格的温度T_in均小于共晶线温度T_E。

本发明所作出的改进是基于最新理论而提出的。最新的理论研究表明，石墨片可以通过二维搭接方式进行长大从而长成球状石墨。相比于基于螺旋位错石墨长大模型，二维搭接石墨长大模型可以再现凝固过程中球状石墨连续生长状态。因此本发明的球墨铸铁铸锭石墨球尺寸数值预测的方法，采用更符合物理实际的连续形核的二维搭接石墨长大模型，同时采用三维模型计算能量传输，二维模型计算动量传输，研究冷却速度与石墨球尺寸之间的关系，这在明晰不同铸造工艺下球墨铸铁铸锭石墨球尺寸变化特点方面具有重要意义。

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

对X米×Y米×Z米的铸造系统进行宏观尺度网格剖分，X方向、Y方向和Z方向分别采用Δx、Δy、Δz作为网格剖分步长，Δx米＝Δy米＝Δz米，Δx、Δy和Δz的取值范围1×10^-3米～4×10^-3米，计算网格的标号为(i，j，k)_char，其中i，j和k均为整数，i的取值范围是1～L，j的取值范围是1～M，k的取值范围是1～N， 下角标char＝2表示铸型网格，下角标char＝0表示铸锭网格，下角标char＝21表示铸锭中平行于重力方向的中心截面网格，下角标char＝4、5、6和7分别表示内冷铁网格、冒口套网格、保温材料网格和绝热材料网格；铸造系统在X轴、Y轴、Z轴方向上的最小值分别为X_min、Y_min、Z_min(米)，在X轴、Y轴、Z轴方向上的最大值分别为X_max、Y_max、Z_max(米)。

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

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤二中，针对所有计算网格(i，j，k)_char只计算能量守恒方程，即计算温度场，获得三维方向铸造系统温度场分布。

步骤二具体为：

若计算网格(i，j，k)_char的下角标char＝0，则表明该计算网格为铸锭网格但不为中心截面网格，采用下式计算温度场，获得三维方向铸造系统温度场分布：

如果计算网格(i，j，k)_char的下角标char＝21，则表明该计算网格为铸锭网格且为中心截面网格，采用下式计算温度场：

如果计算网格(i，j，k)_char的下角标char不为0且不为21，表明该网格不为铸锭网格，采用下式计算温度场：

[H]_m-char＝ρ_m-charc_m-charT_m-char

其中下角标m-char表示非铸锭材质；char可取值为2，4，5，6，7；c_m-char为比热(J/kgK)，ρ_m-char为密度(kg/m³)，λ_m-char为导热系数(W/m K)，T_m-char为温度(℃)，[H]_m-char为热焓(J/m³)，t为时间(s)。下角标in表示球墨铸铁，c_in为比热(J/kg K)，ρ_in为密度(kg/m³)，λ_in为导热系数(W/m K)，Ti_n为温度(℃)，[H]_in为热焓(J/m³)，L_in-heat为合金潜热(J/kg)，为二维方向合金液流动速度的合速度(m/s)，T_L为液相线(K)，T_E为共晶温度(K)；为哈密顿算子 表示对T_in应用哈密顿算子，即

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

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

针对所有下角标char＝21的计算网格(i，j，k)_char＝21，计算动量守恒方程，获得该计算网格中的金属液流动速度：

其中f_l为液相分数，U_z和U_y为二维截面上Z方向和Y方向上液体流动速度且0s时的值为0m/s，P为液相压强(Pa)，μ_l为液相粘度(Pa·s)，为重力加速度(m/s²)，β_T为温度膨胀系数(1/℃)，K_per为糊状区渗透率(m²)，λ_c为枝晶臂间距(m)。

通过本实施方式的各个公式可以对U_z和U_y进行求解，求解得到的值即为该计算网格中的金属液流动速度。

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

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤四中，针对所有下角标char＝0和21的计算网格(i，j，k)_char，计算石墨球尺寸。具体为：

步骤四一、针对所有下角标char＝0和21的计算网格(i，j，k)_char，计算过冷度(ΔT)下的形核率J_a。如果J_a＞0，则形核已完成，不计算形核率。如果J_a＝0，则采用下式进行计算。

其中ΔT为过冷度(℃)；J_a为形核率(m^-3s^-1)，0s时计算网格的形核率均为0；ε为石墨核心与液相浸润程度，为0.1～1之间的一个随机数。

步骤四二、计算球状石墨生长半径。如果J_a＞0，则

N^gra＝J_a·Δt

dR_G＝V_growth·Δt

其中N^gra为形核密度(m^-3)，V_growth为石墨球生长速度(m/s)，ΔT_growth为石墨球生长所需界面动力学过冷(℃)，R_G为石墨球半径(μm)其初始值为0.1μm，dR_G为Δt内石墨球半径变化值(μm)，g^gra为石墨球所占体积分数且初始值为0。上角标中的t+Δt和t分别表示当前时刻和上一时刻。即g^gra-(t+Δt)代表当前时刻的g^gra，g^gra-t代表上一时刻的g^gra。

当计算网格(i，j，k)_char所对应的形核率J_a＞0且表明石墨核心即将进入生长状态，初始状态时界面动力学过冷ΔT_growth为10^-3℃；当J_a＞0且表明石墨核心已经处于生长状态，则ΔT_growth由计算且为上一时刻所得值。当计算网格(i，j，k)_char所对应的g^gra大于或等于1时，表明该网格石墨球生长结束。

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

<实施例>

使用表1和表2提供的参数进行实验：

表1

表2

本实施例的铸造系统三维示意图如图1所示。图1中铸造系统尺寸为：其中Z_min＝0m，Z_max＝1.1m；Y_min＝0m，Y_max＝0.656m；X_min＝0m，X_max＝0.656m。网格剖分步长为0.004m，即Δx＝Δy＝Δz＝0.004m。

铸锭三维示意图如图2(a)所示，铸锭中心截面图如图2(b)所示，P1和P2是在中心界面上选取的两个相同高度的点。P1(X＝0.33m，Y＝0.184m，Z＝0.74m)，P2(X＝0.33m，Y＝0.31m，Z＝0.74m)。

图3(a)为铸锭温度场二维模拟、液体流动仅针对中心截面二维模拟时的P1点和P2点的冷却曲线，图3(b)为铸锭温度场三维模拟、液体流动仅针对中心截面二维模拟时P1点和P2点的冷却曲线。对比图3(a)和图3(b)，图3(b)中针对温度场进行三维计算、仅针对中心二维截面计算流场，P1和P2点处冷却速度快于图3(a)中P1和P2点冷却速度。P1点冷却曲线出现明显斜率变化。图3(a)中P1点0s～1492s，平均冷却速度为0.04℃/s；1492s～15865s，平均冷却速度为0.006℃/s。P2点0s～15865s，平均冷却速度为0.004℃/s。图3(b)中P1点0s～1817s，平均冷却速度为0.06℃/s；1817s～15865s，平均冷却速度为0.0095℃/s。P2点0s～15865s，平均冷却速度为0.0067℃/s。由此可见三维温度场计算所得冷却速度快于二维温度场计算所得冷却速度，这一点已经被实验证实。

图4为铸锭温度场二维模拟、液体流动仅针对中心截面二维模拟时的铸锭中心截面温度场分布图，图4(a)为4805s时刻的温度场分布，图4(b)为12005s时刻的温度场分布。

图5为铸锭温度场三维模拟、液体流动仅针对中心截面二维模拟时的铸锭中心截面温度场分布图，其中图5(a)为4805s时刻的温度场分布，图5(b)为12005s时刻的温度场分布。

对比图4和图5，由温度场分布图可以看出三维温度场计算条件下相同时刻下中心截面等温线所对应温度低于二维温度场计算所得结果。这一点符合实验研究和理论分析。

图6(a)为铸锭温度场二维模拟、液体流动仅针对中心截面二维模拟时P1点和P2点处石墨球半径随时间变化曲线，P1点和P2点处实验测量结果。图6(b)为铸锭温度场三维模拟、液体流动仅针对中心截面二维模拟时P1点和P2点处石墨球半径随时间变化曲线，P1点和P2点处实验测量结果。

对比图6(a)和图6(b)，相同时刻下图6(b)中石墨球半径低于图6(a)中石墨球半径。图6(a)中15865s时，P1点处和P2点处的石墨球半径分别为8986um和10798um。图6(b)中15865s时，P1点处和P2点处的石墨球半径分别为7836um和9344um。图6(b)中计算值与实验测量吻合较好。

由以上实施例可以看出，本发明的仿真计算结果与实验结果非常相近，具有较高的预测准确率。

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

Claims (5)

1.一种球墨铸铁铸锭石墨球尺寸数值预测的方法，其特征在于，包括：

步骤一、对铸造系统进行网格划分；

步骤二、对所有网格计算能量守恒方程，获得三维方向铸造系统温度场分布；

步骤三、对于铸锭中平行于重力方向的中心截面网格计算动量守恒方程，获得该计算网格中的金属液流动速度；

步骤四、对于铸锭网格和铸锭中平行于重力方向的中心截面网格，计算石墨球尺寸；

步骤五、重复步骤二至四，直至所有铸锭网格和铸锭中平行于重力方向的中心截面网格的温度T_in均小于共晶线温度T_E。

2.根据权利要求1所述的球墨铸铁铸锭石墨球尺寸数值预测的方法，其特征在于，步骤一具体为：

对X米×Y米×Z米的铸造系统进行宏观尺度网格剖分，X方向、Y方向和Z方向分别采用Δx、Δy、Δz作为网格剖分步长，Δx米＝Δy米＝Δz米，Δx、Δy和Δz的取值范围为1×10^-3米～4×10^-3米，计算网格的标号为(i，j，k)_char，其中i，j和k均为整数，i的取值范围是1～L，j的取值范围是1～M，k的取值范围是1～N， 下角标char＝2表示铸型网格，下角标char＝0表示铸锭网格，下角标char＝21表示铸锭中平行于重力方向的中心截面网格；铸造系统在X轴、Y轴、Z轴方向上的最小值分别为X_min、Y_min、Z_min，在X轴、Y轴、Z轴方向上的最大值分别为X_max、Y_max、Z_max。

3.根据权利要求2所述的球墨铸铁铸锭石墨球尺寸数值预测的方法，其特征在于，步骤二具体为：

若计算网格(i，j，k)_char的下角标char＝0，则表明该计算网格为铸锭网格但不为中心截面网格，采用下式计算温度场：

如果计算网格(i，j，k)_char的下角标char＝21，则表明该计算网格为铸锭网格且为中心截面网格，采用下式计算温度场：

如果计算网格(i，j，k)_char的下角标char不为0且不为21，表明该网格不为铸锭网格，采用下式计算温度场：

[H]_m-char＝ρ_m-charc_m-charT_m-char

其中下角标m-char表示非铸锭材质；c_m-char为比热，ρ_m-char为密度，λ_m-char为导热系数，T_m-char为温度，[H]_m-char为热焓，t为时间；下角标in表示球墨铸铁，c_in为比热，ρ_in为密度，λ_in为导热系数，T_in为温度，[H]_in为热焓，L_in-heat为合金潜热，为二维方向合金液流动速度的合速度，T_L为液相线，T_E为共晶温度；为哈密顿算子。

4.根据权利要求3所述的球墨铸铁铸锭石墨球尺寸数值预测的方法，其特征在于，步骤三具体为：

针对所有下角标char＝21的计算网格(i，j，k)_char＝21，计算动量守恒方程，获得该计算网格中的金属液流动速度：

其中f_l为液相分数，U_z和U_y为二维截面上Z方向和Y方向上液体流动速度且0s时的值为0m/s，P为液相压强，μ_l为液相粘度，为重力加速度，β_T为温度膨胀系数，K_per为糊状区渗透率，λ_c为枝晶臂间距。

5.根据权利要求4所述的球墨铸铁铸锭石墨球尺寸数值预测的方法，其特征在于，步骤四具体为：

步骤四一、针对所有下角标char＝0和21的计算网格(i，j，k)_char，计算过冷度ΔT下的形核率J_a；若J_a＞0，则形核已完成，不计算形核率；若J_a＝0，则采用下式进行计算形核率J_a：

其中ΔT为过冷度；J_a为形核率；ε为石墨核心与液相浸润程度；

步骤四二、计算球状石墨生长半径，若J_a＞0，则

N^gra＝J_a·Δt

dR_G＝V_growth·Δt

其中N^gra为形核密度，V_growth为石墨球生长速度，ΔT_growth为石墨球生长所需界面动力学过冷，R_G为石墨球半径，dR_G为Δt内石墨球半径变化值，g^gra为石墨球所占体积分数；上角标中的t+Δt和t分别表示当前时刻和上一时刻；

当计算网格(i，j，k)_char所对应的形核率J_a＞0且表明石墨核心即将进入生长状态；当J_a＞0且表明石墨核心已经处于生长状态，则ΔT_growth由计算且为上一时刻所得值；当计算网格(i，j，k)_char所对应的g^gra大于或等于1时，表明该网格石墨球生长结束。