CN102652963B - 一种中厚板轧后超快冷过程温度场的耦合控制方法 - Google Patents

Abstract

一种中厚板轧后超快冷过程温度场的耦合控制方法，属于轧钢自动控制技术领域，过程自动化系统采集钢板的PDI参数，采用温度场的耦合控制算法，对每一个水流密度层别分别进行温度场仿真，从而设定冷却规程，并下发给基础自动化系统执行该冷却规程，同时将自动采集到的轧线上仪表的实测数据上传给过程自动化系统，触发过程自动化系统进行冷却规程的修正和自学习。本发明方法实现了比热与温度场的耦合，有效解决了传统有限元法进行超快冷过程温度场仿真时所表现出的振荡问题；成功应用该温度场耦合控制方法进行中厚板超快冷过程的在线温度场解析，实现了冷却工艺的高精度设定。

Description

一种中厚板轧后超快冷过程温度场的耦合控制方法

技术领域

本发明属于轧钢自动控制技术领域，特别是涉及到一种中厚板轧后超快冷过程温度场的耦合控制方法。

背景技术

在中厚板轧后的冷却过程中，高精度的温度场模型不仅仅是实现终冷温度、冷却速度等单一控制目标的基础，更可以为目前许多研究者提出的冷却路径等新思想的生产应用奠定坚实的基础。在中厚板轧后冷却过程中进行温度场仿真的数值模拟方法主要有解析法、有限差分法和有限元法。由于中厚板轧后冷却过程中工业现场环境因素的时变性，使得基于解析法而构建的温度场模型的适应范围较小，不能满足模型在线应用的要求。而基于有限差分法而构建的温度场模型只注重了节点的作用，忽略了把节点联系起来的单元属性，即忽略了单元属性对各节点温度的“贡献”。而有限元法在继承有限差分法进行离散的核心思想的同时，还注重单元属性对节点温度的“贡献”，其求解精度更高。因此，人们越来越倾向于利用有限元法进行温度场的仿真。

为了解决传统层流冷却技术所存在的冷却能力不足、冷却均匀性差等问题，作为TMCP工艺主要组成部分的轧后冷却技术，已由传统的层流冷却技术向以超快速冷却为核心的新一代轧后冷却控制技术发展。首钢、鞍钢等中厚板生产线的超快冷工业实践表明，超快速冷却具有超常的冷却能力，约为普通层流冷却的2~5倍。超快冷过程中，钢板表面直接与冷却介质接触，由于钢板本身导热热阻的存在，使得钢板内部会产生很大的温度梯度，而单元属性直接与温度相关，因此倘若超快冷过程中仍忽略单元属性对单元节点温度的“贡献”而采用有限差分法进行温度场的仿真，势必会影响模型的计算精度，迫切需要基于有限元法开发出适合于超快冷过程的温度场控制模型，然而采用传统的有限元法所构建的温度场模型在应用过程中时常表现出振荡问题。虽然人们进行了大量卓有成效的工作，如细分网格、时间步长动态变化和集中热容矩阵等，在一定程度上解决了有限元求解温度场过程中的振荡问题，但效果都不太理想。

发明内容

针对现有技术存在的不足，本发明提供一种中厚板轧后超快冷过程温度场的耦合控制方法，有效解决有限元法进行温度场仿真时所产生的振荡问题，实现中厚板超快冷过程的在线温度场解析和冷却工艺的高精度设定。

本发明方法具体按如下步骤进行：

步骤1：在轧机倒数第二道次抛钢时刻，超快冷过程自动化系统接收轧机发送的PDI（primary data input）数据：包括钢板ID、钢板的厚度、钢板的长度、钢板的宽度、钢板的终轧温度、钢板的终冷温度、钢板的化学成分、冷却速度以及钢板的冷却模式（是否需要进行水冷）；

步骤2：对接收到的PDI数据进行校核，判断所接收的PDI数据是否满足该厂的工艺参数限定范围；若满足，执行步骤3；若不满足，则输出错误信息，此时轧钢现场工作人员需根据所产生的错误信息查找产生错误的原因并提出解决方案，超快冷过程自动化系统将当前钢板按不需要水冷处理，执行步骤11；

步骤3：根据PDI数据，利用温度场耦合控制算法对每一个水流密度层别分别进行温度场仿真；

所述温度场耦合控制算法，具体按如下步骤执行：

步骤3.1：建立中厚板轧后超快冷过程的有限元温度场模型；

本发明方法是把求解超快冷过程中钢板温度随时间和空间连续分布的问题转化为在时间领域与空间领域内求解有限个离散点的温度值问题，用这些离散点的温度值去逼近连续的温度分布。因此结合中厚板生产的特点，在空间领域内将钢板的厚向离散化为E个单元，n个节点，在时间域内将水冷时间进行离散化，则离散化后的有限元模型如图1所示，并对单元和节点进行编号；

步骤3.2：确定有限元温度场模型的边界条件；

中厚板超快冷过程中，其水冷换热系数模型为：

α＝k·k_w·k_p·a·Q^b·e^-cT    (1)

式中：k为自适应系数；k_w为水温补偿系数；k_p为压力补偿系数；a，b，c为模型回归系数；Q为超快冷集管的水流密度；T为钢板表面温度；α为水冷换热系数。该水冷换热系数模型即为有限元温度场模型的边界条件；

步骤3.3：构建单元温度插值函数；

对于离散后的每个单元而言，单元内部的温度可认为是线性变化的，因此对任意单元e假设其单元插值函数T^(e)为：

T^(e)=a₁+a₂x

式中：a₁，a₂为常数，可以用节点温度来表示。对于由节点i和节点j组成的单元，设节点i处的温度为T_i，x的值为0；在节点j处的温度为T_j，x的值为L，即单元长度为L，则有：

$T^{\left(e\right)}=\left[\begin{array}{cc}1-\frac{x}{L}& \frac{x}{L}\end{array}\right]\left[\begin{array}{c}{T}_i\\ T_j\end{array}\right]---\left(2\right)$

令：形函数矩阵[N]=[N_i,N_j]，其中

步骤3.4：构建单元比热插值函数；

由于连续分布的温度场已经离散化，形成了由一系列单元组成的温度场，即温度场的自由度已经由无限变成了有限，在对温度场的求解过程中若不对单元的比热进行合理的处理，则会使仿真的温度场结果产生振荡；由于大量钢板的实测比热温度曲线表明，钢板的比热与其温度间呈线性关系；因此，设任意单元e的比热插值函数

为：

$c_p^{\left(e\right)}=c_0+c_{1}T$

式中，c₀，c₁为常数，可用单元节点比热来表示。对于节点i，j组成的单元，设节点i的比热为c_pi，节点j的比热为c_pj，则有下式成立：

$c_p^{\left(e\right)}=(1-\frac{x}{L})c_{\mathrm{pi}}+\frac{x}{L}{c}_{\mathrm{pj}}$

