连铸过程中板坯凝固的数值模拟(1)

1、.连铸过程中板坯凝固的数值模拟数学模型的建立水冷钢液面结晶器控制方程凝固壳图5-1 铸坯凝固及坐标系从结晶器弯月面处沿板坯中心取一长为dx、高为dy的微元体，把微元体定位在二维坐标系中，与板坯一起向下运动。对微元体做热平衡，根据控制微元体内能量守恒关系，可推导出描述热量传输的基本方程：（5-1）式中，温度；时间；导热系数；密度；比热容；内热源。定解条件(1)初始条件：(浇注温度)(2)边界条件：结晶器：在结晶器冷却区，结晶器器壁的热流密度是拉坯速度和距钢水液面距离的函数，这一函数是由savage和pritchard提出并由Davies实验确定。(5-2)式中：热流密度；拉速；距钢水液面距离。二

2、冷区： ()(5-3)式中：对流换热系数；板坯表面温度；冷却水温度；空冷区：()(5-4)式中，物体黑度；韦斯蒂芬-波尔兹曼常数；外部环境温度。本项目中涉及的板坯无空冷区。内热源项的消除与凝固潜热的等效对于同一铸件液相的内能大于固相的内能，因此当合金凝固由液相转变为固相时，必然要产生的内能变化，这个内能变化量（通常用表示）成为凝固潜热。潜热释放是凝固过程区别于一般导热过程的显著特点。由于潜热释放，会明显降低铸件的冷却和凝固速度，傅立叶方程由于潜热的释放,实际变成了具有内热源的温度场，如图5-1所示。在实际使用的铸造合金中，通常不产生平滑的凝固界面，几乎所有的铸造合金都是树枝状结晶，在结晶过程中

3、要不断释放出结晶潜热。因此，在铸件凝固过程中，可将由于固相率增加所释放出的潜热看作相当于散失的热量，这样一来再处理合金凝固问题的数值求解时就可使问题变得简单多了。假定单位体积、单位时间内固相率的增量为，因此潜热放出的发热量为，此时的传热方程变为：(5-5)式中，密度，；比热容，；潜热，；固相率；温度，；时间，；导热系数，；二维坐标，。对于实用多元和金，要确定固相率和温度的关系，通常可以先采用热分析法求出凝固开始温度（液相线温度）和结束温度（固相线温度），假定如下：(1) 假定为线性分布时(5-6)(2) 假定为2次函数分布时(5-7)在实际处理潜热问题中，知道了固相率和温度关系之后，可采用以下

4、几种常用的潜热处理方法：等价比热容法、温度回升法和热焓法。等价比热容法比热容是指单位质量物体降低单位温度所释放的热量。单位质量金属在凝固温度范围内降低单位温度所释放的热量也可以理解成比热容，实际上这个比热容包括两部分即物体的真正比热容和潜热引起的比热容。将代入式(5-5)可得整理可得(5-8)=，即为等价比热容或有效比热容(亦称当量比热容)。在实际计算中假设潜热均匀释放，则（固相率采用线性分布函数）在温度区间内用置换，再根据固相率和温度的关系，将式(5-8)按非稳态到热处理。该方法适用于处理结晶温度范围较宽的合金，能满足绝大多数铸造合金的要求，目前应用较为广泛。温度回升法凝固开始的一段时间内，

5、固相不断增多，但温度基本上保持在熔点附近，这是由于所释放的潜热补偿了传热所引起的温度的下降。由于热量的多少常以单元体的温度变化来表示，因此可将这部分热量折算成所能补偿的温度下降，加到温度计算中去。假定某个领域（体积）的液态金属凝固时固相率增加，其释放出的潜热（被夺走的热量）可表示为(5-9)温度回升法中，先不考虑潜热的释放进行温度计算，求出微小时间内的温度降低量，为计算初始温度。如果，就产生凝固，由于释放潜热，使温度回升到（假定没有过冷），因此下式成立： (5-10)由式(5-9)与式(5-10)等量置换可得(5-11)此法采用固相率的增加来代替前热的放出，如果，则表明该领域的凝固结束。热焓法

