CN115455760A - 铸坯在线热状态跟踪方法 - Google Patents

Abstract

本发明提供一种铸坯在线热状态跟踪方法，包括如下步骤：按照区块划分规则，对不对称的异形坯进行区域块划分，分别得到铸坯中心区域块和铸坯表面区域块；对铸坯中心区域块进行结构化的粗网格划分，对铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格，并对结构化网格的网格节点进行索引描述，将每个结构化块的网格数据以及索引对应存储；计算每个结构化块的网格单元的传热方程系数矩阵和凝固潜热释放源项；依次通过不完全矩阵LU分解和SIP算法耦合计算结构化块内的网格单元的温度场。利用本发明能够解决现有技术不适合用于不对称的异性坯与板坯两用铸机中，系统可靠性降低以及受时间步长限制等问题。

Description

铸坯在线热状态跟踪方法

技术领域

本发明涉及连铸加工技术领域，更为具体地，涉及一种铸坯在线热状态跟踪方法。

背景技术

针对复杂截面形状的异形连铸坯，不能采用结构化的正交网格进行温度场计算。早先，本项技术的发明者提出了采用显式格式的间断伽略金有限元方法进行铸坯热跟踪的方法，取得了发明专利(专利号，ZL201811231595.6)。在某钢厂实际应用中，取1/4异形连铸坯，采用三角形有限元网格进行连铸坯温度场计算。异形连铸坯在线热跟踪在国内首次投入生产应用。当时，可对两流异形坯进行生产控制。后来钢厂扩大产能，从两个铸流增加到三个铸流；但二级过程控制计算机并没有升级改造，仍然使用原来的二级控制软件，仅仅通过修改配置文件，就能适应三流异形坯在线热跟踪和动态配水的任务。但是，数值计算占用的CPU资源增加了，系统冗余不超过60％，虽然到目前使用中还没有发生故障，但降低了系统的可靠性。

该钢厂连铸机的另外一个特点是异形坯与板坯两用铸机，既能生产板坯也能生产异形坯。对于二级过程控制软件来说，就要能够同时计算板坯和异形坯的温度场。当时，采用了两种不同的方法，对于板坯采用结构化的正交网格，用有限体积方法进行数值计算；对于异形坯则采用非结构化的网格，用有限元方法进行计算。有限元方法计算程序本身就较为复杂，再加上有限体积方法，使得二级过程控制软件较为复杂，不便于后期维护和升级。随着钢铁工业竞争日益激烈，一机多用成为降低成本的可行方法。从该角度考虑，需要开发新的、统一的铸坯热跟踪的方法。

另外，该厂提出了用非对称异形坯轧制城市轨道(云轨)的设想，为此，本公司为该厂进行了连铸机设计、制定了生产工艺。考虑到非对称异形坯的几何结构和冷却状态，原来的异形连铸坯在线热跟踪系统将不再能使用。其主要原因有以下两个方面：一是，计算域扩大了，不能取1/4断面，至少要取 1/2。二是，非对称异形坯冷却更加不均匀，容易产生裂纹，需要更好地控制二次冷却，这就要求对铸坯划分更细的网格，对于显式格式的算法，计算时间步长受到稳定性条件的限制，网格细化后，计算效率将大幅度地降低，因为时间步长与网格尺寸的平方成反比。针对非对称异形坯采用什么计算网格和温度场计算方法需要仔细分析。

考虑到钢铁生产中常见的板坯和方坯，由于几何形状是对称的，若冷却条件也是对称的，则可以取1/4断面进行研究。目前，针对这类连铸坯,大多数温度场在线计算模型普遍采用显式格式的有限体积方法，取1/4断面，采用结构化的正交网格进行计算。这种计算方法，程序较为简单。对于大断面板坯，特别是一些特殊钢种，生产中容易出现表面和角部裂纹。考虑到铸坯表面喷淋冷却的不均匀性，需要采用更细密的网格，仔细刻画铸坯冷却边界条件，这就要对整个铸坯进行温度场计算，才能更好地控制铸坯冷却。与早期的模型相比较，计算域是原来的4倍。另外，若在铸坯表面附近和铸坯角部采用更细密的网格，那么，满足稳定性条件的计算时间步长会更小，早期模型的计算效率会受到很大影响。

因此，综上所述，目前现有技术中的铸坯在线热状态跟踪方法不适合用于不对称的异性坯与板坯两用铸机中，并且，现有技术中采用的网格划分方法由于无法具体考虑不对称异形坯的内部传热分布，导致系统可靠性降低；并且现有技术中采用的显式格式的有限体积方法计算连铸坯的温度场的时存在受时间步长限制的问题。

发明内容

鉴于上述问题，本发明的目的是提供一种铸坯在线热状态跟踪方法，以解决目前现有技术中的铸坯在线热状态跟踪方法不适合用于不对称的异性坯与板坯两用铸机中，并且，现有技术中采用的网格划分方法由于无法具体考虑不对称异形坯的内部传热分布，导致系统可靠性降低；并且现有技术中采用的显式格式的有限体积方法计算连铸坯的温度场的时存在受时间步长限制等问题。

本发明提供一种铸坯在线热状态跟踪方法，包括如下步骤：

