基于无网格伽辽金法的连铸坯凝固计算方法

王 宁，王旭东✉，蔡来强，姚 曼

1) 大连理工大学材料科学与工程学院，大连 116024 2) 辽宁省凝固控制与数字化制备技术重点实验室，大连 116024

结晶器内钢液温度几乎无法直接测量，数值模拟是研究结晶器内铸坯凝固行为的重要途径[1].随着研究问题的深入，基于网格划分的数值方法无法精确重构凝固坯壳形貌和液/固相区，对于裂纹和大变形等网格复杂变化的问题甚至无法求解.无网格方法用离散的节点代替网格结构，在自适应、求解裂纹和大变形等问题体现出显著优势，已经应用于流体力学、电磁学、爆炸等[2−5]多个领域，也有少数研究将其应用于连铸过程.Zhang等[6]用有限点法和无网格局部彼得罗夫伽辽金方法分析了方坯凝固过程中的热弹塑性问题；Vertnik与Šarler[7]利用局部径向基函数配点法预测了钢液的湍流流动和传热行为；Yamasaki等[8]用基于粒子的MPS法模拟二冷区复杂的喷洒水流；Vaghefi等[9]用MLPG法分析了坯壳的热弹塑性情况；Hostos等[10]用无网格伽辽金方法(element-free Galerkin method，EFG)计算圆坯的传热和力学行为.但是与有限元、有限差分等模拟方法相比，无网格方法在连铸过程的研究仍然较少，且已有研究均采用均匀布置的节点离散计算区域，暂无针对无网格方法节点布置的灵活性及其与传统模拟方法相互验证的研究，而无网格方法最大的优势在于其节点布置的灵活性，因此探讨多种节点布置方法的计算精度是十分必要的.

本文依据EFG法基本原理和铸坯传热凝固控制方程，建立了基于EFG法的小方坯凝固过程数值模型，运用C++语言自行设计和开发了面向对象的铸坯凝固过程无网格数值计算软件.探讨多种节点布置方式下的无网格与有限元方法在计算精确度、自适应性、网格依赖性等方面的差异，为无网格方法求解连铸过程的传热、凝固以及后续的应力、多相耦合与界面追踪等问题提供基础.

1 铸坯传热/凝固模型的建立

依据铸坯的凝固特点，建立模型前作如下假设：沿拉坯方向传递的热量只占总热量的3%～6%，忽略铸坯纵向传热，简化为二维非稳态传热过程[11]；铸坯顶部和底部绝热；忽略弯月面以上保护渣吸收的热量；忽略结晶器的锥度和弧度以及铸坯凝固产生的收缩变形；铸坯热物性参数各向同性.

1.1 传热控制模型

连铸坯凝固过程二维非稳态传热方程[12]：

式中，k为导热系数，W·m−1·K−1；T和Tin分别为温度和初始浇铸温度，℃；T˙为温度对于时间的导数；ρ为密度，kg·m−3；c为比热容，J·kg−1·K−1；nx和ny分别为沿x和y方向边界外法线的方向余弦；q为结晶器/铸坯界面上的热流密度，W·m−2，采用Savage与Pritchard提出的经验公式[13]：,t表示时间，s，A和B为与结晶器有关的系数；

式中，Q[14]为内热源，J·m−3·s−1，fs为固相率，L为钢液的凝固潜热，J·kg−1.

1.2 变分处理

取温度变分δT为任意函数，利用分部积分法，建立能量泛函，得到等效积分弱形式为[15]：

将求解域Ω离散为有限的N个节点，求解域内t时刻任意位置(x,y)处的温度T(x,y,t)可用邻近的n个节点的温度与形函数的乘积加和来表示，得到式(4)：

式中，ΦI为节点I的形函数，将式（4）代入式(3)，整理得到式(5)：

通过求解式（5）可以得到连铸坯二维横截面的温度分布，沿浇铸方向，采用以下无条件稳定的向后差分法进行离散求解[15]，其中Δt表示时间步长：

2 无网格伽辽金法

2.1 最小二乘近似

1994年，Belytschko等在他们的重要文章里提出了EFG法[16]，其将移动最小二乘近似方法应用于伽辽金弱形式而产生一组代数方程，由于采取了伽辽金弱式，故需要一套与节点无关的背景网格进行数值积分.在现有的几十种无网格方法中，EFG法是应用最广泛的方法之一，相比于其他无网格方法，EFG法具有数值稳定、精度高的优点.它与有限元法的本质区别是：有限元法采用预定义的单元构造形函数，而EFG法利用如图1所示支持域中的节点和移动最小二乘法构造形函数[12].

图1 支持域、背景网格、积分点和场节点示意Fig.1 Support domains, background cells, evaluation points, and nodes

移动最小二乘法中，场变量u(x)近似表示为：

