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

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

背景技术

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

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

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

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

技术实现要素：

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

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

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

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

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

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

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

附图说明

图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轴的选取可以依照实际情况而定。坐标系的选择不同不会影响预测结果。

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

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

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

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

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

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

对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轴方向上的最小值分别为xmin、ymin、zmin(米)，在x轴、y轴、z轴方向上的最大值分别为xmax、ymax、zmax(米)。

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

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤二中，针对所有计算网格(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-charcm-chartm-char

其中下角标m-char表示非铸锭材质；char可取值为2，4，5，6，7；cm-char为比热(j/kgk)，ρm-char为密度(kg/m³)，λm-char为导热系数(w/mk)，tm-char为温度(℃)，[h]m-char为热焓(j/m³)，t为时间(s)。下角标in表示球墨铸铁，cin为比热(j/kgk)，ρin为密度(kg/m³)，λin为导热系数(w/mk)，tin为温度(℃)，[h]in为热焓(j/m³)，lin-heat为合金潜热(j/kg)，为二维方向合金液流动速度的合速度(m/s)，tl为液相线(k)，te为共晶温度(k)；为哈密顿算子表示对tin应用哈密顿算子，即

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

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

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

其中fl为液相分数，uz和uy为二维截面上z方向和y方向上液体流动速度且0s时的值为0m/s，p为液相压强(pa)，μl为液相粘度(pa·s)，为重力加速度(m/s²)，βt为温度膨胀系数(1/℃)，kper为糊状区渗透率(m²)，λc为枝晶臂间距(m)。

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

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

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

针对所有下角标char＝0和21的计算网格(i，j，k)char，计算石墨球尺寸：

其中和分别为t+δt时刻和t时刻下的石墨球半径(μm)，t为0s时，δt为时间步长(s)，ρgra为石墨密度(kg/m³)，ggra为石墨球生长系数(m/s)，为液相-石墨界面处碳成分(wt％)，为热力学平衡状态下液相-石墨界面处碳成分(wt％)，为远处液相中碳成分(wt％)，为碳在液相中扩散系数(m²/s)，kpgra为石墨与液相之间碳元素分配系数。

采用迭代法计算δt时间步长内的和首先，

由r(ni)^*可得

其中ni为迭代次数；ni＝1表示δt时间步长内第一次迭代，此时ni＞1，r(ni+1)＝r(ni)^*。

如果的绝对值小于3×10，则表明计算收敛，用该值计算

如果的绝对值大于或等于3×10^-5，则表明计算不收敛，令r(ni+1)＝r(ni)^*，再次计算r(ni+1)^*、重复上述迭代计算，直到计算收敛；记录(t+δt)时刻下的最终获得随时间变化曲线。

针对下角标char＝21的计算网格(i，j，k)char，记录每个时刻t下的针对下角标char＝0的计算网格(i，j，k)char，记录每个时刻t下的

如果实验中所取试样位置为中心截面，则rg～时间曲线可以与实验结果进行对比。如果实验中所取试样位置为非中心截面，则o.8rg～时间曲线可以与实验结果进行对比。而在实验中，通常会在中心截面取试样，因为中心截面具有代表性，非中心截面不具备代表性。本发明在图6中给出的实验结果来自中心截面，因此图6中示出的是rg～时间曲线。如果进行了非中心截面的实验，也可以按照相似的方法绘制0.8rg关于时间的曲线。

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

<实施例>

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

表1

表2

本实施例的铸造系统三维示意图如图1所示。图1中铸造系统尺寸为：其中zmin＝0m，zmax＝1.1m；ymin＝0m，ymax＝0.656m；xmin＝0m，xmax＝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点处的石墨球半径分别为9199um和11000um。图6(b)中15865s时，p1点处和p2点处的石墨球半径分别为7923um和9466um。图6(b)中计算值与实验测量吻合较好。

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

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