S1、按照预设的区块划分规则，对不对称的异形坯进行区域块划分，分别得到铸坯中心区域块和铸坯表面区域块；

S2、对所述铸坯中心区域块进行结构化的粗网格划分，对所述铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格，并对所述结构化网格的节点进行索引描述，将每个结构化网格的网格数据以及索引对应存储；

S3、根据所述网格数据以及索引，分别对所述铸坯中心区域块内的每个结构化网格和所述铸坯表面区域块内的每个结构化网格的传热方程系数矩阵和凝固潜热释放源项进行计算；

S4、分别对所述铸坯中心区域块和所述铸坯表面区域块的传热方程系数矩阵进行不完全矩阵LU分解，并结合所述凝固潜热释放源项计算残差，分别得到所述铸坯中心区域块内的结构化网格的残差向量和所述铸坯表面区域块内的结构化网格的残差向量；

S5、根据所述铸坯中心区域块内的结构化网格的残差向量和所述铸坯表面区域块内的结构化网格的残差向量，采用SIP算法耦合计算所述铸坯中心区域块与所述铸坯表面区域块的温度场，得到多块的结构化网格的温度场，并保证区域块界面上的热平衡，通过迭代算法得到满足预设收敛条件的温度场；

S6、将根据所述满足预设收敛条件的温度场获取的每个结构化网格的温度值添加至所述网格数据中，并根据预设时间步长，返回步骤S3，向前推进时间，将每个结构化网格在不同时间的温度值与其对应的结构化网格的索引进行对应存储。

此外，优选的方案是，所述区块划分规则包括：使所述不对称的异形坯从铸坯表面指向铸坯内部的方向上保持厚度均匀，且铸坯中心区域块的边界由折线或圆弧曲线段围成。

此外，优选的方案是，所述对所述铸坯中心区域块进行结构化的粗网格划分，对所述铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格，并对多块的所述结构化网格的节点进行索引描述，将每个所述结构化网格的网格数据以及索引对应存储包括：

按照第一预设网格尺寸对所述铸坯中心区域块进行结构化的粗网格划分，按照第二预设网格尺寸对所述铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格；

标注所述结构化网格的节点的位置坐标以通过索引描述的方式确定所述结构化网格的位置，并将每个所述结构化网格的网格数据以及索引对应存储。

此外，优选的方案是，按照O-type或C-type的数据结构存储所述铸坯表面区域块的细网格划分后得到的结构化网格的网格数据，并建立所述铸坯中心区域块和所述铸坯表面区域块之间边界两侧的结构化网格的对应关系，将所述对应关系进行存储。

此外，优选的方案是，所述凝固潜热释放源项的计算公式为：

其中，

ρ为铸坯密度，g_l为液相率，L为凝固潜热，t为时间，T为温度；S_P为传热计算中凝固潜热释放源项经过线性化处理后的比例项，S_C为常数项，S_P和 S_C均为钢种凝固特性参数。

此外，优选的方案是，在分别对所述铸坯中心区域块内的每个结构化网格和所述铸坯表面区域块内的每个结构化网格的凝固潜热释放源项进行计算之前，还包括：

建立预设温度区间内各铸坯温度对应的钢种凝固特性参数的数据表格；在计算所述凝固潜热释放源项时，通过查找所述数据表格获取与当前铸坯温度所对应的凝固特性参数。

此外，优选的方案是，分别对所述铸坯中心区域块内的每个结构化网格和所述铸坯表面区域块内的每个结构化网格的传热方程系数矩阵进行计算，包括：

根据不包含所述铸坯中心区域块和所述铸坯表面区域块之间边界的结构化网格的内部节点的传热离散方程和所述包含所述铸坯中心区域块和所述铸坯表面区域块之间边界的内部边界传热离散方程，确定传热离散方程组；

根据所述传热离散方程组，确定传热方程系数矩阵。

此外，优选的方案是，所述不包含所述铸坯中心区域块和所述铸坯表面区域块之间边界的结构化网格的内部节点的传热离散方程为：

其中，

下标l为相邻的结构化网格的索引，l的取值范围为东E、西W、南S和北N四个相邻的结构化网格，下标P为当前结构化网格，A_P为系数，A_l为与当前结构化网格P相邻的结构化网格的矩阵系数，φ为求解变量；Q_p为载荷向量；

所述包含所述铸坯中心区域块和所述铸坯表面区域块之间边界的内部边界传热离散方程为：

对于铸坯中心区域块，所述内部边界传热离散方程为：

对于铸坯表面区域块，所述内部边界传热离散方程为：

其中，

P为包含铸坯中心区域块和铸坯表面区域块的内部边界的结构化网格，下标l的取值范围为不涉及内部边界的相邻单元，l的取值为西W、南S和北N，下标L和R分别表示内部边界两侧的左和右邻接单元，数组A_R和A_L分别为方程式的系数。

此外，优选的方案是，所述结构化网格的传热方程系数矩阵的不完全矩阵LU分解的公式为：

Aφ＝b；

LU＝A+N＝M；其中，A为传热方程系数矩阵，M为A的一个近似矩阵， N为M与A的差值矩阵，N满足Nφ≈0的条件。