式中，pT(x)为基函数，m为其项数，对于二维计算域，线性基pT(x)=(1x y)，m=3，平方基pT(x)=(1x y x2x y y2)，m=6，函数a(x)是x的系数向量，aT(x)={a1(x)a2(x)···am(x)}.

利用式(7)构造加权残量泛函

式中，xi是x支持域内的节点，n表示支持域中的节点数，Wi(x)=W(x−xi)为权函数，它对移动最小二乘近似的性能至关重要，其大小与节点xi与采样点x之间的距离成反比例.

2.2 形函数

为了得到a(x)，对J求极值，即计算∂J/∂a=0得到：

将式（9）代入式（7），得到：

式中形函数[17]

2.3 高斯积分离散

在EFG法中，为了求解式(5)，将问题域离散成nc个背景单元，再利用高斯积分法求解每个背景单元中ng个积分点的数值积分[17]，即：

式中，wi是第i个积分点的加权因子，是对背景单元k第i个积分点处进行域内积分的雅克比矩阵，为对边界l第i个积分点处曲线积分的雅克比矩阵.nct和ngt分别为离散边界的曲线单元数和子曲线上的高斯点数.

3 模型验证与结果分析

3.1 实验条件和计算参数

以某钢厂断面为150 mm×150 mm的连铸小方坯为研究对象，工艺参数和铸坯热物性参数见表1，其中Kc为对流放大因子，取值范围为5～10；Tliq为液相线温度.

3.2 计算结果与讨论

3.2.1 节点均匀布置

有限元法的收敛性是指：①当单元尺寸不断减小时，有限元解逐渐趋近于精确解；②或当单元尺寸固定时，每个单元的自由度数越多，有限元解越趋近于精确解.由于无法测量结晶器内高温铸坯的实际温度，且难以计算铸坯非稳态传热模型的解析解，本文依据收敛性①，定义极细密单元分布下的有限元解为参考解[18].如图2所示，以距离弯月面367 mm处的铸坯横向切片为对象，计算不同单元数目下小方坯表面中心在浇铸10 s时的温度值，发现随着单元数量的增多，温度变化由开始的快速下降变为在1326.6 ℃附近缓慢降低，因此本文取90000 (300×300)个单元的ANSYS计算结果作为参考解，以验证后续不同无网格节点和有限元单元设置条件下计算结果的准确性和精度.

表1 工艺参数和铸坯热物性参数Table 1 Casting conditions and thermophysical properties of slab

为了考察基于EFG法的小方坯凝固传热模型的正确性和精确度，在小方坯横截面上用31×31个均匀分布的节点离散二维计算域，将其数值解与参考解和30×30个单元的ANSYS有限元数值解进行对比，对应的EFG法和有限元法的离散模型如图3所示.

图4为距离弯月面300 mm和结晶器出口处的温度分布，将EFG法数值解与相同节点数的ANSYS软件的数值解和参考解作对比，从温度分布看，EFG法的计算结果与有限元和参考解基本吻合，证明了EFG法计算小方坯凝固传热的可行性.

图2 单元数量对温度的影响Fig.2 Influence of the number of elements on temperature

计算并分析小方坯的角部节点和表面中心节点温度随浇铸方向的变化，比较EFG法和有限元法的精度，如图5所示，与相同单元数目的有限元法相比，EFG法的计算结果更加接近参考解，且其变化趋势也与参考解高度吻合.

3.2.2 节点加密布置

EFG法和有限元法中，节点和单元的数目越多，越有助于分析细节的变化，但节点和单元的增加也会导致计算成本的增加，因此综合考量计算精度和时间，可在高梯度区域插入较多节点和单元，以同时满足计算的效率和准确性.铸坯边界处温度梯度高，凝固行为、坯壳厚度将直接影响铸坯的质量，因此为了进一步考察铸坯的凝固进程，可在边界附近加密布置节点和单元.

有限元法及其他基于网格的方法需要重新划分网格，在处理疏密网格的界面时需要节点间的连接信息，前处理阶段耗费时间较长，处理繁琐.与之相比，无网格法如EFG法无需进行网格划分，无需节点间的连接信息，只需要在高梯度区域插入独立的离散点，即可得到精确度更高的数值模型.

图3 区域离散示意.（a） EFG （31×31节点）; （b） 有限元 （30×30单元）Fig.3 Schematics of discrete area: (a) EFG (31×31 nodes); (b) FEM (30×30 elements)

图4 不同高度铸坯横截面温度分布.（a） EFG （31×31节点）; （b） 有限元（30×30单元）；（c） 参考解 （300×300单元）Fig.4 Temperature distribution of billet at different heights: (a) EFG (31×31 nodes); (b) FEM (30×30 elements); (c) reference solution (300×300 elements)

图5 EFG法与有限元法数值解的对比.（a） 铸坯角部; （b） 表面中心Fig.5 Comparison between the numerical solutions obtained by EFG and FEM: (a) corner of slab; (b) midpoint of surface

