CN101178748A - 一种有限元求解轧制过程温度场的集中热容矩阵方法 - Google Patents

Abstract

一种有限元求解轧制过程温度场的集中热容矩阵方法，通过将有限元线性方程组中热容矩阵的同行或同列元素相加代替对角线元素进行计算，在不影响计算效率的情况下有效地克服了传统有限元方法计算瞬态温度场时产生的振荡现象，保证了计算的稳定性；该方法计算的中厚板轧制过程中的钢板表面温度值与实测温度的结果对比表明，该方法保证了计算瞬态温度场时具有较高的计算精度。

Description

一种有限元求解轧制过程温度场的集中热容矩阵方法

技术领域

本发明属于轧制技术领域，特别涉及一种有限元求解轧制过程温度场的集中热容矩阵方法。

背景技术

在现代工程技术领域中，有限元法作为一种有效的数值分析方法已广泛应用于结构、力和热等的分析过程中；近年来，许多研究人员采用有限元法分析热轧和淬火等过程中瞬态温度场。在用有限元法分析瞬态温度场时，通常采用的方法是在空间上进行结构离散、在时间上采用有限差分格式进行求解。该方法从初始时刻的温度场开始，每隔一个时间步长进行一次迭代计算，进而求出各个时刻的温度场；其优点是节省计算机内存，但是往往产生时间和空间上的振荡现象，影响了计算的稳定性和精度。

对于上述振荡现象，研究人员从差分格式、网格划分方式、加权变系数、迭代变步长、在实数子空间上求解和用Norsette法代替差分格式计算等多方面进行了深入的研究，得出了一些降低或消除振荡的方法，但效果都不太理想。

发明内容

本发明的目的就是通过采用集中热容矩阵的方法克服上述采用有限元法求解瞬态温度场时常常产生的时间和空间上的振荡现象，保证计算的稳定性，提高的计算精度。

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

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

轧制参数：轧制时间，轧件宽度，轧件厚度，初始时间，初始温度，轧件周围介质温度，时间步长。

材料热物性参数：热传导系数，黑度，比热，密度。

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

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

单元和节点编号沿厚度方向和宽度方向逐渐增加，图1中，i为单元编号，j为节点编号，H为厚度，W为宽度。确定边界条件：AB和AD边界绝热，BC和CD边界换热。以A点坐标为零，计算各节点坐标，在宽度方向和厚度方向上单元均匀划分。

③根据不同轧制过程，确定边界换热系数h。

在板带轧制过程中，不同轧制阶段，具有不同换热系数计算模型。热轧板带

在空冷过程中，其表面换热方式主要为辐射和自然对流，辐射系数表述为：

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。

热轧板带在高压水除磷过程中，主要换热方式为强迫对流和侧面辐射，辐射计算方法同上，对流系数表达式为：

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

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

在轧制过程中，板带与轧辊之间接触换热是主要热损失方式。接触换热系数与氧化铁皮厚度、氧化铁皮导热率和轧制压力有关。接触换热系数表达式为：

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

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

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

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

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

$k(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{\∂y^2})-\mathrm{\ρc}\frac{\∂T}{\∂t}=0---\left(4\right)$

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

$K_{\mathrm{Tij}}={\∫\∫}_{S}k(\frac{{\∂N}_i}{\∂x}\·\frac{{\∂N}_j}{\∂x}+\frac{{\∂N}_i}{\∂y}\·\frac{{\∂N}_j}{\∂y})\mathrm{dS}+{\∫}_{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)}={\∫\∫}_{S^e}k(\frac{{\∂N}_i}{\∂x}\·\frac{{\∂N}_j}{\∂x}+\frac{{\∂N}_i}{\∂y}\·\frac{{\∂N}_j}{\∂y})\mathrm{dS}$

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

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

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

式中：

T-瞬时温度(K)

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

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

t-时间(s)

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

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

$I=\frac{1}{2}{\∫\∫}_S\left[k\right[{\left(\frac{\∂T}{\∂x}\right)}^2+{\left(\frac{\∂T}{\∂y}\right)}^2]+2\mathrm{\ρc}\frac{\∂T}{\∂t}T]\mathrm{dS}+\frac{1}{2}{\∫}_{l}h(T-T_{\∞})\mathrm{dl}---\left(5\right)$

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

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

式中：

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

K_3ij＝∫∫_SρcN_iN_jdS

p_i＝∫_LhT_∞N_idL

其中：

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

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

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

h-换热系数；

N-形函数；

i，j-节点编号；

⑥利用二点向后差分格式，将系统方程转化为瞬态温度场求解的线性方程组。将系统方程(6)中的温度对时间偏导数表示为二点向后差分格式：

$\frac{\∂T}{\∂t}=\frac{1}{\mathrm{\Δt}}(T_t-T_{t-\mathrm{\Δt}})---\left(7\right)$

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

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

⑦集中热容矩阵

式(8)中，[K₃]是由各单元质量阵和比热的乘积组装成的矩阵，称为协调质量热容矩阵(简称热容矩阵)，该矩阵为n(n为节点总数)阶对称带状矩阵。由于并非只有对角元素为零，由(1)式中各方程可以看出流入某节点的热量不仅与该节点的温度变化有关还与其周围节点的温度变化有关，这是不合理的。

对于该问题一个合理的解决办法是将单元的质量集中到各节点上去。一种可行的做法是将热容矩阵的同行或同列元素相加代替对角线元素，新的热容矩阵只有对角线元素有值，其余元素均为零，如(3)式所示。