此外，优选的方案是，在所述根据所述铸坯中心区域块内的结构化网格的残差向量和所述铸坯表面区域块内的结构化网格的残差向量，采用SIP算法耦合计算所述铸坯中心区域块与所述铸坯表面区域块的温度场，得到多块的结构化网格的温度场，并保证区域块界面上的热平衡，通过迭代算法得到满足预设收敛条件的温度场的过程中，

根据残差向量迭代计算温度场的公式为：(A+N)δⁿ⁺¹＝rⁿ；其中，

δⁿ⁺¹为温度场的修正值，rⁿ为第n次计算的残差向量，A为传热方程系数矩阵，N为差值矩阵；

SIP算法的公式为：(LU)δⁿ⁺¹＝rⁿ。

从上面的技术方案可知，本发明提供的铸坯在线热状态跟踪方法，通过将不对称的异形坯划分为铸坯中心区域块和铸坯表面区域块，再分别对铸坯中心区域块进行结构化的粗网格划分和对铸坯表面区域块进行结构化的细网格划分，解决了复杂形状的异形坯网格剖分问题，获得高质量的结构化网格，在铸坯表面附近传热剧烈的部分采用细网格，而在内部采用较粗的网格，可以提高计算效率，这种网格能够适用于所有连铸坯形，开发出的铸坯热跟踪技术能够满足多功能连铸机的使用要求，降低二级控制软件的维护和升级成本；选用隐式格式的有限体积方法计算连铸坯的温度场，解决了显式格式受时间步长限制的问题；对于特殊钢宽厚板坯和非对称异形坯，需要细化网格、更仔细地描述铸坯的传热和凝固，更严格控制铸坯均匀冷却，采用本发明提供的方法更为有利。

为了实现上述以及相关目的，本发明的一个或多个方面包括后面将详细说明的特征。下面的说明以及附图详细说明了本发明的某些示例性方面。然而，这些方面指示的仅仅是可使用本发明的原理的各种方式中的一些方式。此外，本发明旨在包括所有这些方面以及它们的等同物。

附图说明

通过参考以下结合附图的说明，并且随着对本发明的更全面理解，本发明的其它目的及结果将更加明白及易于理解。在附图中：

图1为根据本发明实施例的铸坯在线热状态跟踪方法的流程示意图；

图2为根据本发明实施例的不对称异形坯断面网格划分示意图；

图3为根据本发明实施例的板坯网格划分示意图；

图4为根据本发明实施例的一个区域块中的计算单元的有限体积计算网格示意图；

图5为根据本发明实施例的两个区域块之间边界左右两侧的内部边界单元的有限体积计算网格示意图；

图6为根据本发明实施例1的非结构化的FEM(有限单元法)计算网格的示意图；

图7为根据本发明实施例1的结构化的Multi-blocks FVM(多块有限体积) 网格的示意图；

图8为根据本发明实施例1的结晶器出口温度云图对比图，左侧是FEM，右侧是FVM；

图9为根据本发明实施例1的二冷区结束温度对比图，左侧是FEM，右侧是FVM；

图10为根据本发明实施例1的铸坯冷却曲线对比图。

在附图中相同的标号指示相似或相应的特征或功能。

具体实施方式

在下面的描述中，出于说明的目的，为了提供对一个或多个实施例的全面理解，阐述了许多具体细节。然而，很明显，也可以在没有这些具体细节的情况下实现这些实施例。

针对前述提出的技术中的铸坯在线热状态跟踪方法不适合用于不对称的异性坯与板坯两用铸机中，并且，现有技术中采用的网格划分方法由于无法具体考虑不对称异形坯的内部传热分布，导致系统可靠性降低；并且现有技术中采用的显式格式的有限体积方法计算连铸坯的温度场的时存在受时间步长限制等问题，提出了一种铸坯在线热状态跟踪方法。

以下将结合附图对本发明的具体实施例进行详细描述。

为了说明本发明提供的铸坯在线热状态跟踪方法，图1示出了根据本发明实施例的铸坯在线热状态跟踪方法的流程；图2示出了根据本发明实施例的不对称异形坯断面网格划分；图3示出了根据本发明实施例的板坯网格划分；图4示出了根据本发明实施例的一个区域块中的计算单元的有限体积计算网格；图5示出了据本发明实施例的两个区域块之间边界中的内部边界单元的有限体积计算网格；图6示出了根据本发明实施例1的非结构化的FEM (有限单元法)计算网格；图7示出了根据本发明实施例1的结构化的 Multi-blocksFVM(多块有限体积)网格；图8示出了根据本发明实施例1的结晶器出口温度云图对比，左侧是FEM，右侧是FVM；图9示出了根据本发明实施例1的二冷区结束温度对比，左侧是FEM，右侧是FVM；图10示出了根据本发明实施例1的铸坯冷却曲线对比，用粗实线绘制FEM的结果，用细实线绘制FVM的结果。

如图1至图10共同所示，本发明提供的铸坯在线热状态跟踪方法，包括如下步骤：

S1、按照预设的区块划分规则，对不对称的异形坯进行区域块划分，分别得到铸坯中心区域块和铸坯表面区域块。

作为本发明的一个优选方案，区块划分规则包括：使不对称的异形坯从铸坯表面指向铸坯内部的方向上保持厚度均匀，且铸坯中心区域块的边界由折线或圆弧等曲线段围成。

