板带热轧过程中S型变步长法预测瞬态温度场方法

摘要

板带热轧过程中S型变步长法预测瞬态温度场方法，S型变步长的时间步数学模型为：Δt＝a－b×exp(－c×td)，利用S型变步长方法进行瞬态温度场有限元求解不仅提高了计算精度，满足了初始温度场设定要求，而且大大缩短了计算时间。这种方法在有限元法求解瞬态温度场时解决了计算精度低、速度慢的问题，提高了计算效率。

说明书

技术领域

本发明属于轧制技术领域，特别涉及一种板带热轧过程S型变步长法预测瞬态温度场方法。

背景技术

在板带热轧制过程中，从板坯出炉到层流冷却结束，轧件瞬态温度场分布计算对轧制过程轧件组织演变、力学性能和板形等有重要影响。过去经常采用现场测试和数学模型进行确定，由于受实际条件和模型限制影响，降低了温度预测精度和温度分布信息。有限元法作为一种有效的数值分析方法已成为预测板板坯空冷过程瞬态温度场的重要手段。目前在有限元法预测板带热轧过程的瞬态温度场中，往往采用较多单元划分或者较小的定时间步长来提高计算精度，获得更为详尽的温度分布信息。然而单元划分过多和时间步长过小增加了计算时间，降低了计算效率。

发明内容

为了克服上述方法计算速度慢，计算时间和计算精度不能同时满足的缺点，本发明提出一种S型变步长法预测瞬态温度场方法，其目的是保证板带温度预测精度的情况下，缩短计算时间，提高计算效率。

实现本发明目的技术解决方案如下：

1、采集轧制过程数据，包括：轧制参数，材料热物性参数，单元划分信息。

轧制参数：初始时间，轧制时间，轧件宽度，轧件厚度，初始温度，轧件周围介质温度材料热物性参数：热传导系数，黑度，比热，密度

单元划分信息：宽度单元数和厚度单元数

2、根据单元划分数据、轧件宽度和厚度尺寸建立有限元分析模型(见图1)，然后进行单元节点编号、确定换热边界和计算节点坐标。

3、根据轧制过程实际条件和轧制阶段确定换热系数h和内热源强度。热轧过程包括空冷阶段，除磷阶段，轧制阶段。整个轧制过程不同阶段的换热系数和内热源强度计算如下：

(1)热轧板带在空冷过程中，值为零；其表面换热方式主要为辐射和自然对流，换热系数h通过式(1)和式(2)计算：

HR＝σ·ε·(T+T_air)(T²+T_air²)    (1)

其中：HR为辐射系数，σ为Stefan-Boltzman常数，σ＝5.67×10^-8W/(m²·K⁴)；ε为黑度系数，ε与温度的关系式为ε＝0.125(T/1000)²-0.38(T/1000)+1.1。

热轧板带在空冷过程的对流方式为自然对流，其表达式为：

其中：T(K)为板带表面温度；T_air为环境温度；b表示板带宽度。

(2)热轧板带在高压水除磷过程中，值为零；表面换热方式为强迫对流和侧面辐射，换热系数h通过式(1)和式(3)计算：

HC_W＝124.7×w^0.663×10^{-0.00147(T-273.16)}    (3)

其中w(L/min·m²)为水流密度；T(K)板带表面温度。

(3)在轧制过程中，板带与轧辊之间接触换热是主要热损失方式，忽略塑性变形和摩擦做功，$\stackrel{·}{q}=0.$接触换热系数与轧制压力有关。换热系数h通过式(4)计算：

IHTC＝695p_m-34400(W/m²K)    (4)

式中：p_m(MPa)-轧制压力

4、利用有限元基本原理，计算四边形等参单元的形函数N、B矩阵和雅克比矩阵J。

5、以二维热传导基本方程为基础，利用欧拉方程建立等效泛函，确定温度场求解的系统方程。

以热力学第一定律为依据建立无内热源强度的二维热传导的微分方程为：

$k(\frac{{\partial }^{2}T}{{\partial x}^2}+\frac{{\partial }^{2}T}{{\partial y}^2})-\mathrm{\rho c}\frac{\partial T}{\partial t}=0---\left(5\right)$

其中：

T-瞬时温度(K)

ρ-材料密度(kg/m³)

c-材料比热(J/(kg·K))

t-时间(s)

k-热传导系数(W/(m·K))

利用欧拉方程将二维热传导问题方程(5)变为等效泛涵：