即： $c_p^{\left(e\right)}=N_i{c}_{\mathrm{pi}}+N_j{c}_{\mathrm{pj}}=\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\left[\begin{array}{c}{c}_{\mathrm{pi}}\\ c_{\mathrm{pj}}\end{array}\right]---\left(3\right)$

由比热插值函数的形函数N_i,N_j的表达式可知，形函数有如下性质：

(a)形函数是节点坐标的线性函数，与温度插值函数的类型一致；

(b)在节点处，当形函数的标号与该节点的标号相同时，其值为“1”，否则为“0”，如在节点i处形函数N_i值为“1”，而N_j值为“0”；

(c)在单元内部任意点处，其形函数之和恒为“1”，即N_i+N_j=1；

步骤3.5：构建超快冷过程的单元刚度方程；

以一维无内热源导热微分方程为基础，利用欧拉方程建立等效泛函，对于任意单元e其等效泛函I^e为：

$I^e=\frac{1}{2}{\∫\∫\∫}_{v^e}[k{\left(\frac{\∂T}{\∂x}\right)}^2+2\·{\mathrm{\ρc}}_p\frac{\∂T}{\∂\mathrm{\τ}}T]\mathrm{dv}+\frac{1}{2}\underset{s^e}{\∫\∫}\mathrm{\α}{(T-T_{\∞})}^2\mathrm{ds}---\left(4\right)$

式中：v^e为中厚板超快冷过程所选取的微元体的体积，s^e为该微元体与冷却介质直接接触的面积；ρ,k,c_p,τ分别为钢板的密度、热导率、比热及时间；T_∞为环境温度；

对式(4)求一阶偏导数并置零，则得到如下单元刚度方程：

$\left[k_{1\mathrm{ij}}^e\right]{\left\{T\right\}}^e+\left[k_{3\mathrm{ij}}^e\right]{\left\{\stackrel{\·}{T}\right\}}^e={\left\{p_{1\mathrm{ij}}\right\}}^e---\left(5\right)$

式中：

为单元e的刚度矩阵；{T}^e为单元e的节点温度列阵；

为单元e的变温矩阵；

单元e内各节点温度对时间的偏导数列阵；{p_1ij}^e为单元e的温度载荷列阵；其中：

$\left[k_{1\mathrm{ij}}^e\right]=\underset{v^e}{\∫\∫\∫}{\left[B\right]}^T\left[k\right]\left[B\right]\mathrm{dv}+\underset{s^e}{\∫\∫}\mathrm{\α}\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\mathrm{ds}$

$\left[k_{3\mathrm{ij}}^e\right]=\underset{v^e}{\∫\∫\∫}\mathrm{\ρ}\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\left[\begin{array}{c}{c}_{\mathrm{pi}}\\ c_{\mathrm{pj}}\end{array}\right]\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\mathrm{dv}$

${\left\{p_{1\mathrm{ij}}\right\}}^e=\underset{s^e}{\∫\∫}\mathrm{\α}{T}_{\∞}\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\mathrm{ds}$

式中：[B]为单元的温度应变矩阵， $\left[B\right]=\left[\begin{array}{cc}\frac{{\∂N}_i}{\∂x}& \frac{{\∂N}_j}{\∂x}\end{array}\right];$

现以钢板厚度方向上任意相邻的三个节点为例来说明，如图2所示，设节点i,j组成的单元为①，节点j,k组成的单元为②，节点i,j,k对应的温度分别为T_i,T_j,T_k，节点i,j,k对应的比热分别为：c_pi,c_pj,c_pk，单元①的长度为L₁，单元②的长度为L₂，其中节点i为外部单元，节点j，k均为内部单元，在节点i处，由于节点i与温度为T_∞的冷却介质直接接触发生对流换热，其对流换热系数为α；

对节点j而言，由于单元①和单元②都通过节点j进行热量传递，为了建模的方便，不妨将单元①中的节点j标记为

将单元②中的节点j标记为

对于单元①而言，

根据式(5)建立单元刚度方程：

$\left[B\right]=\left[\begin{array}{cc}-\frac{1}{L_1}& \frac{1}{L_1}\end{array}\right]---\left(6\right)$

$\left[k_{1\mathrm{ij}}^e\right]=\underset{v^e}{\∫\∫\∫}{\left[B\right]}^T\left[k\right]\left[B\right]\mathrm{dv}+\underset{s^e}{\∫\∫}\mathrm{\α}\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\mathrm{ds}=A\·\left[\begin{array}{cc}\frac{k}{L_1}+\mathrm{\α}& -\frac{k}{L_1}\\ -\frac{k}{L_1}& \frac{k}{L_1}\end{array}\right]---\left(7\right)$

$\left[k_{3\mathrm{ij}}^e\right]=\underset{v^e}{\∫\∫\∫}\mathrm{\ρ}\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\left[\begin{array}{c}{c}_{\mathrm{pi}}\\ c_{\mathrm{pj}}\end{array}\right]\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\mathrm{dv}=A\·\left[\begin{array}{cc}\frac{L_1}{12}(3{\mathrm{\ρc}}_{\mathrm{pi}}+{\mathrm{\ρc}}_{\mathrm{pj}})& \frac{L_1}{12}({\mathrm{\ρc}}_{\mathrm{pj}}+{\mathrm{\ρc}}_{\mathrm{pj}})\\ \frac{L_1}{12}({\mathrm{\ρc}}_{\mathrm{pi}}+{\mathrm{\ρc}}_{\mathrm{pj}})& \frac{L_1}{12}({\mathrm{\ρc}}_{\mathrm{pi}}+3{\mathrm{\ρc}}_{\mathrm{pj}})\end{array}\right]---\left(8\right)$

${\left\{p_{1\mathrm{ij}}\right\}}^e=\underset{s^e}{\∫\∫}\mathrm{\α}{T}_{\∞}\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\mathrm{ds}=A\·\left[\begin{array}{c}\mathrm{\α}{T}_{\∞}\\ 0\end{array}\right]---\left(9\right)$

其中A为单位面积

整理得：

同理可得，单元②的刚度方程为：

步骤3.6：构建整体刚度方程；

利用扩大阶数法构建钢板超快冷过程整体的刚度矩阵方程，由于

因此得到如下整体刚度方程：

$\left[K_1\right]\left\{T\right\}+\left[K_3\right]\left\{\stackrel{\·}{T}\right\}=\left\{P\right\}---\left(12\right)$

式中：[K₁]为温度刚度矩阵，

[K₃]为变温矩阵，

{T}为节点温度列阵；{P}为节点载荷列阵，

为节点温度对时间的偏导数列阵， ${\left[\stackrel{\·}{T}\right]}^T=\left[\begin{array}{cccc}{\∂T}_1/\∂\mathrm{\τ}& {\∂T}_2/\∂\mathrm{\τ}& \·\·\·& {\∂T}_n/\∂\mathrm{\τ}\end{array}\right];$ E为单元总数；e表示单元编号为e的单元；T₁,T₂,…,T_n为节点温度；n为节点总数；