以截面最复杂的不对称异形坯为例，说明如何划分网格。如图2所示，该异形坯横断面左右对称，上下不对称，可以取其1/2作为研究对象。隐去该铸坯的断面尺寸，仅按照比例绘制断面。为了表述清楚，用粗网格划分铸坯截面。图2是划分的Multi-Blocks网格，若沿着腹板部分的中心线AB和翼缘中心部分的BCDEFB分界线将图形剪开，去除翼缘中心部分BCDEFB包围的区域，再将剩余的部分拉直，可以得到逻辑上的结构化网格(Block A)，这种网格称之为C-type网格。而另一区域则是另一个结构化的网格(Block B)。这样划分网格，可以获得结构化的网格。

为了说明所用技术的通用性，图3给出了完整的板坯横断面网格划分示例。铸坯表面附近的区域构成Block A，这种网格从逻辑上讲与圆环的网格划分一致，称之为O-type；铸坯中心部分为Block B。需要说明的是：所开发的有限体积温度场计算程序是通用的，既能处理非正交网格，当然也能处理正交网格；既能处理一个块，也能处理多个块。对于板坯、方坯可以像早期的程序一样只用一个块的正交结构化网格即可。

S2、对铸坯中心区域块进行结构化的粗网格划分，对铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格，并对结构化网格的节点进行索引描述，将每个结构化网格的网格数据以及索引对应存储。

作为本发明的一个优选方案，对所述铸坯中心区域块进行结构化的粗网格划分，对铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格，并对多块的结构化网格的节点进行索引描述，将每个结构化网格的网格数据以及索引对应存储包括：

按照第一预设网格尺寸对铸坯中心区域块进行结构化的粗网格划分，按照第二预设网格尺寸对铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格；

标注结构化网格的网格节点的位置坐标以通过索引描述的方式确定结构化网格的位置，并将每个结构化网格的网格数据以及索引对应存储。

作为本发明的一个优选方案，按照O-type或C-type的数据结构存储铸坯表面区域块的细网格划分后得到的结构化网格的网格数据，并建立铸坯中心区域块和铸坯表面区域块之间边界两侧的结构化网格的对应关系，将对应关系进行存储。

对每个区域块，Block A和Block B的节点用索引(i,j)来描述，网格数据存储十分方便。除此之外，还有两个显而易见的好处。其一，容易控制网格的疏密，合理划分网格。例如，在铸坯表面附近，由于温度变化较大，采用更细密的网格；而在铸坯内部，温度梯度较小，可采用较粗的网格。Block A 的选取(或者说，BCDEFB闭合曲线的选取)原则(区块划分规则)是：尽可能使得从铸坯表面指向内部的方向上保持厚度均匀，与铸坯凝固规律相一致。Block B的划分原则，不限于选择折线围成闭合曲线，也可选择圆弧等曲线段。特别是如图2所示，在CD段、FE段和CBF段，选择适当的圆弧曲线，使得划分出的网格更接近于正交网格，用专业术语来说就是偏斜度更小。总之，这样划分网格，可以得到更高质量的计算网格，控制整个传热离散控制方程的规模，优化和提高计算速度；其二，计算结果的显示，数据后处理。例如，若要显示异形坯的表面温度分布，分析铸坯冷却的均匀性，在铸坯表面部分采用C-type网格，就很方便，只需将第一层网格上的温度值取出，按照网格坐标绘制图形和曲线即可。若要显示铸坯表皮下几个mm的温度分布，只要按照预先划分好的网格取出对应层上的温度数据。然而，若要用非结构化的有限元网格，就不那么方便，需要插值计算。

S3、根据网格数据以及索引，分别对铸坯中心区域块内的每个结构化网格和铸坯表面区域块内的每个结构化网格的传热方程系数矩阵和凝固潜热释放源项进行计算。

作为本发明的一个优选方案，凝固潜热释放源项的计算公式为：

其中，

ρ为铸坯密度，g_l为液相率，L为凝固潜热，t为时间，T为温度；S_P为传热计算中凝固潜热释放源项经过线性化处理后的比例项，S_C为常数项，S_P和 S_C均为钢种凝固特性参数。

由于在拉坯方向上(Z轴方向)温度梯度非常小，相对于铸坯移动造成的传热，可以忽略该方向导热影响，仅考虑铸坯横截面上另两个(X和Y)方向的导热。若采用放大液态金属导热系数的方法处理钢液流动造成的传热效应，则常用的薄片移动法传热控制方程为：

上式中，κ_eff是当量导热系数，对于紊流区的液体取mκ_l，κ_l是钢液的导热系数，κ_l取常数，m是放大倍数，按经验取值；对于凝固后的坯壳，导热系数为固相导热系数κ_s，随温度变化。两相区的导热系数按照κ_eff＝(1-g_l)κ_s+g_lκ_l计算，g_l为液相率。C_p为比热。与导热系数的处理方法类似，两相混合物的比热按C_p＝(1-g_l)C_ps+g_lC_pl计算，其中，C_ps和C_pl分别是固相和液相的比热。(1) 式右端第二项为液相率随时间变化造成的潜热释放。

关于凝固潜热的处理，常用的方法是将源项线性化处理。