6、凝固过程金属的焓可定义为(5-12)上式对温度求导，可得(5-13)将式(5-13)代入式(5-11)即得(5-14)这种方法与等价比热容法类似，适用于有一定结晶温度范围的合金。对比上述的三种潜热处理方法可知：温度回升法的思路比较明确，实现也比较简单，但是在程序设计要分配大量的存储空间，如果希望程序有断点续算功能则会生成大量的临时文件，而且温度回升法适应的范围小。热焓法和等价比热容法适应面要广得多，但是热焓法引入了另一个物理量(焓)，计算复杂，所以等价比热容法相对来说更经济，物理意义更明确。凝固液芯对流换热的导热等效凝固液芯对流换热的导热等效考虑液芯区对流换热对板坯传热的影响，采用提高液芯区内

7、钢液导热系数值，使液芯内的导热量计算值高于实际值，多出来的导热量作为对流换热的等效补偿，为了简化计算，实际的导热系数可看作是静止钢液导热系数的整数倍，表达式为： ()（5-15)式中，补偿对流换热的等效导热系数； 经验常数，钢液导入区，； 静止钢液的导热系数。热物性参数的数值处理方法由于导热系数、比热和密度是随温度变化的，其处理方法一般有常数法、线形函数法和插值法等，在不同的温度区间内用二次曲线拟合的方法得到温度-导热系数、温度-比热之间关系的曲线，来确定导热系数、比热值。对于板坯密度的取值，可在不同相区内取不同的值。(1) 固相区：；(2) 液相区：；(3) 固、液两相区：。钢的液相线温度和

8、固相线温度取决于化学成分，与C、Si、Mn、P、S、Cu、Cr、Al等元素含量有关。钢的比热容与钢种、温度有关，一般来说，比热容随温度升高而增大。(1) 固相区：；(2) 液相区；(3) 固、液两相区：利用前面提到的潜热处理方法将凝固潜热转化作液固两相等效比热容：(5-16)钢的导热系数与钢种、温度有关。(1) 对固相区，导热系数一般视为常数：；对低碳钢，；(2) 对液相区，流动的钢液采用前面提到的凝固液心对流换热的导热等效的方法，一般相当于静止钢液导热系数的48，取6；(3) 对于固、液两相区采用：(5-17)式中，液相导热系数；固相导热系数；两相区等效导热系数。数值处理方法网格的划分内部节

9、点对于非稳态导热方程(5-1)的数值处理方法，可将板坯进行网格剖分，获得许多矩形网格单元控制体，采用内节点法，使节点位于控制体的中心，控制体内值的变化可以用相关节点的值代替，本项目采用有限容积法来进行离散化。图5-2 内节点法离散区域坐标变换本项目采用适体坐标的方法生成网格，适体坐标的网格生成问题，可以看成是一个边值问题。边值问题的求解是偏微分方程领域中的一个经典课题。现在常用椭圆形系统生成网格，已知：(1) 计算平面上方向的节点总数及节点位置。在计算平面上网格总是均匀划分的，一般取(或0.1或其它方便的数值)。(2) 物理平面计算区域边界上的节点设置，这种节点设置方式应反映出我们对网格疏密布

10、置的要求。所以我们需要找出计算平面上求解区域内的一点与物理平面上一点之间的对应关系，如图5-3所示。图5-3 适体坐标示意图从物理平面上来看，把看成是物理平面上被求解的因变量，则就构成了物理平面上的一个边值问题：即已经知道物理平面上与边界点相应的，要求出与内部一点对应的。在数学上描写边值问题的最简单的椭圆形方程就是Laplace方程。根据Laplace方程解的唯一性原理，可以把看作为是物理平面上Laplace方程的解：(5-18)这样在物理平面的求解区域边界上规定、的取值方法，就形成了物理平面上的第一类边界条件的Laplace问题。虽然对于这类问题已研究得很成熟，但是由于物理平面上是个不规则区