本文探讨了EFG法和有限元法在边界加密条件下的数值解精度，小方坯二维计算域的离散模型如图6所示，EFG法在边界处加密1层，在非边界处稀疏布置，总体上仍然使用31×31个节点，有限元法采用ANSYS软件设定的边界单元加密方式，网格间的连续性使其无法添加独立的网格，因此边界附近的单元整体变动较大，共使用2056个单元.

图7展示了基于EFG法和ANSYS软件得到的小方坯表面温度分布，由角部附近的温度等高线位置可以看出，在边界处布置加密节点的情况下，基于EFG法的铸坯表面温度数值解低于有限元法.

为了进一步分析计算结果的差异，表2详细列出了距离弯月面不同高度处，在边界处布置加密节点的EFG法和有限元法的数值解与参考解的差值.可以看出，两种方法越接近结晶器出口处的数值解精度越高，同时EFG法数值解的精度明显高于有限元法，其与参考解的差值均在0.4 K以内.但是有限元法的单元数目远高于EFG法的总节点数961.需要指出的是，EFG法仅在边界处加密了节点，节点数目并未发生变化.以上结果说明，相比于有限元法，EFG法在节点灵活性、自适应以及解的精确度方面有显著优势.

图6 区域加密离散示意.(a) EFG (31×31节点); (b) 有限元 (2056单元)Fig.6 Distribution of local refinement near the boundary: (a) EFG (31×31 nodes); (b) FEM (2056 elements)

图7 铸坯表面温度分布.（a） EFG (31×31节点); （b） 有限元 (2056单元)Fig.7 Temperature distribution of billet surface: (a) EFG (31×31nodes); (b) FEM (2056 elements)

表2 两种方法数值解的对比Table 2 Comparison of the numerical results of two methods

3.2.3 节点随机布置

裂纹是造成铸坯报废的最主要缺陷，其产生与扩展均与结晶器内钢水的凝固行为密切相关，基于网格的数值模拟方法难以精确重构坯壳形貌和裂纹发展轨迹，而EFG法节点布置的灵活性使其可实现多种布点模式，如随机布置节点，可以减少由人为设置网格尺寸带来的误差，这在其他基于网格的数值方法中是难以实现的.为了验证EFG法对非规则布点模型的适应性，采用如图8所示的31×31和45×45个节点的随机布点方式计算小方坯的凝固行为.

图8 随机布点示意.（a） 31×31节点; （b） 45×45节点Fig.8 Distribution of random nodes: (a) 31×31 nodes; (b) 45×45 nodes

图9是利用EFG法布置不同数量随机节点条件下的结晶器出口处铸坯表面温度计算结果.可以看出，随机31×31节点得到的计算结果仍能保持一定的数值的稳定性和精度，但在个别位置光滑性较差，精度略低于ANSYS均匀30×30单元的数值解.与之相比，随机45×45节点模型数值解的光滑性有明显改善，与参考解吻合得较好，说明在节点有限加密的情况下，EFG法对非规则的随机布点的适应性较强，这种节点随机设置的方式在基于网格模拟方法中是难以做到的.

图9 结晶器出口处铸坯表面温度分布Fig.9 Temperature distribution of billet surface at mold exit

为进一步分析节点随机布置的精度，分别提取基于EFG法和ANSYS有限元软件的距弯月面不同高度铸坯角部的温度值，与参考解作差值并对比，结果如表3所示.节点数目为45×45的EFG的结果精确度不仅高于31×31的有限元法，还明显高于与其节点数目相近2056个单元的有限元边界加密模型的数值解精度.因此，可通过适度增加随机节点数目的方式来提高EFG法非规则模型数值解的精度，其结果显著优于单元数目相近的有限元模型的计算精度.

4 结论

（1）以结晶器内铸坯的传热和凝固过程为对象，建立并开发了基于EFG法的二维非稳态数值计算模型和专用计算程序.

（2）分别用EFG法和有限元ANSYS软件计算小方坯凝固过程，通过与细密单元的参考解对比，证实EFG法计算铸坯的凝固过程时具有以下优点：相同均匀分布的单元数目下，EFG法温度计算结果的精度比有限元法高；EFG法易于局部加密，仅需少数节点即可得到比有限元法更高精确度的数值解，适于跟踪梯度较大区域的数值变化情况；EFG法可通过适度增加随机节点的数目来提高模拟结果的光滑性与精度，回避了网格模拟方法在布置非均匀网格时网格过渡区域复杂的数据传递，在非规则随机节点设置方面显示出极大的灵活性.

表3 EFG 45×45随机布点与有限元法数值解的对比Table 3 Comparison of EFG 45×45 random nodes and FEM numerical results

（3）相比于有限元基于网格的数值模拟方法，EFG法回避了网格依赖性，具有良好的精度和自适应性，为进一步探索坯壳重构、液/固界面追踪、裂纹形成等网格模拟方法难以甚至无法处理的问题，提供了高效、可行的数值计算方法.