作为源项出现在方程(1)中，为了进行计算必须知道液相率g_l随温度T的变化规律，假定g_l与T满足线性关系：

式中T_s和T_l分别为固相线和液相线温度。将(2)式代入(1)式，则源项可以写为：

上式中T⁰和Δt分别为前一时刻的温度和计算时间步长。注意到：在将凝固潜热源项离散化的过程中采用了线性处理方法，并且在(3)式中的比例系数为负数。这样，可以保证在离散方程组系数矩阵中主对角元素占优，满足方程正定条件。

作为本发明的一个优选方案，在分别对铸坯中心区域块内的每个结构化网格和铸坯表面区域块内的每个结构化网格的凝固潜热释放源项进行计算之前，还包括：

建立预设温度区间内各铸坯温度对应的钢种凝固特性参数的数据表格；

在计算凝固潜热释放源项时，通过查找数据表格获取与当前铸坯温度所对应的凝固特性参数。在计算过程中根据表格编号索引快速查找凝固参数Sp 和Sc。

在假定液相率与温度服从线性关系，给出了凝固潜热的处理方法。实际上，铸坯凝固过程中，它们之间的关系不是线性的。根据连铸生产的钢种凝固特点，可以假设液相率与温度服从分段线性关系，在每个分段上，可以用上述方法进行处理。凝固潜热的释放给传热方程带来强烈的非线性特征，如果不能恰当地处理凝固潜热，那么，会造成隐式格式的方程组求解困难、甚至失败。将液相率与温度的曲线分段线性化，可以很好地处理凝固潜热的释放。但是，分段线性化处理就必须查找温度区间，知道当前的温度处于哪个凝固区间。若每次迭代都要判断分段区间，再计算出线段的斜率等参数，则会非常浪费CPU时间。为了加快计算，为每个计算节点分配一个整型量，用来记录凝固区间分段的位置索引i。参照公式(3)，根据钢种凝固规律将每个分段上的Sp和Sc等常量提前计算出来，建成数据表，在以后的计算中直接使用，避免重复计算。