11、域，于是在物理平面上解这一问题又碰到了不规则边界的困难。换个角度来说，如果从计算平面上来看，从计算平面上的边值问题出发来考虑，则情况就大为改观，因为在计算平面上可以永远取成一个规则区域。所谓计算平面上的边值问题，就是指在计算平面的矩形边界上规定的取值方法，然后通过求解微分方程来确定计算区域内部各点的值，即找出与计算平面求解区域内各相应点的物理平面上的坐标。实际上用椭圆形方程来生成网格时都是通过求解计算平面上的边值问题来进行的。为此需要把物理平面上的Laplace方程转换到计算平面上以为变量的方程。所以把方程(5-18)转化为：(5-19)其中参数的定义为：(5-20)图5-4 板坯在物理平面上

12、的网格划分300´50式(5-2)连同计算平面求解区边界上已知的与之间的关系就构成了计算平面上第一类边界条件的边值问题，这样就变得比较容易求解。图5-5 物理平面上垂直段、圆弧段、水平段网格的局部放大300´50图5-6 计算平面上的网格划分（）控制方程的转换及其离散控制方程的转换本项目中只考虑扩散（导热）问题，不考虑对流问题，液态钢水的对流对温度场的影响通过凝固液芯对流换热的导热等效的方法给予简化。因此，二维非稳态导热的控制方程如方程(5-1)所示：在将物理平面的控制方程转换到计算平面的过程中，需要利用函数的导数与其反函数导数间的数学关系。设物理平面上坐标与间的关系为，则

13、物理平面上对于的偏导数与计算平面上对的偏导数之间的关系可用以下形式表示之：(5-21)其中J称为雅克比因子，它代表了计算空间中控制容积的膨胀程度。利用以上关系式可以得到计算平面上与式(5-1)相应的非稳态对流-扩散方程：(5-22)其中：如定义式(5-20)。计算平面的离散在计算平面的网格上离散通用控制方程(5-22)。计算平面上的同位网格如图5-6所示，在向，节点间的距离分别为，界面之间的距离分别表示为，对控制容积作积分，为了简化起见，令，并且用近似的代替，从而得到最终的离散形式为：(5-26)其中 系数取决于所采用的格式，本程序采用的格式有：(1)混合格式；(2)中心差分格式；(3)迎风格

14、式；(4) 乘方格式。但因为连铸是一个板坯连续凝固过程，有拉速()的存在。因此，我们在整个温度场计算平面上所有节点处都加入了一项横向(方向)的对流项，速度存在于温度场的计算中，只起到一个能量传递的作用，速度作为已知量处理。图5-7 计算平面上的网格计算平面上的速度需要由板坯的拉速转换而生成。（5-27）因此，引入所产生的对流项为：，（5-28a），(垂直拉坯方向速度为零)（5-28b）扩散项为：，（5-29ab），（5-29cd）为归结为源项的部分，写成如下形式：（5-30）上式中。计算方法采用有限差分法计算铸坯各离散点的温度和凝固状况。离散后的方程。其中空间离散采用有限容积法，时间项可由全隐

15、式差分表达，离散后得到的线性方程组采用交替方向隐式迭代法(ADI)结合解三对角矩阵的三对角矩阵算法 (TDMA)求解，全隐式条件下的物性参数通过欠松弛方法获得收敛解。当计算区域的边界条件为第2或第3类边界条件时,边界上的温度是未知的。为使内部节点的温度代数方程组得以封闭,采用附加源项法将第2类或第3类边界条件所规定的进入或导出的热量作为与边界相邻节点的当量源项。连铸二冷制度的优化本项目采用的方法是目标表面温度反算法。首先考虑钢种的高温力学性能和前面所述冶金准则以及二冷配水的工艺要素，确定二冷区目标表面温度模型，然后以目标温度分布作为连板坯凝固传热模型中的二冷区边界条件，利用传热模型反算出二冷区

