基于无网格法的连续铸造凝固与热应力数值模拟张磊 申厚发 荣一鸣 黄天佑摘要：伴随着无网格法凝固模型的发展,一个完整的热机械分析系统逐渐成长起来。这种无网格法凝固模型基于有限点法（FPM）和以弹塑性分析模型为基本原理的无网格局部Petrov-Galerkin法（MLPG）。此系统用于模拟连续铸造时模中钢坯的凝固过程及热应力场。模拟结果与实际测量具有很好的一致性，这很好的说明了应力应变的分布特点和非拐角区域缺陷形成机制。当使用快速计算与动态节点分布时其计算精度与FPM相同。此处举了一个使用拉格朗日类型的MLPG法解决大变形的实例。结果显示连续铸造过程中不正确的冷却缩减规范是潜在裂纹问题的主要原因。实验显示无网格法将会成为研究连续铸造过程的一个有力分析工具，尤其对于机械收缩产生的大变形与开裂问题。关键词：数值模拟 凝固 热弹塑性 连续铸造 大变形一引言 连续铸造是世界上最引人瞩目的钢铁生产方法，它在诸多国家的钢铁行业中扮演着领军角色。连续铸造过程中，热量传递、金属凝固和凝固壳中应力应变发展将直接影响产品质量，甚至会加助裂纹的形成1，2。连续铸造通常在高于钢熔点的高温下进行，因此几乎无法测得温度值与应力值。而进行实物模拟不仅困难而且成本高。 凝固过程中的热传输现象和应力应变的发展可用方程组形式的数学模型描述（一般使用一系列不同或部分相同的等式表达）。模型中，应力和温度是待求未知量，采用有限点（FPM）或有限元（FEM）数值方法求解2，3。 在这些方法中，当定义一个局部不同的控制方程，通常把此局部区域划分为网格。此处的网格是指由相连接点以一定规则构成的网格之间的空间或间隙。通过在相邻的节点恰当的划分，网格可以用一系列不同的代数方程近似表达，然后联立所有领域中的代数方程就构成了整个系统的代数表达式4。 网格是有限元FEM和有限节点FDM法的基础，但它同时也限制了它们的发展空间。例如，FDM要求网格划分的均匀，而这在复杂边界条件的问题上非常困难。相反的，FEM在处理复杂边界问题上比FDM轻松但却不适用于大变形时网格畸变及裂纹扩展随机路径带来的网格重叠问题。而在连续铸造过程中经常产生大变形和裂纹问题。因此在设计时需要精确预测。在这种背景下，为成功模拟连续铸造过程，已有20年研究历史的无网格法再次受到关注。无网格法构建于节点近似理论基础上，它不需要划分网格，由于不需要描述节点之间的联系也无需节点连接方法，需要时只需在指定区域增加或减少节点数目即可。由于在解决大变形和裂纹扩展问题上具有很大优越性，无网格法在连续模拟具有很广阔的发展空间。 无网格法用在区域和节点随机分布的节点来描述问题区域及边界。首先根据节点分布创建近似的未知函数，再采用插值法（最小二乘法）使近似值和实际值之间的误差最小化，未知函数就可以用因具有不同控制方程与边界条件而具有不同残量的近似方程来表达。最后通过使整个区域及其边界的加权残量强解与弱解最小化解决问题。 到目前为止，已有十余种无网格法问世（见表1），传统的无网格法通常基于构造未知域近似函数也就是选择合适的函数来表达与解决。表1主要无网格法小结4，5名称待求系统方程解形式近似方案背景网格积分漫射元法弱形式最小二乘法配点，Galerkin法有需要无单元Galerkin法弱形式最小二乘法配点，Galerkin法有需要无网格局部Petrov-Galerkin法局部弱形式最小二乘法配点，Petrov-Galerkin法无需要有限点法强形式最小二乘法，特定元代求法（泰勒式）无不需要光滑质点流体动力学法强形式全局求法无不需要重构核粒子法弱形式或强形式全局求法（重构核粒子法）有需要Hp-云团法弱形式配点法，移动最小二乘法有需要FEM配点法弱形式配点法，移动最小二乘法有需要 在这些无网格法中，当需要加入背景网格时，问题的解决难度与复杂程度会大幅增加。SPH法不需要划分背景网格在解决无边界问题是比较有效，FPM和MLPG只需要节点是真正的无网格法。FPM相对较容易实现，但却存在数值稳定的问题。对于非自身接连问题，例如在流体机械中需要采用一定的措施确保数值近似相同；在传热问题中纽曼带边界条件属于区域对流，当使用加权函数代替非对称有限区域温度场时可能导致不精确及不稳定等问题，因此需要施加一些特殊的稳定化措施。但同样的稳定化方程并不适用于热弹塑性问题。与FPM相比，MLPG不容易实现，但却更精确稳定。在互不相连的模型中，MLPG可如FPM一样使用。此外由FPM法得到的温度可直接加在相应的热度的载荷上来解决弹塑性问题。本文研究的即是一个全面无网格热机械分析系统。它采用以FPM为基础的凝固模型和以MLPG为基础的高温弹塑模型来模拟热传递凝固和高温下连续铸造铸模的应力应变转化。系统预测结果与测量结果具有很好的一致性并能合理的表达应力应变分布特征。同时对于非拐角处裂纹的形成机制也能有一定的预测。在其他一些模拟上FEM达不到同等的精度。二、模型开发 一个完整的无网格热机械分析系统已经建立起来。 凝固过程无网格法分析模型基于FPM理论建立。Neumann边界条件采用Onate稳定化方法进行调整。弹塑性热应力无网格模拟基于MLPG理论，此系统用来模拟凝固过程和钢坯铸造的热应力分布情况。2.1形函数及其导数形函数用最小二乘法11构建。首先，将未知函数方程近似值定义为 式中pI（x）为以m为项数的多项式基函数。aI（x）为空间坐标函数x=x1，x2，x3T的未知参量。通过使用最小二乘法，a（x）通过使加权离散L2误差量最小化来确定： （2） 式中wI（x）为以节点I为紧支9，12的权函数，uI为未知节点值。矩阵P，W定义为： 且有 对函数J（x）的a求导并使其导数为0，得到J（x）的极小值： 解之得： 式中 A=PTWP， B=PTW （8）将a代入式（1）得 式中为MLS近似形函数。形函数的一阶倒数为13 式中A，i-1=-A-1A-1，kA-1，二阶导数为 式中 最小二乘近似式（9）可用于各种无网格法（如FPM，MLPG，EFG等）中来代替相应的未知域函数。MLS近似形函数及其导数可由式（10.1）(10.2)计算而得。对未知域函数的求解转化为对a的求解。在传热问题中，未知域函数为温度方程，在应力应变问题中未知域函数为位移方程。解平衡方程u并将其代入式（9）就得到未知域函数。2.2 凝固模型在对连续钢坯铸造过程进行传热与凝固分析时，假定系统处于持续稳定状态，Fig.1为连续钢坯铸造截面图。沿浇铸方向可将其简化为二维平面图。在不同的高度，边界上的通量大小由水冷能力决定。Fig.1b为相应的二维计算机模型，x为厚度方向，y为宽度方向，z为充型方向。若材料的焓H为已知，连续铸造过程的热平衡方程可表示为 式中T为温度，v为充型速度，为密度，k为热导率，H为焓函数在表面边缘也就是Fig.1b中所示的溶解域边界，热通量边界条件为 式中n为边界标准矢量。qn为已知边界热通量，由模中水冷能力决定。在模中弯液面初始温度为 其中Tpour为浇注温度。FPM用来解连续铸造中具有特殊边界条件的热平衡方程。溶解域被离散为节点（见Fig.1b）也就是一系列节点随机分布于溶解域内部及边界上（见Fig.2b），将这些点代入方程（11）-（13）并求解可得到以下方程组： 式中N为节点总数，为溶解域，为热通量边界。将MLS近似函数及其导数（见式（10.1）（10.2）代入式（4）中，将稳定项8加到Neumann边界条件上，则连续铸造凝固过程的最终FPM公式变为 式中T1h为T1的最小二乘近似函数（9），h为特征长度7，8.与其它无网格法相比，FPM方案更容易实现，当采用一定稳定化措施时刻获得很高精度。Fig.3为一个由上述FPM方案解决的凝固问题。结果证明FPM模拟方案非常高效并且达到的精度更高。2.3 变形与应力模型 连续铸造过程中凝固壳的变形和应力应变分布可以由热弹塑性应力应变关系及基于小变形的边界条件获得。在2D图的边缘薄片处，应力应变的发展由以下几个方程控制 力边界条件： ij，jnj-ti=0 在上 （16.2） 位移边界条件： ui=ui 在u上 （16.3）式中ij为应力张量，bi为域内已知体力，（），i代表()xi，ti为受拉边界t上给定的拉力，ui为位移边界u上给定位移，ni为边界上标准单元外表面。根据Mises屈服准则及Prandtl-Reuss流动条件来确定应力-应变关系。在弹性变形尺度内，应力增量与应变增量的关系为 =De（-T） （17.1）在塑性变形范围 =Dep（-T）+T （17.2）式中为全应力增量，为全应变增量，De为弹性刚度矩阵，Dep为弹-塑性刚度矩阵，T为热应变增量，T为热应力增量14.上述方程可采用无网格局部Petrov-Galerkin法解。在域内取一点x，以x为圆心做圆得到域的子域s，在此子域中取方程（16.1）全体的局部弱形式，并取1来保证本质边界条件，用离散定理来保证自然边界条件则方程（16）所代表的MLPG公式可变为 式中u为试探函数，t为检验函数，v代表由检验函数所得的应变矩阵，代表由试探函数的应力矢量；su表示s在位移边界的部分，具有指定本质边界条件；s为求解域边界部分具有指定自然边界条件（见Fig.4）此时可应用节点模型，将式（9）中所有节点取近似值u则应力应变的热弹塑性关系式（17）可加到（18）上，得连续铸造中凝固壳的热弹塑性MLPG离散系统： 式中 其中，I,J=1，2，N，D为刚度矩阵在变形相应范围内分别以De，Dep表示，N为单元表面标准矩阵，BJ为MLS近似方程的导数矩阵。S为边界条件矩阵9，10. 表2 制造参数模型尺寸（mm）115115圆角半径（mm）8弯液面（mm）100模长（mm）700充型速度（mmin）2.3钢含碳量（wt%）0.16浇注温度15203.结论及讨论 连续钢坯铸造热传递，凝固，凝固壳应力应变发展均可用无网格模型模拟。Fig.2a为模型几何图。在溶解域随机分布一系列节点，取空间步长1.2mm，共2395个节点，见Fig.2b，表2为对应制造参数，在凝固表面上分布7个节点，2个位于拐角处，2个位于表面上（见Fig2a）来检测温度分布。检测结果见Fig5，由图知凝固前沿冷却速度明显高于拐角处。在凝固前沿，热量迅速扩散，表面凝固为一固体壳，Fig6给出了距弯液面不同距离处凝固层厚度，并将之与FEM模拟结果及实际测量值进行了对比15。