为了减少查找凝固区间而花费的CPU时间，可以利用连铸坯凝固的特点。一般，铸坯温度随时间逐渐降低，可以利用这个特点优化查找的过程。假定温度T^*(Tⁱ＜T^*≤Tⁱ⁺¹)对应的液相率为gl^*(glⁱ＜gl^*≤glⁱ⁺¹)，属于凝固曲线上的第i个区段。为了查找T^*对应的gl^*分段数据表，若从头开始循环查找T^*所属的区段，则需要进行i次比较运算，非常浪费。可以利用凝固区间分段的位置索引i变量，若已经记录了计算节点前一时刻(或前一轮迭代)的温度T'^*在关系曲线上对应的位置i'，先判断当前温度T^*是否还处在第i'区段(T^i'＜T^*≤T^{i '+1})，如果满足条件，就可以直接引用数据表，这也是计算中经常出现的情况。若不满足条件，优先向温度降低的方向比较(T^i'-1＜T^*≤T^i')，倘若不满足温度降低的条件，则向相反的方向查找，这样只需经过少量的判断，就能完成查表。预先计算凝固参数数据表(Sp、Sc与T关系数表)，并利用凝固区间分段的位置索引i，采用优化路线查找参数，可以大幅度地提高计算效率。

作为本发明的一个优选方案，分别对铸坯中心区域块内的每个结构化网格和铸坯表面区域块内的每个结构化网格的传热方程系数矩阵进行计算，包括：

根据不包含铸坯中心区域块和铸坯表面区域块之间边界的结构化网格的内部节点的传热离散方程和包含铸坯中心区域块和铸坯表面区域块之间边界的内部边界传热离散方程，确定传热离散方程组；

根据传热离散方程组，确定传热方程系数矩阵。

作为本发明的一个优选方案，不包含铸坯中心区域块和铸坯表面区域块之间边界的结构化网格的内部节点的传热离散方程为：

其中，

下标l为相邻的结构化网格的索引，l的取值范围为东E、西W、南S和北N四个相邻的结构化网格，下标P为当前结构化网格，A_P为系数，A_l为与当前结构化网格P相邻的结构化网格的矩阵系数，φ为求解变量；Q_p为载荷向量；

包含铸坯中心区域块和铸坯表面区域块之间边界的内部边界传热离散方程为：

对于铸坯中心区域块，内部边界传热离散方程为：

对于铸坯表面区域块，内部边界传热离散方程为：

其中，

P为包含铸坯中心区域块和铸坯表面区域块的内部边界的结构化网格，下标l的取值范围为不涉及内部边界的相邻单元，l的取值为西W、南S和北N，下标L和R分别表示内部边界两侧的左和右邻接单元，数组A_R和A_L分别为方程式的系数。

如图4结合图5所示，在有限体积方法中，必须计算单元界面上的参数φ_e和其梯度值，可以采用线性插值近似计算。例如，界面上e点的参数φ_e，按照以下公式计算：

φ_e≈φ_e'＝φ_Eλ_e+φ_P(1-λ_e)

其中，插值系数λ_e根据几何参数按照以下公式进行计算，其中r_e，r_E和r_P是几何点矢量：

当计算网格偏斜严重时，e和e'点之间存在较大的差距，界面上e点的参数φ_e需要按照以下公式进行修正计算：

φ_e≈φ_e'+(gradφ)_e'·(r_e-r_e')

为了避免数值震荡，界面上e点的梯度按照下式计算：

其中

取前一个时间步长的计算结果，这一项作为源项出现在离散方程式中。根据有限体积方法，可以用单元界面上的变量φ和其梯度参数，推导得出离散方程式：

为了方便，对于块界面上的网格，不用E,W,S和N方式进行编号，用L和 R的方式编号，记录相邻关系，用A_R和A_L替代原来的A_E等来表示离散方程式系数，专门保存在特别的数组中。A_R和A_L的计算方法与A_E等离散方程式系数的计算方法相同。

S4、分别对铸坯中心区域块和铸坯表面区域块的传热方程系数矩阵进行不完全矩阵LU分解，并结合凝固潜热释放源项计算残差，分别得到铸坯中心区域块内的结构化网格的残差向量和铸坯表面区域块内的结构化网格的残差向量。

作为本发明的一个优选方案，结构化网格的传热方程系数矩阵的不完全矩阵LU分解的公式为：

Aφ＝b；

LU＝A+N＝M；其中，A为传热方程系数矩阵，M为A的一个近似矩阵， N为M与A的差值矩阵，N满足Nφ≈0的条件。

S5、根据铸坯中心区域块内的结构化网格的残差向量和铸坯表面区域块内的结构化网格的残差向量，采用SIP算法耦合计算铸坯中心区域块与铸坯表面区域块的温度场，得到多块的结构化网格的温度场，并保证区域块界面上的热平衡，通过迭代算法得到满足预设收敛条件的温度场。

作为本发明的一个优选方案，在根据铸坯中心区域块内的结构化网格的残差向量和铸坯表面区域块内的结构化网格的残差向量，采用SIP算法耦合计算铸坯中心区域块与所述铸坯表面区域块的温度场，得到多块的结构化网格的温度场，并保证区域块界面上的热平衡，通过迭代算法得到满足预设收敛条件的温度场的过程中，

根据残差向量迭代计算温度场的公式为：(A+N)δⁿ⁺¹＝rⁿ；其中，

δⁿ⁺¹为温度场的修正值，rⁿ为第n次计算的残差向量，A为传热方程系数矩阵，N为差值矩阵；

SIP算法的公式为：(LU)δⁿ⁺¹＝rⁿ。

具体的，根据SIP算法在结构化网格上求解温度场的迭代方程包括：

(A+N)φ＝(A+N)φ+(b-Aφ)

在方程得到收敛解时，(b-Aφ)＝0，因而，可以采用下式迭代：

(A+N)φⁿ⁺¹＝(A+N)φⁿ+(b-Aφⁿ)；即：

(A+N)δⁿ⁺¹＝rⁿ

φⁿ⁺¹和φⁿ是第n+1和第n次迭代的结果，其差值δⁿ⁺¹是计算过程中的中间变量，是迭代计算的修正量。rⁿ是第n次计算的残差向量。

根据SIP算法有：

(LU)δⁿ⁺¹＝rⁿ

Uδⁿ⁺¹＝L^-1rⁿ＝Rⁿ

修正量δⁿ⁺¹按照下式计算：

δⁿ⁺¹＝U^-1Rⁿ

在SIP算法中，充分利用了上面对角带状矩阵的特点，计算效率高。

采用结构化网格除了存储数据简便外，也可以简化计算。例如，关于i方向上两个相邻的计算单元i和i+1之间的界面，对i单元来说是西界面，对 i+1单元来说是东界面，穿过此界面的热流率，对这两个计算单元来说数值相等，仅相差一个负号。因此，关于离散方程式系数，i单元的系数A_W与i+1 单元的系数A_E就相差一个负号，只需计算一次。另外，为了减少CPU计算量，关于界面插值计算。在计算程序中，采用了与计算流体力学SIMPLE方法类似的循环。这里仅做简要说明：在一个时间步长的推进中，有两个循环，外循环和内循环。在外循环中，通过计算生成离散方程式系数矩阵，并对系数矩阵进行不完全LU分解。在内循环中，通过SIP迭代求解温度场。因为，外循环所做的LU矩阵分解在内循环中保持不变，只做一次LU分解，故而计算效率很高。

S6、将根据满足预设收敛条件的温度场获取的每个结构化网格的温度值添加至网格数据中，并根据预设时间步长，返回步骤S3，向前推进时间，将每个结构化网格在不同时间的温度值与其对应的结构化网格的索引进行对应存储。

以下示例将对本发明予以进一步的说明，以便本领域人员更好的理解本发明的优点和特征。

实施例1：

为了证明本发明提供方法的可靠性，首先做程序验证，用有解析解的简单算例验证程序。给出非稳态半无限一维导热问题的解析解：初始温度场均匀，从0时刻开始起，分别限定第1、2和3类边界条件，求内部温度随时间的变化，解析解是由误差函数组成的超越函数。数值计算结果与理论解一致，绘制的曲线重合，计算表明，采用2至3mm的网格，就能满足连铸工程精度要求。

为了验证多块的非正交网格计算，采用稳态一维圆柱坐标下厚壁圆管导热算例进行验证，若圆管内、外表面温度保持恒定，那么圆管内部的温度分布服从自然对数函数。人为地将计算域化成两个圆环，对每个区域采用O-type 网格进行剖分，形成Multi-blocks结构化网格，数值计算与理论解一致。限于篇幅，不展示这些简单算例的结果。

下面以图6所示的连铸非对称异形坯为例做应用说明，计算所用参数非常接近真实工况条件。铸坯物性参数如下：导热系数30(W.m^-1.℃^-1)，比热 650.0(J.kg^-1.℃^-1)，密度7200.0(Kg.m^-3)，凝固潜热309000(J.kg^-1)，液相线温度1520(℃),固相线温度1480.0(℃)。冷却区划分和综合热交换系数见表1。连铸生产速度为0.6m/min，浇铸温度为1550℃。

表1

根据有解析解的简单算例与数值方法的比较，认为2-3mm的网格就能满足工程精度要求的实事，作为比较的基准，采用较细的(略小于3mm)均匀三角形网格进行有限元计算。如图6所示，共15242个节点，29823个单元。由于计算网格很多，使用直接解方程组的方法进行计算，计算工作量很大；因此，在CPU为I9的台式计算机上完成30分钟的连铸坯仿真计算，大约需要7小时。

若按照本发明提供的方法，将铸坯划分成如图7所示的网格，铸坯表层是Block A，靠近铸坯表面，厚度方向最小网格尺寸2mm，向铸坯内部按照等比级数逐渐加粗网格，最粗网格大约6mm，共计130×15个网格，加粗网格的比例系数范围为1.03-1.08。铸坯翼缘中心部分为Block B，划分成6×30个网格，最大网格尺寸约12mm。网格单元总数为2130(130×15+6×30＝2130)。相比于FEM计算，网格数大为降低，相差2个数量级。完成相同的计算，花费的CPU时间仅4分钟。

为了检查计算结果，将结晶器出口(距弯月面0.7m位置)和二冷区出口 (距弯月面10.06m位置)的温度场绘制出来进行比较，如图8和图9所示，温度分布一致。

如图10所示，取腹板表面中心(A点)，大R角(B点)，侧面中心(C 点)，最后凝固的热节点(D点)，共4个特征点，绘制铸坯冷却曲线，做进一步的对比，见图10。由于两个网格点不完全一致，凝固热节点(D点)的冷却曲线存在一些差距，这是因为FVM网格较粗；其余点数据重合。采用如前文描述的FVM网格，可以获得理想的计算结果，计算CPU时间大幅度降低，能够满足在线温度场的计算要求。

通过上述具体实施方式可看出，本发明提供的铸坯在线热状态跟踪方法，通过将不对称的异形坯划分为铸坯中心区域块和铸坯表面区域块，再分别对铸坯中心区域块进行结构化的粗网格划分和对铸坯表面区域块进行结构化的细网格划分，解决了复杂形状的异形坯网格剖分问题，获得高质量的结构化网格，在铸坯表面附近传热剧烈的部分采用细网格，而在内部采用较粗的网格，可以提高计算效率，这种网格能够适用于所有连铸坯形，开发出的铸坯热跟踪技术能够满足多功能连铸机的使用要求，降低二级控制软件的维护和升级成本；选用隐式格式的有限体积方法计算连铸坯的温度场，解决了显式格式受时间步长限制的问题；对于特殊钢宽厚板坯和非对称异形坯，需要细化网格、更仔细地描述铸坯的传热和凝固，更严格控制铸坯均匀冷却，采用本发明提供的方法更为有利。

如上参照附图以示例的方式描述了根据本发明提出的铸坯在线热状态跟踪方法。但是，本领域技术人员应当理解，对于上述本发明所提出的铸坯在线热状态跟踪方法，还可以在不脱离本发明内容的基础上做出各种改进。因此，本发明的保护范围应当由所附的权利要求书的内容确定。

Claims (10)

1.一种铸坯在线热状态跟踪方法，其特征在于，包括如下步骤：

S1、按照预设的区块划分规则，对不对称的异形坯进行区域块划分，分别得到铸坯中心区域块和铸坯表面区域块；

S2、对所述铸坯中心区域块进行结构化的粗网格划分，对所述铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格，并对所述结构化网格的节点进行索引描述，将每个结构化网格的网格数据以及索引对应存储；

S3、根据所述网格数据以及索引，分别对所述铸坯中心区域块内的每个结构化网格和所述铸坯表面区域块内的每个结构化网格的传热方程系数矩阵和凝固潜热释放源项进行计算；

S4、分别对所述铸坯中心区域块和所述铸坯表面区域块的传热方程系数矩阵进行不完全矩阵LU分解，并结合所述凝固潜热释放源项计算残差，分别得到所述铸坯中心区域块内的结构化网格的残差向量和所述铸坯表面区域块内的结构化网格的残差向量；

S5、根据所述铸坯中心区域块内的结构化网格的残差向量和所述铸坯表面区域块内的结构化网格的残差向量，采用SIP算法耦合计算所述铸坯中心区域块与所述铸坯表面区域块的温度场，得到多块的结构化网格的温度场，并保证区域块界面上的热平衡，通过迭代算法得到满足预设收敛条件的温度场；

S6、将根据所述满足预设收敛条件的温度场获取的每个结构化网格的温度值添加至所述网格数据中，并根据预设时间步长，返回步骤S3，向前推进时间，将每个结构化网格在不同时间的温度值与其对应的结构化网格的索引进行对应存储。

2.根据权利要求1所述的铸坯在线热状态跟踪方法，其特征在于，所述区块划分规则包括：

使所述不对称的异形坯从铸坯表面指向铸坯内部的方向上保持厚度均匀，且铸坯中心区域块的边界由折线或圆弧曲线段围成。

3.根据权利要求1所述的铸坯在线热状态跟踪方法，其特征在于，所述对所述铸坯中心区域块进行结构化的粗网格划分，对所述铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格，并对多块的所述结构化网格的节点进行索引描述，将每个所述结构化网格的网格数据以及索引对应存储包括：

按照第一预设网格尺寸对所述铸坯中心区域块进行结构化的粗网格划分，按照第二预设网格尺寸对所述铸坯表面区域块进行结构化的细网格划分，得到多块的结构化网格；

标注所述结构化网格的节点的位置坐标以通过索引描述的方式确定所述结构化网格的位置，并将每个所述结构化网格的网格数据以及索引对应存储。

4.根据权利要求3所述的铸坯在线热状态跟踪方法，其特征在于，按照O-type或C-type的数据结构存储所述铸坯表面区域块的细网格划分后得到的结构化网格的网格数据，并建立所述铸坯中心区域块和所述铸坯表面区域块之间边界两侧的结构化网格的对应关系，将所述对应关系进行存储。

5.根据权利要求1所述的铸坯在线热状态跟踪方法，其特征在于，所述凝固潜热释放源项的计算公式为：

其中，

ρ为铸坯密度，g_l为液相率，L为凝固潜热，t为时间，T为温度；S_P为传热计算中凝固潜热释放源项经过线性化处理后的比例项，S_c为常数项，S_P和S_C均为钢种凝固特性参数。

6.根据权利要求5所述的铸坯在线热状态跟踪方法，其特征在于，在分别对所述铸坯中心区域块内的每个结构化网格和所述铸坯表面区域块内的每个结构化网格的凝固潜热释放源项进行计算之前，还包括：

建立预设温度区间内各铸坯温度对应的钢种凝固特性参数的数据表格；在计算所述凝固潜热释放源项时，通过查找所述数据表格获取与当前铸坯温度所对应的凝固特性参数。

7.根据权利要求3所述的铸坯在线热状态跟踪方法，其特征在于，分别对所述铸坯中心区域块内的每个结构化网格和所述铸坯表面区域块内的每个结构化网格的传热方程系数矩阵进行计算，包括：

根据不包含所述铸坯中心区域块和所述铸坯表面区域块之间边界的结构化网格的内部节点的传热离散方程和所述包含所述铸坯中心区域块和所述铸坯表面区域块之间边界的内部边界传热离散方程，确定传热离散方程组；

根据所述传热离散方程组，确定传热方程系数矩阵。

8.根据权利要求7所述的铸坯在线热状态跟踪方法，其特征在于，所述不包含所述铸坯中心区域块和所述铸坯表面区域块之间边界的结构化网格的内部节点的传热离散方程为：

其中，

下标l为相邻的结构化网格的索引，l的取值范围为东E、西W、南S和北N四个相邻的结构化网格，下标P为当前结构化网格，A_P为系数，A_l为与当前结构化网格P相邻的结构化网格的矩阵系数，φ为求解变量；Q_p为载荷向量；

所述包含所述铸坯中心区域块和所述铸坯表面区域块之间边界的内部边界传热离散方程为：

对于铸坯中心区域块，所述内部边界传热离散方程为：

对于铸坯表面区域块，所述内部边界传热离散方程为：

其中，

P为包含铸坯中心区域块和铸坯表面区域块的内部边界的结构化网格，下标l的取值范围为不涉及内部边界的相邻单元，l的取值为西W、南S和北N，下标L和R分别表示内部边界两侧的左和右邻接单元，数组A_R和A_L分别为方程式的系数。

9.根据权利要求8所述的铸坯在线热状态跟踪方法，其特征在于，所述结构化网格的传热方程系数矩阵的不完全矩阵LU分解的公式为：

Aφ＝b；

LU＝A+N＝M；其中，A为传热方程系数矩阵，M为A的一个近似矩阵，N为M与A的差值矩阵，N满足Nφ≈0的条件。

10.根据权利要求9所述的铸坯在线热状态跟踪方法，其特征在于，在所述根据所述铸坯中心区域块内的结构化网格的残差向量和所述铸坯表面区域块内的结构化网格的残差向量，采用SIP算法耦合计算所述铸坯中心区域块与所述铸坯表面区域块的温度场，得到多块的结构化网格的温度场，并保证区域块界面上的热平衡，通过迭代算法得到满足预设收敛条件的温度场的过程中，

根据残差向量迭代计算温度场的公式为：(A+N)δⁿ⁺¹＝rⁿ；其中，

δⁿ⁺¹为温度场的修正值，rⁿ为第n次计算的残差向量，A为传热方程系数矩阵，N为差值矩阵；

SIP算法的公式为：(LU)δⁿ⁺¹＝rⁿ。

Images