16、传热系数的分布，最后根据传热系数与喷雾水流密度的关系式计算出二冷区的水量分布。下面以为板坯为例，介绍该方法的应用。首先确定目标表面温度分布，满足两方面的要求：缩小表面温度回升和保证板坯表面温度范围。普碳钢或低合金钢的高温力学性能主要与钢中碳含量有关，因此通过不同钢种的高温热力学特性曲线找出板坯碳含量与塑性温度的对应关系，此温度即是板坯在二冷区终点处的表面控制温度。（5-31）（5-32）式中,为钢中碳的质量分数。板坯出结晶器后表面温度回升超过100/m会有内裂发生。根据该区域内板坯温度的上升与二冷区长度应为直线关系，该区域二冷各段末端板坯表面温度的计算公式如下:(，)式中，为结晶器出口处板坯表

17、面温度；为各段末端温度。设二冷区长度为，共段,，各段末端距结晶器出口处的距离为，二冷区最高温度点距结晶器出口处的距离为。表面温度下降区域温度控制在200/m以下,并与二冷区长度成直线关系，而且与碳含量有关。（5-33）( )（5-34）( )板坯表面其它点的目标温度以抛物线拟合形式给出。第二步求解方坯二维传热模型，其二冷区边界条件为:（5-35）利用数值方法求解出小方坯温度场。根据计算出的温度场，由下式求解出满足目标温度的二冷区传热系数：（5-36）式中:导热系数； 表面温度； 冷却水温度。 最后由M.Ishiguio等人推导出的与水流密度经验公式的关系式：（5-37）式中，换热系数，；水流密

18、度，；冷却水温度，。求出二冷区喷雾冷却的水流密度。由以上算法计算获得二冷区相对应的水流密度分布，再根据各段喷射面积，获得各段的喷水量。程序的编制本项目中板坯导热温度场模拟运算程序采用VisualC+6.0高级语言编制。程序采用交替方向隐式迭代法(ADI)结合解三对角矩阵的三对角矩阵算法(TDMA)进行迭代计算。图5-8为流程图。图5-8 流程图板坯凝固算例比较分析算例外形尺寸：结晶器钢液面高度：800mm；板坯尺寸： 板坯宽度：1200mm；板坯厚度：225mm；二冷区垂直段长度：1,996mm；二冷区弯曲段弧长：1,661mm；二冷区圆弧段弧长：12,585mm；二冷区矫直段弧长：4,340

19、mm；二冷区水平段长度：8,375mm；二冷区弯曲段等效半径：9,896.28mm；二冷区圆弧段半径：10,000mm；二冷区矫直段等效半径：9,224.9mm。弯曲段、矫直段外缘尺寸采用三次方曲线绘制。表5.1 钢种成分钢种/成分0.170.230.170.370.901.200.0350.0351.11.40.20.30.420.500.170.370.500.800.0350.035000.620.680.200.551.201.600.0400.04000表5.2 物性参数名称符号单位：数值钢的潜热(比内能)272,000固相区钢的导热系数34.886液相区静止钢液的导热系数29.1固相区钢密度7600固、液两相区钢密度7200固相区钢的比热容750.0液固相区钢的比热容800.0钢的液相线温度1514钢的液相线温度（凝固温度）1442中间罐内钢水过热度18.5浇注温度1532.5表5.3 其他连铸参数名称符号数值单位：二冷区冷却水温度30.拉速1.08选取钢用目标温度法进行模拟，计算出温度场以及铸坯表面中心温度分布如图5-9，5-10，5-11所示。图5-9 模拟出的温度场分布图5-10 温度场的局部放大图5-11 铸