在整体温度刚度方程中，变温矩阵[K₃]具有如下的形式：

$\left[K_3\right]={\left[\begin{array}{cccccc}{a}_1& b_1& 0& 0& 0& 0\\ c_1& a_2& b_2& 0& 0& 0\\ 0& c_2& \·& \·& 0& 0\\ 0& 0& \·& \·& b_{n-2}& 0\\ 0& 0& 0& c_{n-2}& a_{n-1}& b_{n-1}\\ 0& 0& 0& 0& c_{n-1}& a_n\end{array}\right]}_{n\×n}---\left(13\right)$

其中矩阵元素：

$b_i=c_i=\frac{L_i}{12}(c_{\mathrm{pi}}+c_{\mathrm{pi}+1}),$

$a_1=\frac{L_1}{12}(3c_{p1}+c_{p2}),$

$a_i=\frac{L_{i-1}}{12}(c_{\mathrm{pi}-1}+{3c}_{\mathrm{pi}})+\frac{L_i}{12}({3c}_{\mathrm{pi}}+c_{\mathrm{pi}+1}),(i=\mathrm{2,3}\·\·\·n-1)$

$a_n=\frac{L_{n-1}}{12}(c_{\mathrm{pn}-1}+{3c}_{\mathrm{pn}})$

式中：n为钢板离散化时的节点总数；L_i为第i个单元的长度，c_pi为第i个节点当前温度下的比热；

从其表达式可以看出，本发明的变温矩阵[K₃]具有如下特征：

(a)[K₃]是一个稀疏阵，整体变温矩阵中绝大多数元素为零，非零元素集中在主对角线附近呈带状分布，其半带宽为2；

(b)[K₃]为三对角对称矩阵，且矩阵的非零元素均为正值；

(c)[K₃]与温度刚度矩阵[K₁]具有相同的形式；

(d)[K₃]表示的物理意义十分明确：对于任意一个节点，其温度的变化只与包含该节点的单元有关，即只有包含该节点的单元才对节点的温度变化有“贡献”，而其他单元对该节点的“贡献”为零；

步骤3.7：时域问题的处理及温度场的求解；

由于上述温度场刚度方程中涉及到时域问题，因此常采用Galerkin差分格式进行处理，则：