$I=\frac{1}{2}{\int \int }_S\left[k\right[{\left(\frac{\partial T}{\partial x}\right)}^2+{\left(\frac{\partial T}{\partial y}\right)}^2]+2\mathrm{\rho c}\frac{\partial T}{\partial t}T]\mathrm{dS}+\frac{1}{2}{\int }_{l}h(T-T_{\infty })\mathrm{dl}---\left(6\right)$

根据热传导问题的变分原理，对泛函式(6)求一阶偏导数并置零，得到温度求解的系统方程：

$\left[K_T\right]\left\{T\right\}+\left[K_3\right]\left\{\frac{\partial T}{\partial t}\right\}=\left\{p\right\}---\left(7\right)$

式中：[K_T]-温度刚度矩阵，$\left[K_T\right]=\underset{e=1}{\overset{E}{\Sigma }}\left(\right[K_1^{\left(e\right)}]+[K_2^{\left(e\right)}\left]\right);$[K₃]-变温矩阵，$\left[K_3\right]=\underset{e=1}{\overset{E}{\Sigma }}\left[K_3^{\left(e\right)}\right];${p}-常数项列式，$\left\{p\right\}=\underset{e=1}{\overset{E}{\Sigma }}\left\{p^{\left(e\right)}\right\};${T}-温度列式；E-单元总数；上标e表示每个单元。

$K_{\mathrm{Tij}}={\int \int }_{S}k(\frac{{\partial N}_i}{\partial x}·\frac{{\partial N}_j}{\partial x}+\frac{{\partial N}_i}{\partial y}·\frac{{\partial N}_j}{\partial y})\mathrm{dS}+{\int }_{L}h{N}_i{N}_j\mathrm{dL}$

K_3ij＝∫∫_SρcN_iN_jdS

p_i＝∫_LhT_∞N_idL

对每个单元来说，刚度矩阵、变温矩阵和常数项可以通过下式求解：

$K_{1\mathrm{ij}}^{\left(e\right)}={\int \int }_{S^e}k(\frac{\partial N_i}{\partial x}·\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}·\frac{\partial N_j}{\partial y})\mathrm{dS}$

$K_{2\mathrm{ij}}^{\left(e\right)}={\int }_{L^e}h{N}_i{N}_j\mathrm{dL}$

$K_{3\mathrm{ij}}^{\left(e\right)}={\int \int }_{S^e}\mathrm{\rho c}{N}_i{N}_j\mathrm{dS}$

${\left\{p_i\right\}}^{\left(e\right)}={\int \int }_{S^e}\stackrel{·}{q}{N}_i\mathrm{dS}++{\int }_{L^e}h{T}_{\infty }{N}_i\mathrm{dL}$

其中：

k-热传导系数(W/(m·K))；

ρ-材料密度(kg/m³)；

c-材料比热(J/(kg·K))；

h-换热系数；

N-形函数；

i，j节点编号。

6、利用二点向后差分格式，将系统方程转化为瞬态温度场求解的线性方程组。

将系统方程(7)中的温度对时间偏导数表示为二点向后差分格式：

$\frac{\partial T}{\partial t}=\frac{1}{\mathrm{\Delta t}}(T_t-T_{t-\mathrm{\Delta t}})---\left(8\right)$

将时间向后差分格式(8)带入系统方程式(7)得到温度场求解的线性方程组：

$\left(\right[K_T]+\frac{1}{\mathrm{\Delta t}}[K_3\left]\right){\left\{T\right\}}_t=\frac{1}{\mathrm{\Delta t}}\left[K_3\right]{\left\{T\right\}}_{t-\mathrm{\Delta t}}+\left\{p\right\}---\left(9\right)$

7、利用S型变步长法设定线性方程组(9)中的时间步长Δt，S型变时间步长的基本原理和数学模型。得到最终的温度场求解的线性方程组：

$\left(\right[K_T]+\frac{1}{(a-b\times \exp (-c\times t^d))}[K_3\left]\right){\left\{T\right\}}_t=\frac{1}{(a-b\times \exp (-c\times t^d))}\left[K_3\right]{\left\{T\right\}}_{t-\mathrm{\Delta t}}+\left\{p\right\}---\left(10\right)$

8、采用一维变带宽存储法求解线性方程组(10)可以获得板带轧制过程瞬态温度场。

9、计算时间在初始时间基础上加上时间步长，根据计算时间和设定轧制时间判断程序计算是否结束，如果计算时间大于轧制时间退出，否则继续返回计算。根据轧制时间判断求解过程是否结束，并将S型变步长法应用到瞬态温度场的有限元求解过程。