集中热容矩阵后(1)式中各方程的物理意义是十分明确且合理的，即流入某节点的热量与该节点温度的变化相平衡。

⑧求解线性方程组，获得温度场

⑨根据迭代次数判断程序计算是否结束，如果迭代次数大于设定迭代次数则退出程序，否则继续迭代计算。

计算流程如图2所示。

本发明在传统算法的基础上通过采用集中热容矩阵的方法，在不影响计算效率的情况下有效地克服了振荡现象，能够保证了计算的稳定性，提高计算精度。

附图说明

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

图2本发明方法单元划分图，

图3本发明方法计算流程图，

图4实施本发明方法之前温度场的空间分布图，

图5实施本发明方法之前温度随时间变化图，

图6实施本发明方法之后温度场的空间分布图，

图7实施本发明方法之后温度随时间变化图。

图中：i为单元编号，j为节点编号，H为厚度，W为宽度，1为换热边界，P，Q，R为有限元分析模型上任意选择的测温点，H/2为半厚度，B/2为半宽度。

具体实施方式

以中厚板出加热炉后空冷过程中横断面二维瞬态温度场的计算为例加以说明。

(1)计算条件

钢板宽B＝2.0m，厚H＝0.22m，密度ρ＝7800.0kg/m3，比热c＝670.0J/(kg.K)、导热系数k＝30.0W/(m.K)，与空气之间的换热系数视为温度的函数，钢板的初始温度均匀分布(T0＝1150℃)，环境温度T∞＝25℃。

(2)有限元网格

假设轧件对称，取断面四分之一作为研究对象，以其对称中心为原点坐标，采用四边形等参单元，宽度方向单元数m＝30，厚度方向单元数n＝10，单元划分如图3所示。

集中热容矩阵前的计算结果

取迭代时间步长为0.5s，采用向后差分格式迭代四步后，钢板横断面温度的计算结果如图4所示，边界附近节点的温度值超出了合理的范围(T0)，产生了空间上的振荡现象。

图5为振荡区域中点P(66，99)、Q(967，99)和R(967，22)的温度值随时间的变化情况；由图4可以看出温度值在迭代开始阶段都有一个违背传热规律的上升的过程，产生了时间上的振荡现象。空间和时间上的振荡影响了计算的稳定性和精度，必须加以消除。

集中热容矩阵后的计算结果

图6和图7为对上述算例采用集中热容矩阵法后的计算结果。可见集中热容矩阵后，空间上的温度振荡得到很好的抑制；边界振荡较剧烈区中的P、Q和R点的温度值在迭代开始阶段的上升现象也基本消失。

经计算可知，采用集中热容矩阵后，用有限元法计算瞬态温度场可以获得较为稳定精确的计算结果。

Claims (1)

1.一种有限元求解轧制过程温度场的集中热容矩阵方法，其特征在于该方法包括以下步骤：

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

轧制参数：初始时间，轧制时间，轧件宽度，轧件厚度，初始温度，轧件周围介质温度，时间步长

材料热物性参数：热传导系数，黑度，比热，密度

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

(2)建立有限元分析模型，进行单元节点编号、确定换热边界和计算节点坐标

(3)根据不同轧制过程，确定边界换热系数h

热轧板带，在空冷过程中，辐射系数表述为：

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

式中：HR为辐射系数，σ＝5.67×10^-8W/(m²·K⁴)

ε＝0.125(T/1000)²-0.38(T/1000)+1.1

热轧板带在高压水除磷过程中，辐射系数HR，对流系数表达式为：

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

式中：w为水流密度，T板带表面温度

在轧制过程中，接触换热系数表达式为：

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

式中：p_m-轧制压力

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

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

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

$k(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{{\∂y}^2})-\mathrm{\ρc}\frac{\∂T}{\∂t}=0$

其中：

T-瞬时温度，ρ-材料密度，c-材料比热，t-时间

k-热传导系数

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

$I=\frac{1}{2}{\∫\∫}_S\left[k\right[{\left(\frac{\∂T}{\∂x}\right)}^2+{\left(\frac{\∂T}{\∂y}\right)}^2]+2\mathrm{\ρc}\frac{\∂T}{\∂t}T]\mathrm{dS}+\frac{1}{2}{\∫}_{l}h(T-T_{\∞})\mathrm{dl}$

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

$\left[K_T\right]\left\{T\right\}+\left[K_3\right]\left\{\frac{\∂T}{\∂t}\right\}=\left\{p\right\}$

式中：

[K_T]-温度刚度矩阵

[K₃]-变温矩阵

{p}-常数项列式

k-热传导系数

ρ-材料密度

c-材料比热

h-换热系数

N-形函数

i，j节点编号

(6)进行集中热容矩阵，求解瞬态温度场

[K₃]热容矩阵，该矩阵为n(n为节点总数)阶对称带状矩阵，将热容矩阵的同行或同列元素相加代替对角线元素，新的热容矩阵只有对角线元素有值，其余元素均为零

(7)利用二点向后差分格式，将系统方程转化为瞬态温度场求解的线性方程组二点向后差分格式 $\frac{\∂T}{\∂t}=\frac{1}{\mathrm{\Δt}}(T_t-T_{t-\mathrm{\Δt}})$

温度场求解的线性方程组

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

t-Δt温度场已知，求出t时刻的温度场，将所得温度作为初始条件，反复迭代求解下去，得出任意时刻温度场。

Images