$\left(\frac{2}{3}\mathrm{\Δ\τ}\right[K_1]+[K_3\left]\right\{T\}{|}_{\mathrm{\τ}}=\{\left[K_3\right]-\frac{1}{3}\mathrm{\Δ\τ}\left[K_1\right]\left\}\right\{T\}{|}_{\mathrm{\τ}-\mathrm{\Δ\τ}}+\mathrm{\Δ\τ}\{P\}---\left(14\right)$

式中：{T}|_τ为待解的τ时刻的节点温度列阵；{T}|_τ-Δτ为τ-Δτ时刻已知的节点温度列阵；Δτ为时间步长；

为了便于说明，式(14)可进一步简化为如下形式：

A{T}|_τ=b    (15)

其中，A为节点温度的系数矩阵，

b为常数项列阵， $b=\left(\right[K_3]-\frac{1}{3}\mathrm{\Δ\τ}[K_1\left]\right)\left\{T\right\}{|}_{\mathrm{\τ}-\mathrm{\Δ\τ}}+\mathrm{\Δ\τ}\left\{P\right\};$

由于式(15)为对称正定方程组，对A进行Cholesky分解A=GG^T，则将式(15)转换为如下的方程组：

$\left\{\begin{array}{c}\mathrm{Gy}=b\\ G^T\left\{T\right\}{|}_{\mathrm{\τ}}=y\end{array}\right.$

则求得节点温度列阵：{T}|_τ=(G^T)^-1×(G^-1×b)

式中：G为进行Cholesky分解获得的下三角矩阵；G^T为矩阵G的转制；y为中间变量；(G^T)^-1为矩阵G^T的逆矩阵；G^-1为矩阵G的逆矩阵；

步骤3.8：判断循环终止条件；

本发明的有限元温度场模型中可以以终冷温度作为循环终止条件，也可以以钢板冷却时间作为循环终止条件，应根据工艺的需要进行选择；若满足循环终止条件，则结束循环；否则以步骤3.7计算出的温度场作为初始温度场，开始下一次迭代计算，直到满足循环终止条件为止；

步骤4：根据冷却工艺要求，进行冷却规程设定，包括超快冷集管的流量、开启组数及集管组态，并下发给超快冷的基础自动化系统；

利用温度场耦合控制算法的计算结果，可直接获得各水流密度层别对应的水冷时间，进而可以获得该超快冷钢板各水流密度层别条件下对应的冷却速度。此后再根据工艺要求所设定的冷却速度进行线性插值确定达到目标冷却速率所需要的水流密度以及所需要的水冷时间，进而确定出集管的开启数量及集管的组态。所设定的冷却规程将下发给基础自动化系统执行具体操作；

每一个水流密度层别对应的冷却速度C_r由下式计算：

$C_r=\frac{T_{\mathrm{start}}-T_{\mathrm{stop}}}{t}---\left(16\right)$

式中：T_start为开冷温度，T_stop为终冷温度；t为冷却时间，由温度场仿真结果直接获得；

步骤5：在轧机末道次抛钢时刻，超快冷过程自动化系统再次接收轧机发送的PDI数据，重复步骤2-4，进行超快冷钢板冷却工艺的修正；

步骤6：当钢板尾部离开精轧机机后测温仪时，超快冷过程自动化系统接收轧机基础自动化系统所采集的钢板轧后实测温度，触发超快冷过程自动化系统进行钢板超快冷过程的辊道速度设定，并将其下发给超快冷基础自动化系统，再由超快冷基础自动化系统转发给轧机基础自动化系统执行；

步骤7：冷却过程中，超快冷基础自动化系统不断采集钢板的实测温度，包括超快冷设备入口处测温仪的实测温度以及超快冷设备出口处测温仪的实测温度，并将该温度上传给超快冷过程自动化系统；

步骤8：钢板尾部离开开启的最后一组集管时，超快冷基础自动化系统顺序关闭开启的所有集管；

步骤9：当钢板尾部离开超快冷设备出口测温仪时，超快冷过程自动化系统利用步骤3所述温度场耦合控制算法进行模型的自学习；

步骤10：自学习结束后，生成冷却结果报表并上传给轧机过程自动化系统，便于对冷却结果进行查询；

步骤11：超快冷过程自动化系统上传完冷却结果报表后，进入等待阶段，当下一块钢板轧制到倒数第二道次抛钢时，重复步骤1-11。

执行上述控制方法所采用的控制系统，包括超快冷基础自动化系统、超快冷过程自动化系统、超快冷设备、侧喷装置和检测仪表，其中，超快冷基础自动化系统采用1套PLC系统作为超快冷设备的主控单元，通过Profibus DP网与远程I/O模块通讯，远程I/O模块通过电缆和现场仪表直接连接；超快冷过程自动化系统采用一台服务器作为超快冷过程机，通过工业以太网与超快冷基础自动化系统进行数据通讯，用于工艺过程控制模型及模块的计算处理，实现合理的超快冷工艺参数设定；超快冷设备由若干集管组成，在每组集管的供水管路上都安装调节阀和流量计，其中流量计用于检测各集管的水流量，而调节阀则根据所检测的水流量不断调节阀门的开闭程度，以使各集管的实际水流量与过程自动化系统所下发的冷却工艺一致。为了清除钢板上表面残水，在超快冷设备上安装有侧喷装置，在安装侧喷装置的供水管路上安装有开闭阀，用于开启和关闭侧喷装置；为了检验控制精度，检测仪表选用测温仪，实时测量钢板表面温度，测温仪安装在超快冷设备的入口和出口处。

上述的温度场耦合控制方法嵌入超快冷过程自动化系统中；超快冷过程自动化系统根据所采集到的PDI参数自动进行冷却工艺的设定，并下发给超快冷基础自动化系统。超快冷基础自动化系统接收并执行所下发的冷却工艺，同时将自动采集到的轧线上仪表的实测数据上传给超快冷过程自动化系统，触发超快冷过程自动化系统进行冷却工艺的修正和自学习。超快冷基础自动化系统采集轧线上仪表的实测数据包括：调节阀的开闭程度、流量计检测的流量、测温仪检测的钢板表面温度。

有益效果：本发明方法实现了比热与温度场的耦合控制，将单元属性对节点温度的“贡献”添加到了有限元温度场模型中；所构建的变温矩阵具有与刚度矩阵相同的形式，都是半带宽为2的对称矩阵；有效解决了传统有限元法进行超快冷过程温度场仿真时所表现出的振荡问题；该温度场耦合控制方法成功应用于中厚板超快冷过程的在线温度场解析，实现了冷却工艺的高精度设定。

附图说明

图1为本发明方法有限元网格划分模型；

图2为本发明方法一维有限单元模型；

图3为本发明实施例控制方法的流程图；

图4为本发明实施例控制方法中温度场耦合控制算法流程图；

图5为本发明实施例控制方法中设定超快冷钢板冷却规程的流程图；

图6为本发明实施例一本发明方法的仿真结果图；

图7为本发明实施例一传统有限元法的仿真结果图；

图8为本发明实施例二控制系统结构示意图；

图9为本发明实施例二温度场耦合控制方法的在线应用效果。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

本发明的方法流程如图3所示，具体按如下步骤进行：

步骤1：在轧机倒数第二道次抛钢时刻，过程自动化系统接收轧机发送的PDI（primary datainput）数据：包括钢板ID、钢板的厚度、钢板的长度、钢板的宽度、钢板的终轧温度、钢板的终冷温度、钢板的化学成分、冷却速度以及钢板的冷却模式（是否需要进行水冷）；

步骤2：对接收到的PDI数据进行校核，判断所接收的PDI数据是否满足该厂的工艺参数限定范围；若满足，执行步骤3；若不满足，则输出错误信息，此时轧钢现场工作人员需根据所产生的错误信息查找产生错误的原因并提出解决方案，超快冷过程自动化系统将当前钢板按不需要水冷处理，执行步骤11；

步骤3：根据PDI数据，利用温度场耦合控制算法对每一个水流密度层别分别进行温度场仿真；

所述温度场耦合控制算法，流程图如图4所示，具体按如下步骤执行：

步骤3.1：建立中厚板轧后超快冷过程的有限元温度场模型；

本发明方法是把求解超快冷过程中钢板温度随时间和空间连续分布的问题转化为在时间领域与空间领域内求解有限个离散点的温度值问题，用这些离散点的温度值去逼近连续的温度分布。因此结合中厚板生产的特点，在空间领域内将钢板的厚向离散化为E个单元，n个节点，在时间域内将水冷时间进行离散化，则离散化后的有限元模型如图1所示，并对单元和节点进行编号；

步骤3.2：确定有限元温度场模型的边界条件；

中厚板超快冷过程中，其水冷换热系数模型为：

α＝k·k_w·k_p·a·Q^b·e^-cT    (1)

式中：k为自适应系数；k_w为水温补偿系数；k_p为压力补偿系数；a，b，c为模型回归系数；Q为超快冷集管的水流密度；T为钢板表面温度；α为水冷换热系数，水冷换热系数模型即为有限元温度场模型的边界条件；

步骤3.3：构建单元温度插值函数；

对任意单元e假设单元插值函数T^(e)为：

T^(e)=a₁+a₂x

式中：a₁，a₂为常数，可以用节点温度来表示，对于由节点i和节点j组成的单元，设节点i处的温度为T_i，x的值为0；在节点j处的温度为T_j，x的值为L，即单元长度为L，则有：

$T^{\left(e\right)}=\left[\begin{array}{cc}1-\frac{x}{L}& \frac{x}{L}\end{array}\right]\left[\begin{array}{c}{T}_i\\ T_j\end{array}\right]---\left(2\right)$

令：形函数矩阵[N]=[N_i,N_j]，其中

步骤3.4：构建单元比热插值函数；

设任意单元e的比热插值函数

为：

$c_p^{\left(e\right)}=c_0+c_{1}T$

式中，c₀，c₁为常数，可用单元节点比热来表示。对于节点i，j组成的单元，设节点i的比热为c_pi，节点j的比热为c_pj，则有下式成立：

$c_p^{\left(e\right)}=(1-\frac{x}{L})c_{\mathrm{pi}}+\frac{x}{L}{c}_{\mathrm{pj}}$

即： $c_p^{\left(e\right)}=N_i{c}_{\mathrm{pi}}+N_j{c}_{\mathrm{pj}}=\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\left[\begin{array}{c}{c}_{\mathrm{pi}}\\ c_{\mathrm{pj}}\end{array}\right]---\left(3\right)$

由比热插值函数的形函数N_i,N_j的表达式可知，形函数有如下性质：

(a)形函数是节点坐标的线性函数，与温度插值函数的类型一致；

(b)在节点处，当形函数的标号与该节点的标号相同时，其值为“1”，否则为“0”，如在节点i处形函数N_i值为“1”，而N_j值为“0”；

(c)在单元内部任意点处，其形函数之和恒为“1”，即N_i+N_j=1；

步骤3.5：构建超快冷过程的单元刚度方程；

以一维无内热源导热微分方程为基础，利用欧拉方程建立等效泛函，对于任意单元e其等效泛函I^e为：

$I^e=\frac{1}{2}{\∫\∫\∫}_{v^e}[k{\left(\frac{\∂T}{\∂x}\right)}^2+2\·{\mathrm{\ρc}}_p\frac{\∂T}{\∂\mathrm{\τ}}T]\mathrm{dv}+\frac{1}{2}\underset{s^e}{\∫\∫}\mathrm{\α}{(T-T_{\∞})}^2\mathrm{ds}---\left(4\right)$

式中：v^e为中厚板超快冷过程所选取的微元体的体积，s^e为该微元体与冷却介质直接接触的面积；ρ,k,c_p,τ分别为钢板的密度、热导率、比热及时间；T_∞为环境温度；

对式(4)求一阶偏导数并置零，则得到如下单元刚度方程：

$\left[k_{1\mathrm{ij}}^e\right]{\left\{T\right\}}^e+\left[k_{3\mathrm{ij}}^e\right]{\left\{\stackrel{\·}{T}\right\}}^e={\left\{p_{1\mathrm{ij}}\right\}}^e---\left(5\right)$

式中：

为单元e的刚度矩阵；{T}^e为单元e的节点温度列阵；

为单元e的变温矩阵；

单元e内各节点温度对时间的偏导数列阵；{p_1ij}^e为单元e的温度载荷列阵；其中：

$\left[k_{1\mathrm{ij}}^e\right]=\underset{v^e}{\∫\∫\∫}{\left[B\right]}^T\left[k\right]\left[B\right]\mathrm{dv}+\underset{s^e}{\∫\∫}\mathrm{\α}\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\mathrm{ds}$

$\left[k_{3\mathrm{ij}}^e\right]=\underset{v^e}{\∫\∫\∫}\mathrm{\ρ}\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\left[\begin{array}{c}{c}_{\mathrm{pi}}\\ c_{\mathrm{pj}}\end{array}\right]\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\mathrm{dv}$

${\left\{p_{1\mathrm{ij}}\right\}}^e=\underset{s^e}{\∫\∫}\mathrm{\α}{T}_{\∞}\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\mathrm{ds}$

式中：[B]为单元的温度应变矩阵， $\left[B\right]=\left[\begin{array}{cc}\frac{{\∂N}_i}{\∂x}& \frac{{\∂N}_j}{\∂x}\end{array}\right];$

钢板厚度方向上任意相邻的三个节点i,j,k，如图2所示，设节点i,j组成的单元为①，节点j,k组成的单元为②，节点i,j,k对应的温度分别为T_i,T_j,T_k，节点i,j,k对应的比热分别为：c_pi,c_pj,c_pk，单元①的长度为L₁，单元②的长度为L₂，其中节点i为外部单元，节点j，k均为内部单元，在节点i处，由于节点i与温度为T_∞的冷却介质直接接触发生对流换热，其对流换热系数为α；

将单元①中的节点j标记为

将单元②中的节点j标记为

对于单元①而言，

根据式(5)建立单元刚度方程：

$\left[B\right]=\left[\begin{array}{cc}-\frac{1}{L_1}& \frac{1}{L_1}\end{array}\right]---\left(6\right)$

$\left[k_{1\mathrm{ij}}^e\right]=\underset{v^e}{\∫\∫\∫}{\left[B\right]}^T\left[k\right]\left[B\right]\mathrm{dv}+\underset{s^e}{\∫\∫}\mathrm{\α}\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\mathrm{ds}=A\·\left[\begin{array}{cc}\frac{k}{L_1}+\mathrm{\α}& -\frac{k}{L_1}\\ -\frac{k}{L_1}& \frac{k}{L_1}\end{array}\right]---\left(7\right)$

$\left[k_{3\mathrm{ij}}^e\right]=\underset{v^e}{\∫\∫\∫}\mathrm{\ρ}\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\left[\begin{array}{c}{c}_{\mathrm{pi}}\\ c_{\mathrm{pj}}\end{array}\right]\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\left[\begin{array}{cc}{N}_i& N_j\end{array}\right]\mathrm{dv}=A\·\left[\begin{array}{cc}\frac{L_1}{12}(3{\mathrm{\ρc}}_{\mathrm{pi}}+{\mathrm{\ρc}}_{\mathrm{pj}})& \frac{L_1}{12}({\mathrm{\ρc}}_{\mathrm{pi}}+{\mathrm{\ρc}}_{\mathrm{pj}})\\ \frac{L_1}{12}({\mathrm{\ρc}}_{\mathrm{pi}}+{\mathrm{\ρc}}_{\mathrm{pj}})& \frac{L_1}{12}({\mathrm{\ρc}}_{\mathrm{pi}}+3{\mathrm{\ρc}}_{\mathrm{pj}})\end{array}\right]---\left(8\right)$

${\left\{p_{1\mathrm{ij}}\right\}}^e=\underset{s^e}{\∫\∫}\mathrm{\α}{T}_{\∞}\left[\begin{array}{c}{N}_i\\ N_j\end{array}\right]\mathrm{ds}=A\·\left[\begin{array}{c}\mathrm{\α}{T}_{\∞}\\ 0\end{array}\right]---\left(9\right)$

其中A为单位面积

整理得：

同理可得，单元②的刚度方程为：

步骤3.6：构建整体刚度方程；

利用扩大阶数法构建钢板超快冷过程整体的刚度矩阵方程，由于

因此得到如下整体刚度方程：

$\left[K_1\right]\left\{T\right\}+\left[K_3\right]\left\{\stackrel{\·}{T}\right\}=\left\{P\right\}---\left(12\right)$

式中：[K₁]为温度刚度矩阵，

[K₃]为变温矩阵，

{T}为节点温度列阵；{P}为节点载荷列阵，

为节点温度对时间的偏导数列阵， ${\left[\stackrel{\·}{T}\right]}^T=\left[\begin{array}{cccc}{\∂T}_1/\∂\mathrm{\τ}& {\∂T}_2/\∂\mathrm{\τ}& \·\·\·& {\∂T}_n/\∂\mathrm{\τ}\end{array}\right];$ E为单元总数；e表示单元编号为e的单元；T₁,T₂,…,T_n为节点温度；n为节点总数；

在整体温度刚度方程中，变温矩阵[K₃]具有如下的形式：

$\left[K_3\right]={\left[\begin{array}{cccccc}{a}_1& b_1& 0& 0& 0& 0\\ c_1& a_2& b_2& 0& 0& 0\\ 0& c_2& \·& \·& 0& 0\\ 0& 0& \·& \·& b_{n-2}& 0\\ 0& 0& 0& c_{n-2}& a_{n-1}& b_{n-1}\\ 0& 0& 0& 0& c_{n-1}& a_n\end{array}\right]}_{n\×n}---\left(13\right)$

其中矩阵元素：

$b_i=c_i=\frac{L_i}{12}(c_{\mathrm{pi}}+c_{\mathrm{pi}+1}),$

$a_1=\frac{L_1}{12}(3c_{p1}+c_{p2}),$

$a_i=\frac{L_{i-1}}{12}(c_{\mathrm{pi}-1}+{3c}_{\mathrm{pi}})+\frac{L_i}{12}({3c}_{\mathrm{pi}}+c_{\mathrm{pi}+1}),(i=\mathrm{2,3}\·\·\·n-1)$

$a_n=\frac{L_{n-1}}{12}(c_{\mathrm{pn}-1}+{3c}_{\mathrm{pn}})$

式中：n为钢板离散化时的节点总数；L_i为第i个单元的长度，c_pi为第i个节点当前温度下的比热；

从其表达式可以看出，本发明的变温矩阵[K₃]具有如下特征：

(a)[K₃]是一个稀疏阵，整体变温矩阵中绝大多数元素为零，非零元素集中在主对角线附近呈带状分布，其半带宽为2；

(b)[K₃]为三对角对称矩阵，且矩阵的非零元素均为正值；

(c)[K₃]与温度刚度矩阵[K₁]具有相同的形式；

(d)[K₃]表示的物理意义十分明确：对于任意一个节点，其温度的变化只与包含该节点的单元有关，即只有包含该节点的单元才对节点的温度变化有“贡献”，而其他单元对该节点的“贡献”为零；

步骤3.7：时域问题的处理及温度场的求解；

采用Galerkin差分格式进行处理温度场刚度方程中所涉及到的时域问题，则：

$\left(\frac{2}{3}\mathrm{\Δ\τ}\right[K_1]+[K_3\left]\right\{T\}{|}_{\mathrm{\τ}}=\{\left[K_3\right]-\frac{1}{3}\mathrm{\Δ\τ}\left[K_1\right]\left\}\right\{T\}{|}_{\mathrm{\τ}-\mathrm{\Δ\τ}}+\mathrm{\Δ\τ}\{P\}---\left(14\right)$

式中：{T}|_τ为待解的τ时刻的节点温度列阵；{T}|_τ-Δτ为τ-Δτ时刻已知的节点温度列阵；Δτ为时间步长；

式(14)可进一步简化为如下形式：

A{T}|_τ=b                (15)

其中，A为节点温度的系数矩阵，

b为常数项列阵， $b=\left(\right[K_3]-\frac{1}{3}\mathrm{\Δ\τ}[K_1\left]\right)\left\{T\right\}{|}_{\mathrm{\τ}-\mathrm{\Δ\τ}}+\mathrm{\Δ\τ}\left\{P\right\};$

由于式(15)为对称正定方程组，对A进行Cholesky分解A=GG^T，则将式(15)转换为如下的方程组：

$\left\{\begin{array}{c}\mathrm{Gy}=b\\ G^T\left\{T\right\}{|}_{\mathrm{\τ}}=y\end{array}\right.$

则求得节点温度列阵：{T}|_τ=(G^T)^-1×(G^-1×b)

式中：G为进行Cholesky分解获得的下三角矩阵；G^T为矩阵G的转制；y为中间变量；(G^T)^-1为矩阵G^T的逆矩阵；G^-1为矩阵G的逆矩阵；

步骤3.8：判断循环终止条件；

本发明的有限元温度场模型中可以以终冷温度作为循环终止条件，也钢板冷却时间作为循环终止条件，应根据工艺的需要进行选择；若满足循环终止条件，则结束循环；否则以步骤3.7计算出的温度场作为初始温度场，开始下一次迭代计算，直到满足循环终止条件为止；

步骤4：根据冷却工艺要求，进行冷却规程设定，包括超快冷集管的流量、开启组数及集管组态，并下发给超快冷的基础自动化系统；

所述的冷却规程的设定，其流程图如图5所示，首先利用温度场耦合控制算法的计算结果，获得各水流密度层别对应的水冷时间，进而可以获得该超快冷钢板各水流密度层别条件下对应的冷却速度。此后再根据工艺要求所设定的冷却速度进行线性插值确定达到目标冷却速率所需要的水流密度以及所需要的水冷时间，进而确定出集管的开启数量及集管的组态。所设定的冷却规程将下发给基础自动化系统执行具体操作；

每一个水流密度层别对应的冷却速度C_r由下式计算：

$C_r=\frac{T_{\mathrm{start}}-T_{\mathrm{stop}}}{t}---\left(16\right)$

式中：T_start为开冷温度，T_stop为终冷温度；t为冷却时间，由温度场仿真结果直接获得；

步骤5：在轧机末道次抛钢时刻，超快冷过程自动化系统再次接收轧机发送的PDI数据，重复步骤2-4，进行超快冷钢板冷却工艺的修正；

步骤6：当钢板尾部离开精轧机机后测温仪时，超快冷过程自动化系统接收轧机基础自动化系统所采集的钢板轧后实测温度，触发超快冷过程自动化系统进行钢板超快冷过程的辊道速度设定，并将其下发给超快冷基础自动化系统，再由超快冷基础自动化系统转发给轧机基础自动化系统执行；

步骤7：冷却过程中，超快冷基础自动化系统不断采集钢板的实测温度，包括超快冷设备入口处测温仪的实测温度以及超快冷设备出口处测温仪的实测温度，并将该温度上传给超快冷过程自动化系统；

步骤8：钢板尾部离开开启的最后一组集管时，超快冷基础自动化系统顺序关闭开启的所有集管；

步骤9：当钢板尾部离开超快冷设备出口测温仪时，超快冷过程自动化系统利用步骤3所述温度场耦合控制算法进行模型的自学习；

步骤10：自学习结束后，生成冷却结果报表并上传给轧机过程自动化系统，便于对冷却结果进行查询；

步骤11：超快冷过程自动化系统上传完冷却结果报表后，进入等待阶段，当下一块钢板轧制到倒数第二道次抛钢时，重复步骤1-11。

执行上述控制方法所采用的控制系统，包括超快冷基础自动化系统、超快冷过程自动化系统、超快冷设备、侧喷装置和检测仪表，其中，超快冷基础自动化系统采用1套PLC系统作为超快冷设备的主控单元，通过Profibus DP网与远程I/O模块通讯，远程I/O模块通过电缆和现场仪表直接连接；超快冷过程自动化系统采用一台服务器作为超快冷过程机，通过工业以太网与超快冷基础自动化系统进行数据通讯，用于工艺过程控制模型及模块的计算处理，实现合理的超快冷工艺参数设定；超快冷设备由若干集管组成，在每组集管的供水管路上都安装调节阀和流量计，其中流量计用于检测各集管的水流量，而调节阀则根据所检测的水流量不断调节阀门的开闭程度，以使各集管的实际水流量与超快冷基础自动化系统所下发的冷却工艺一致。为了清除钢板上表面残水，在超快冷设备上安装有侧喷装置，在安装侧喷装置的供水管路上安装有开闭阀，用于开启和关闭侧喷装置；为了检验控制精度，选用测温仪实时测量钢板表面温度，测温仪安装在超快冷设备的入口和出口处。

上述的温度场耦合控制方法嵌入超快冷过程自动化系统中；超快冷过程自动化系统根据所接收到的PDI参数自动进行冷却工艺的设定，并下发给超快冷基础自动化系统。超快冷基础自动化系统接收并执行所下发的冷却工艺，同时将自动采集到的轧线上仪表的实测数据上传给超快冷过程自动化系统，触发超快冷过程自动化系统进行冷却工艺的修正和自学习。超快冷基础自动化系统采集轧线上仪表的实测数据包括：调节阀的开闭程度、流量计检测的流量、测温仪检测的钢板表面温度。

实施例一：

选用钢板ID号为038682000，钢种Q550D，钢板厚度：50mm；开冷温度780℃，水冷时间21S；化学成分：C:0.084%；Mn:1.8%；Cr：0.02%；Nb:0.045%；V：0.051%；Mo：0.002%；Ni：0.01%；Ti：0.016%；Cu：0.01%；水温24℃；密度：7800kg/m3。比热：根据其化学成分从表1中查表获得；热导率：根据其化学成分从表2中查表获得；水冷换热系数取为水冷过程的综合水冷换热系数4000W/(m²·K)；时间步长Δτ=3S；将钢板的厚向按等间距划分为20个有限单元，21个节点；

表1比热表(kJ/kg·K)

表2热传导系数表（W/m·K）

基于上述仿真条件，利用本发明的温度场耦合控制方法和传统的有限元法所进行的温度场仿真结果如表3所示。

表3温度场数值仿真结果的对比

为了便于比较，将其绘制成曲线，其中图6为利用本发明的耦合控制方法所进行的温度场仿真结果，图7为传统有限元法所进行的温度场仿真结果。从图7中可以看出，对于钢板的表面节点而言，采用传统的有限元法进行温度场仿真时出现了振荡现象，具体表现为第二次迭代计算的温度为381.5℃较第一次迭代计算的温度362.2℃要高，这显然违背了热力学的基本规律。而采用本发明的温度场耦合控制方法则很好的解决了上述问题。

实施例二：

本实施例控制方法所采用的控制系统，包括超快冷基础自动化系统、超快冷过程自动化系统、超快冷设备、侧喷装置和检测仪表，其系统结构示意图如图8所示。

本实施例超快冷基础自动化系统采用一套西门子S7-400PLC作为超快冷设备的主控单元，通过Profibus DP网与远程I/O模块通讯，远程I/O模块通过电缆和现场仪表直接连接；其硬件包括：UR1机架一个，电源PS 40720A一个，中央处理器CPU 416-2DP一个，通讯模块CP 443-1一个，计数器模块FM 451FIX SPEED两个，以及AI、AO模块、DI、DO模块若干；其应用软件：采用SIMATIC STEP7应用软件。

本实施例超快冷过程自动化系统选择HP DX2818MT(250W)服务器作为超快冷过程机，其系统软件选用Microsoft Windows Server 2003 Enterprise Edition Service Pack 1，应用软件选用Microsoft SQL Server 2005、Visual C++6.0和Microsoft Visual Studio.NET 2003。该服务器的外部输入设备为键盘终端，用来输入程序、数据和操作指令，输出设备为ViewSonic液晶显示器，用以显示超快冷模型的计算结果。

本实施例超快冷设备有9组集管，在每组集管的供水管路上均安装有流量计和调节阀，此外，在超快冷设备上安装有侧喷装置，在其供水管路上安装有开闭阀，用于开启和关闭侧喷装置；在超快冷设备的入口和出口处安装了测温仪，用于实时测量钢板的表面温度。

本实施例超快冷过程自动化系统根据冷却工艺要求进行冷却规程的设定和辊道速度制度的设定。其中冷却规程的设定包括集管流量、集管开组组数、集管排布方式等；待超快冷过程自动化系统完成冷却规程的设定后，通过工业以太网将该冷却规程下发给超快冷基础自动化系统。超快冷基础自动化系统接收并执行超快冷过程自动化系统所下发的冷却规程。待超快冷过程自动化系统设定完辊道速度制度后，通过工业以太网将其下发给超快冷基础自动化系统，再由超快冷基础自动化系统通过Profibus DP网将该辊道速度制度转发给轧机基础自动化系统执行具体操作。

现以ID号为0392974000的钢板为例来说明本发明的温度场耦合控制方法的在线应用效果，超快冷过程自动化系统在接收到如表4所示的PDI数据后，首先进行PDI数据的校核，所接收的PDI数据须满足表5所示的工艺参数限定范围。超快冷过程自动化系统基于上述PDI参数，自动计算出如表6所示的冷却规程。

表4.UFC生产用板规格

表5.工艺参数限定范围

表6.UFC用板冷却规程

按照上述冷却工艺，超快冷过程自动化系统利用本发明的温度场耦合控制方法获得钢板轧后冷却过程中的温度场演变规律及该钢板的实际冷却效果如图9所示。由图9可知，本发明的系统及方法具有较高的控制精度，其95%实测终冷温度距离目标终冷温度的偏差在±20℃之内。

Claims (2)

1.一种中厚板轧后超快冷过程温度场的耦合控制方法，其特征在于：具体按如下步骤执行： 

步骤1：在轧机倒数第二道次抛钢时刻，超快冷过程自动化系统接收轧机发送的PDI数据：包括钢板ID、钢板的厚度、钢板的长度、钢板的宽度、钢板的终轧温度、钢板的终冷温度、钢板的化学成分、冷却速度以及钢板的冷却模式； 

步骤2：对接收到的PDI数据进行校核，判断所接收的PDI数据是否满足该厂的工艺参数限定范围；若满足，执行步骤3；若不满足，则输出错误信息，此时轧钢现场工作人员需根据所产生的错误信息查找产生错误的原因并提出解决方案，超快冷过程自动化系统将当前钢板按不需要水冷处理，执行步骤11； 

步骤3：根据PDI数据，利用温度场耦合控制算法对每一个水流密度层别分别进行温度场仿真； 

步骤4：根据冷却工艺要求，进行冷却规程设定，包括超快冷集管的流量、开启组数及集管组态，并下发给超快冷的基础自动化系统； 

每一个水流密度层别对应的冷却速度C_r由下式计算： 

式中：T_start为开冷温度，T_stop为终冷温度；t为冷却时间，由温度场仿真结果直接获得； 

步骤5：在轧机末道次抛钢时刻，超快冷过程自动化系统再次接收轧机发送的PDI数据，重复步骤2-4，进行超快冷钢板冷却工艺的修正； 

步骤6：当钢板尾部离开精轧机机后测温仪时，超快冷过程自动化系统接收轧机基础自动化系统所采集的钢板轧后实测温度，触发超快冷过程自动化系统进行钢板超快冷过程的辊道速度设定，并将其下发给超快冷基础自动化系统，再由超快冷基础自动化系统转发给轧机基础自动化系统执行； 

步骤7：冷却过程中，超快冷基础自动化系统不断采集钢板的实测温度，包括超快冷设备入口处测温仪的实测温度以及超快冷设备出口处测温仪的实测温度，并将该温度上传给超快冷过程自动化系统； 

步骤8：钢板尾部离开开启的最后一组集管时，超快冷基础自动化系统顺序关闭开启的所有集管； 

步骤9：当钢板尾部离开超快冷设备出口测温仪时，超快冷过程自动化系统利用步骤3所述温度场耦合控制算法进行模型的自学习； 

步骤10：自学习结束后，生成冷却结果报表并上传给轧机过程自动化系统，便于对冷却结果进行查询； 

步骤11：超快冷过程自动化系统上传完冷却结果报表后，进入等待阶段，当下一块钢板轧制到倒数第二道次抛钢时，重复步骤1-11。 

2.根据权利要求1所述的中厚板轧后超快冷过程温度场的耦合控制方法，其特征在于：步骤3中所述的温度场耦合控制算法，具体按如下步骤执行： 

步骤3.1：建立中厚板轧后超快冷过程的有限元温度场模型； 

在空间领域内将钢板的厚向离散化为E个单元，n个节点，在时间域内将水冷时间进行离散化，并对单元和节点进行编号； 

步骤3.2：确定有限元温度场模型的边界条件； 

中厚板超快冷过程中，其水冷换热系数模型为： 

α＝k·k_w·k_p·a·Q^b·e^-cT   (1) 

式中：k为自适应系数；k_w为水温补偿系数；k_p为压力补偿系数；a，b，c为模型回归系数；Q为超快冷集管的水流密度；T为钢板表面温度；α为水冷换热系数，水冷换热系数模型即为有限元温度场模型的边界条件； 

步骤3.3：构建单元温度插值函数； 

对任意单元e假设其单元插值函数T^(e)为： 

T^(e)=a₁+a₂x 

式中：a₁，a₂为常数，可以用节点温度来表示；对于由节点i和节点j组成的单元，设节点i处的温度为T_i，单元长度方向的坐标x的值为0；在节点j处的温度为T_j，x的值为L，即在节点j处，x的值等于单元长度L，则有： 

令：形函数矩阵[N]=[N_i,N_j]，其中

步骤3.4：构建单元比热插值函数； 

设任意单元e的比热插值函数

为： 

式中，c₀，c₁为常数，可用单元节点比热来表示；对于节点i，j组成的单元，设节点i 的比热为c_pi，节点j的比热为c_pj，则有下式成立： 

即： 

由比热插值函数的形函数N_i,N_j的表达式可知，形函数有如下性质： 

(a)形函数是节点坐标的线性函数，与温度插值函数的类型一致； 

(b)在节点处，当形函数的标号与该节点的标号相同时，其值为“1”，否则为“0”，如在节点i处形函数N_i值为“1”，而N_j值为“0”； 

(c)在单元内部任意点处，其形函数之和恒为“1”，即N_i+N_j=1； 

步骤3.5：构建超快冷过程的单元刚度方程； 

以一维无内热源导热微分方程为基础，利用欧拉方程建立等效泛函，对于任意单元e其等效泛函I^e为： 

式中：v^e为中厚板超快冷过程所选取的微元体的体积，s^e为该微元体与冷却介质直接接触的面积；ρ,k,c_p,τ分别为钢板的密度、热导率、比热及时间；T_∞为环境温度； 

对式(4)求一阶偏导数并置零，则得到如下单元刚度方程： 

式中：

为单元e的刚度矩阵；{T}^e为单元e的节点温度列阵；

为单元e的变温矩阵；单元e内各节点温度对时间的偏导数列阵；{p_1ij}^e为单元e的温度载荷列阵；其中： 

式中：[B]为单元的温度应变矩阵，

钢板厚度方向上任意相邻的三个节点i、j、k，设节点i,j组成的单元为①，节点j,k组成的单元为②，节点i,j,k对应的温度分别为T_i,T_j,T_k，节点i,j,k对应的比热分别为：c_pi,c_pj,c_pk，单元①的长度为L₁，单元②的长度为L₂；其中节点i为外部单元，节点j，k均为内部单元，在节点i处，由于节点i与温度为T_∞的冷却介质直接接触发生对流换热，其对流换热系数为α； 

将单元①中的节点j标记为

将单元②中的节点j标记为对于单元①，根据式（5）建立单元刚度方程： 

其中A为单位面积 

整理得： 

同理可得，单元②的刚度方程为： 

步骤3.6：构建整体刚度方程； 

利用扩大阶数法构建钢板超快冷过程整体的刚度方程，由于因此得到如下整体刚度方程： 

式中：[K₁]为温度刚度矩阵，

[K₃]为变温矩阵，

{T}为节点温度列阵；{P}为节点载荷列阵，

为节点温度对时间的偏导数列阵， 

E为单元总数；e表示单元编号为e的单元；T₁,T₂,…,T_n为节点温度；n为节点总数； 

在整体温度刚度方程中，变温矩阵[K₃]具有如下的形式： 

其中矩阵元素： 

式中：n为钢板离散化时的节点总数；L_i为第i个单元的长度，c_pi为第i个节点当前温度下的比热； 

从其表达式可以看出，变温矩阵[K₃]具有如下特征： 

(a)[K₃]是一个稀疏阵，整体变温矩阵中绝大多数元素为零，非零元素集中在主对角线附近呈带状分布，其半带宽为2； 

(b)[K₃]为三对角对称矩阵，且矩阵的非零元素均为正值； 

(c)[K₃]与温度刚度矩阵[K₁]具有相同的形式； 

(d)[K₃]表示的物理意义是：对于任意一个节点，其温度的变化只与包含该节点的单元有关，即只有包含该节点的单元才对节点的温度变化有“贡献”，而其他单元对该节点的“贡献”为零； 

步骤3.7：时域问题的处理及温度场的求解； 

采用Galerkin差分格式进行处理温度场刚度方程中所涉及到的时域问题，则： 

式中：{T}|_τ为待解的τ时刻的节点温度列阵；{T}|_τ-Δτ为τ-Δτ时刻已知的节点温度列阵；Δτ为时间步长； 

式(14)进一步简化为如下形式： 

A{T}|_τ=b   (15) 

其中，A为节点温度的系数矩阵，

b为常数项列阵， 

由于式(15)为对称正定方程组，对A进行Cholesky分解A=GG^T，则将式(15)转换为如下的方程组： 

则求得节点温度列阵：{T}|_τ=(G^T)^-1×(G^-1×b) 

式中：G为进行Cholesky分解获得的下三角矩阵；G^T为矩阵G的转置；y为中间变量；(G^T)^-1为矩阵G^T的逆矩阵；G^-1为矩阵G的逆矩阵； 

步骤3.8：判断循环终止条件； 

根据需要选择以终冷温度或钢板冷却时间作为循环终止条件，若满足循环终止条件，则结束循环；否则以步骤3.7计算出的温度场作为初始温度场，开始下一次迭代计算，直到满足循环终止条件为止。 

Images