其中S型变步长原理与模型如下：

S型变步长基本思想是在温度场计算初期采用较小的时间步长，然后时间步长虽时间增加逐渐增加，最后成为一恒定值，变时间步长的S型曲线见图3所示。由于热分析初期，边界条件和物体本身温差较大，温度变化较为剧烈，因此采用较小的迭代时间步长，得到瞬时温度变化信息，提高计算精度。随着热传导过程的进行，整个系统逐渐趋于热平衡，采用小的时间步长会导致计算时间急剧增加，因此增大时间步长，以缩短计算时间。时间步长不能一味的增加，时间步长过大，不仅会降低温度求解精度，而且一定程度上不能满足瞬时温度场求解的要求。S型变步长不仅能够保证计算精度，而且能够极大的缩短计算时间。

为了更方便的实现S型变步长在温度场有限元计算中的作用，对S型变步长曲线进行了回归，得出如下模型：

Δt＝a-b×exp(-c×t^d)    (11)

其中：t为时间，a和b主要控制迭代时间步长的初始值和终值的大小，k0＝a-b，k1∝a，c和d主要控制S型曲线的变化率。

本发明在传统算法的基础上通过采用S型变步长法预测瞬态温度场，在不影响计算效率的情况下有效地克服了振荡现象，保证了计算的稳定性。

附图说明

图1本发明方法有限元分析模型图，

图2本发明方法具体操作软件流程图，

图3本发明方法S型变步长曲线图，

图4实施例中第一组方案采用本发明方法与采用定步长计算结果结果比较图，

图5实施例中第二组方案采用本发明方法与采用定步长计算结果结果比较图，

图6实施例中第三组方案采用本发明方法与采用定步长计算结果结果比较图，

图7实施例中第一组方案采用本发明方法与采用定步长计算结果结果比较图。

图中：i为单元编号，j为节点编号，H为厚度，W为宽度，1为换热边界，a为主要控制迭代时间步长的初始值的大小，b为主要控制迭代时间步长的终值的大小，k0＝a-b，k1∝a，A-定步长平均温度，B-定步长表面温度，C-定步长中心温度，D-变步长平均温度，E-变步长表面温度，F-变步长中心温度。

具体实施方式

选择一热轧板坯从出炉到除磷前的空冷过程温度变化进行预测分析，比较定步长法和变步长法对板坯温度计算结果和计算时间的影响。

例：计算条件见表1

表1

热物性参数    数值截面尺寸(宽×厚)热传导率比热密度黑度波尔滋曼常数初始温度环境温度冷却时间    1800mm×220mm    k＝30W/(m·K)    c＝620J/(kg·K)    ρ＝7800kg/m³    0.8    5.67×10^-8    1200℃    30℃    300s

计算分析了四组方案，每组方案的单元划分、定时间步长和变步长分别为：

方案1单元划分：20×10；定时间步长：Δt＝1s；变步长设定：a＝10.0，b＝9.0，c＝0.0001，d＝2.5

方案2单元划分：30×10，定时间步长Δt＝1s；变步长设定：a＝10.0，b＝9.0，c＝0.001，d＝2.5

方案3单元划分为40×20，定时间步长Δt＝0.5s；变步长设定：a＝15.0，b＝14.5，c＝0.0001，d＝2.5

方案4单元划分为50×20，定时间步长Δt＝0.1s；变步长设定：a＝15.0，b＝14.9，c＝0.0001，d＝2.5

图4，5，6，7所示为各个方案下定步长和变步长温度变化比较，可以看出在整个冷却时间内，定步长温度计算结果和变步长温度计算结果基本相同，只有心部温度有一定差别，但是相对误差不超过0.2％，所以变步长和定步长温度计算精度相同。通过表2可以看出，变时间步长的计算时间远远小于定时间步长，这种计算速度的优越性随着单元数目的增加更加明显，变时间步长在保证计算精度的情况下，极大的缩短了计算时间，提高了计算效率。

表2

  方案一方案二方案三方案四  t1(ms)  t2(ms)t1(ms)t2(ms)t1(ms)t2(ms)t1(ms)t2(ms)  287  40564716350490482901475

本发明方法针对板坯出炉空冷过程温度变化进行了求解，变步长计算结果和定步长温度计算基本相同，相对误差不超过0.2％，而变步长法缩短了计算时间，提高了计算速度。