一种快速计算轧件断面温度的有限元和有限差分耦合方法

技术领域

本发明属于轧制技术领域，特别涉及热轧过程中快速计算异形轧件断面温度的有限元与有限差分耦合方法。

背景技术

轧制过程中温度是最重要的参数之一。由于温度直接影响到轧制力和轧后轧件的弯曲，因此精确预报生产过程中各阶段的温度是控制轧制和控制冷却的关键。过去在生产中使用的有限元方法计算时间长同时只适用于简单断面的板带和棒线在线温度计算，难以适应于各种断面复杂的型钢的温度计算。采用有限元在计算复杂断面轧件的温度场时，耗时太多、轧件咬入难，这导致采用有限元法来快速计算复杂断面轧件轧制过程中的温度受到限制。而有限差分法过去主要应用在网格不变的工模具或变形规律简单的板带和棒线材的温度计算中。复杂断面的轧件无法预测轧制过程前后的金属流动规律，而无法应用。有限差分法能快速的计算温度，而有限元法能准确的反映轧制过程金属流动规律，将两种方法耦合在一起在线快速计算复杂断面轧件的温度成为可能，但是目前还没有将有限元和有限差分两种方法耦合在一起应用到复杂断面轧制过程中。

公开号为CN 100495411C的中国专利公开了一种预测热轧过程板带温度场的有限元方法，它包括采集轧制过程数据；对板带断面进行四边形有限单元划分；确定边界换热系数和内热源强度；确定四边形单元的形函数、B矩阵和雅克比矩阵J；确定有限单元的温度刚度矩阵和变温矩阵的组装；一维变带宽存储法求解方程组。该方法可以较快的获得热轧板带轧制过程轧件断面的瞬态温度场。但是该方法只应用在轧件断面形状简单的板带，复杂断面轧件的热轧过程的温度计算无法进行有效快速计算。

公开号为CN101046682 A的中国专利公开了一种预测热轧含Nb带钢组织及力学性能的方法，它用有限差分法建立温度模型，包括辊道上轧件温度模型、粗轧段轧件温度模型、精轧段轧件温度模型、层流冷却段轧件温度模型。该轧件断面温度计算模型也只适合于轧件断面形状简单的板带钢，复杂断面轧件的热轧过程温度计算也无法进行有效快速计算。

发明内容

本发明的目的就是克服计算复杂断面轧件温度计算的有限元方法时间过长、难自动化的缺点，克服有限差分法无法预测轧件轧制前后金属流动规律的缺点，提高轧件断面温度计算的速度和自动化程度，提出一种快速计算轧件断面的有限元和有限差分耦合方法。

本发明所采用的技术方案是：一种快速计算轧件断面温度的有限元和有限差分耦合方法，其特征在于，包括以下步骤：

步骤1：对坯料横断面进行有限差分网格划分，建立有限差分分析模型；

步骤2：采用有限差分模型计算轧件从出炉、除鳞、及第一道次轧制前的断面温度；

步骤3：对第一道次轧制前的轧件进行有限元网格划分，建立三维有限元模型并进行计算，获得第一道次轧制前后轧件的有限元单元及节点数据；

步骤4：建立有限元和有限差分耦合模型，其具体实现包括以下子步骤：

步骤4.1：将第一道次轧制前轧件有限差分网格温度映射到第一道次轧制前轧件有限元网格上；

步骤4.2：将第一道次轧制前轧件有限元节点温度传递到第一道次轧制后轧件的有限元节点上，并将塑性变形热导致的轧件温度升高△T加到第一道次轧制后轧件的所有有限元节点上；

步骤4.3：对预先采用有限元方法计算好的第一道次轧制后的轧件断面进行有限差分网格划分，将第一道次轧制后轧件的有限元节点温度映射到第一道次轧制后轧件的有限差分网格上；

步骤5：采用第一道次轧制后轧件的有限差分网格进行轧制间隙的温度计算；

步骤6：以后各个道次轧制过程都采用步骤1对轧件横断面进行有限差分网格划分和步骤4的耦合模型进行计算，轧制间隙及水雾冷却的温度采用步骤5对轧件横断面进行有限差分温度计算，最终获得轧件各个时刻的断面温度。

作为优选，步骤1中所述的对坯料横断面进行有限差分网格划分，建立有限差分分析模型；其具体实现包括以下子步骤：

步骤1.1：在坯料横断面高度和宽度方向上采用均匀网格，网格大小根据轧件断面的复杂程度来调节；有限差分网格范围为矩形区域，整个有限差分网格覆盖整个坯料横断面，为了便于判断边界点，有限差分网格比坯料横断面大两个节点宽度和高度；

步骤1.2：构建内部节点的二维有限差分方程；

对于二维常物性控制方程(无热源)(热传导平衡方程)

$>\mathrm{\rho c}\frac{\partial T}{\partial \tau }=\lambda (\frac{{\partial }^{2}T}{{\partial x}^2}+\frac{{\partial }^{2}T}{{\partial y}^2})---\left(1\right);$>

将式(1)中对时间的一阶微商用一阶向前差商代替，对坐标的二阶微商用二阶中心差商代替获得内部节点的显式差分方程为：

$>T_{m,n}^{k+1}=\frac{\mathrm{a\Delta \tau }}{{\left(\mathrm{\Delta x}\right)}^2}(T_{m+1,n}^k-2T_{m,n}^k+T_{m-1,n}^k)+\frac{\mathrm{a\Delta \tau }}{{\left(\mathrm{\Delta y}\right)}^2}(T_{m,n+1}^k-2T_{m,n}^k+T_{m,n-1}^k)+T_{m,n}^k---\left(2\right);$>

这里，Δτ为时间步长，如果是均匀离散Δx＝Δy，则有：

$>T_{m,n}^{k+1}=\mathrm{Fo}(T_{m+1,n}^k+T_{m-1,n}^k+T_{m,n+1}^k+T_{m,n-1}^k)+(1-4\mathrm{Fo})T_{m,n}^k---\left(3\right);$>

$>\mathrm{Fo}=\frac{\mathrm{a\Delta \tau }}{{\left(\mathrm{\Delta x}\right)}^2}=\frac{\lambda }{\mathrm{\rho c}}\frac{\mathrm{\Delta \tau }}{{\left(\mathrm{\Delta x}\right)}^2}---\left(4\right);$>

每个时间层k+1的节点(i，j)温度都表达成了上一时间层k的以该节点(i，j)为中心的5个节点温度的显函数，该格式的稳定性条件是：Fo≤1/4；

步骤1.3：构建边界节点的二维有限差分方程；

根据边界节点所处的位置不同，边界节点的二维有限差分方程分为12种；具体如下：

(1)平直边界的右边界：

当边界节点位于右平直边界时，对边界采用热平衡法(即传入节点的热流等于该节点所代表控制的内能变化率)导出二维有限差分边界节点T_m,n热平衡方程为：

$>\mathrm{\lambda \Delta y}\frac{T_{m,n-1}^{\left(k\right)}-T_{m,n}^{\left(k\right)}}{\mathrm{\Delta x}}+\lambda \frac{\mathrm{\Delta x}}{2}\frac{T_{m+1,n}^{\left(k\right)}-T_{m,n}^{\left(k\right)}}{\mathrm{\Delta y}}+\lambda \frac{\mathrm{\Delta x}}{2}\frac{T_{m-1,n}^{\left(k\right)}-T_{m,n}^{\left(k\right)}}{\mathrm{\Delta y}}-h(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})\mathrm{\Delta y}-ϵc_0[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]\mathrm{\Delta y}=\mathrm{\rho c}\frac{\mathrm{\Delta x}}{2}\mathrm{\Delta y}\frac{T_{m,n}^{(k+1)}-T_{m,n}^{\left(k\right)}}{\mathrm{\Delta \tau }}---\left(5\right);$>

取Δx＝Δy得

$>T_{m,n}^{k+1}=F_0(2T_{m,.n-1}^k+T_{m+1,n}^k+T_{m-1,n}^k)+T_{m,n}^k(1-4F_0)-\frac{2\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^k-T_{\mathrm{fluid}})-\frac{2ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^4-T_e^4)---\left(6\right);$>

(2)平直边界的左边界：

$>T_{m,n}^{k+1}=F_0(2T_{m,.n+1}^k+T_{m+1,n}^k+T_{m-1,n}^k)+T_{m,n}^k(1-4F_0)-\frac{2\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^k-T_{\mathrm{fluid}})-\frac{2ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^4-T_e^4)---\left(7\right);$>

(3)平直边界的上边界：

$>T_{m,n}^{k+1}=F_0(2T_{m-1,n}^k+T_{m,n-1}^k+T_{m,n+1}^k)+T_{m,n}^k(1-4F_0)-\frac{2\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^k-T_{\mathrm{fluid}})-\frac{2ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^4-T_e^4)---\left(8\right);$>

(4)平直边界的下边界：

$>T_{m,n}^{k+1}=F_0(2T_{m+1,n}^k+T_{m,n-1}^k+T_{m,n+1}^k)+T_{m,n}^k(1-4F_0)-\frac{2\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^k-T_{\mathrm{fluid}})-\frac{2ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^4-T_e^4)---\left(9\right);$>

(5)凸角边界的右上凸角：

$>T_{m,n}^{(k+1)}=\frac{F_0}{2}(T_{m,n-1}^{\left(k\right)}+T_{m-1,n}^{\left(k\right)})+T_{m,n}^{\left(k\right)}(1-F_0)-\frac{\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]---\left(10\right);$>

(6)凸角边界的右下凸角：

$>T_{m,n}^{(k+1)}=\frac{F_0}{2}(T_{m,n-1}^{\left(k\right)}+T_{m+1,n}^{\left(k\right)})+T_{m,n}^{\left(k\right)}(1-F_0)-\frac{\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]---\left(11\right);$>

(7)凸角边界的左上凸角：

$>T_{m,n}^{(k+1)}=\frac{F_0}{2}(T_{m,n+1}^{\left(k\right)}+T_{m-1,n}^{\left(k\right)})+T_{m,n}^{\left(k\right)}(1-F_0)-\frac{\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]---\left(12\right);$>

(8)凸角边界的左下凸角：

$>T_{m,n}^{(k+1)}=\frac{F_0}{2}(T_{m,n+1}^{\left(k\right)}+T_{m+1,n}^{\left(k\right)})+T_{m,n}^{\left(k\right)}(1-F_0)-\frac{\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]---\left(13\right);$>

(9)凹角边界的右下角：

$>T_{m,n}^{(k+1)}=\frac{4}{3}{F}_0(T_{m,n-1}^{\left(k\right)}+T_{m+1,n}^{\left(k\right)}+\frac{T_{m-1,n}^{\left(k\right)}}{2}+\frac{T_{m,n+1}^{\left(k\right)}}{2})+T_{m,n}^{\left(k\right)}(1-4F_0)-\frac{4\mathrm{h\Delta \tau }}{3\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{4ϵc_0\mathrm{\Delta \tau }}{3\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-T_e^4]---\left(14\right);$>

(10)凹角边界的右上角：

$>T_{m,n}^{(k+1)}=\frac{4}{3}{F}_0(T_{m,n-1}^{\left(k\right)}+T_{m-1,n}^{\left(k\right)}+\frac{T_{m+1,n}^{\left(k\right)}}{2}+\frac{T_{m,n+1}^{\left(k\right)}}{2})+T_{m,n}^{\left(k\right)}(1-4F_0)-\frac{4\mathrm{h\Delta \tau }}{3\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{4ϵc_0\mathrm{\Delta \tau }}{3\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-T_e^4]---\left(15\right);$>

(11)凹角边界的左上角：

$>T_{m,n}^{(k+1)}=\frac{4}{3}{F}_0(T_{m,n+1}^{\left(k\right)}+T_{m-1,n}^{\left(k\right)}+\frac{T_{m+1,n}^{\left(k\right)}}{2}+\frac{T_{m,n-1}^{\left(k\right)}}{2})+T_{m,n}^{\left(k\right)}(1-4F_0)-\frac{4\mathrm{h\Delta \tau }}{3\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{4ϵc_0\mathrm{\Delta \tau }}{3\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-T_e^4]---\left(16\right);$>

(12)凹角边界的左下角：

$>T_{m,n}^{(k+1)}=\frac{4}{3}{F}_0(T_{m,n+1}^{\left(k\right)}+T_{m+1,n}^{\left(k\right)}+\frac{T_{m-1,n}^{\left(k\right)}}{2}+\frac{T_{m,n-1}^{\left(k\right)}}{2})+T_{m,n}^{\left(k\right)}(1-4F_0)-\frac{4\mathrm{h\Delta \tau }}{3\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{4ϵc_0\mathrm{\Delta \tau }}{3\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-T_e^4]---\left(17\right);$>

其中，ρ为密度,单位为kg/m³；c为比热容，单位为J/(kg·℃)；λ为导热率，单位为W/(m·K)；T为温度，单位为K；x、y为横纵坐标，单位为m；Δx、Δy为有限差分网格x、y方向的步距，单位为m；τ为时间，单位为s；Δτ为时间步长，单位为s；a为导温系数，单位为m²/s；为有限差分网格节点(m,n)在k+1时刻的温度，单位为K；Fo为傅里叶数；h为轧件与环境的换热系数，单位为W/(m²·K)；ε为黑度；c₀为黑体的辐射系数；T_e为辐射环境温度，单位为K；T_fluid对流环境温度或水温，单位为K。

作为优选，步骤3中所述的对第一道次轧制前的轧件进行有限元网格划分，建立三维有限元模型并进行有限元计算，采用的是弹塑性有限元理论。

作为优选，步骤2中所述的轧件从出炉、除鳞、及第一道次轧制前的断面温度计算，步骤3中第一道次有限元网格划分及计算，步骤4中所述的建立有限元和有限差分耦合模型进行轧制过程计算，步骤5中所述轧制间隙有限差分温度计算以及步骤6所述后续轧制道次、间隙和水雾冷却温度计算，计算所需要数据包括：坯料参数、物性参数、设备参数、轧制工艺参数、综合影响参数。

作为优选，所述的坯料参数包括坯料初始长度、坯料初始宽度、坯料初始厚度、坯料出炉温度、坯料圆角半径；所述的物性参数包括热传导系数、比热、密度、黑度；所述的设备参数包括各个设备的间距、除鳞及水雾冷却设备参数；所述的轧制工艺参数包括轧制各阶段的轧制时间和间隙时间；所述的综合影响参数包括单元步长、时间步长、环境温度、各阶段的换热系数。

作为优选，所述的坯料参数的确定由操作者给出，所述的物性参数的确定由材料的成分和温度确定，所述的设备参数由具体车间布置确定，所述的轧制工艺参数由具体规格的轧制产品的轧制程序表确定。

作为优选，所述的各阶段的换热系数包括空冷换热模型、除鳞换热系数、水雾冷却换热系数、轧制换热模型。

作为优选，所述的空冷换热模型：

轧件在放置或在输送辊道上运送时发生的热交换主要以辐射和对流的形式而产生的，对流系数为：

h＝η(T₀-T_fluid)^0.25  (18)；

轧件辐射换热热量为：

Φ_r＝εc₀[(T₀)⁴-(T_e)⁴]A   (19)；

轧件的辐射黑度ε为小于1的数，

$>ϵ=\lambda [1.1+\frac{T_0}{1000}(0.125\frac{T_0}{1000}-0.38)]---\left(20\right);$>

所述的除鳞换热系数的确定：

在高压水除鳞过程中，水流密度、水压和轧件表面温度对换热系数的影响比较大，其中表面换热方式主要为强迫对流，对流系数表达式为：

$>h_c=\gamma 107.2w^{0.663}\times {10}^{-0.00147T_0}\times 1.163---\left(21\right);$>所述的水雾冷却换热系数的确定：

轧件局部冷却时，为了防止水的溅射对不需要冷却部位的影响，采用冷却能力较小同时没有溅射的水雾冷却方式；水雾冷却的换热系数表达式为：

lgh_sw＝6.615-0.563lg(H/D)+0.200lgP-0.865lgT₀  (22)；所述的轧制换热模型：

轧制过程作为整体处理，轧制时轧件塑性变形所产生的热量造成温度升高ΔT：

ΔT＝0.183σlnμ  (23)；

根据式(23)可以算出整个轧制过程每个轧制道次的温升，将该温度均匀加到有限差分程序中。

其中，T₀为轧件表面温度，单位为K；η为修正系数；Φ_r为辐射交换热量，单位为W；A为辐射面积，单位为m²；λ为修正系数；h_c为除鳞换热系数，单位为W/(m²·K)；w为水流密度，单位为L/(min·m²)；γ为水压影响系数；h_sw为水雾冷却换热系数，单位为W/(m²·K)；H为喷嘴到钢板距离，单位为mm；D为喷嘴直径，单位为mm，P为水压，单位为MPa；σ为金属变形抗力，单位为MPa，ΔT为轧件在孔型中塑性变形引起的温升，单位为K；μ为延伸系数。

相对于现有技术，本发明具有如下积极效果：

(1)可以快速计算各种断面的轧件轧制过程任何时刻的瞬态温度场；

(2)传统有限元方法计算轧件轧制过程的温度场的时间以天计算，采用本发明计算时间以分钟计算，计算速度提高非常大；

(3)有限差分方法无法计算复杂断面轧件温度；而有限元方法很难连续快速计算复杂断面的轧件断面温度，计算中间需要人工干预，自动化程度不高；本发明计算过程不需要人工干预，自动化程度高。

附图说明

图1：本发明实施例的坯料几何坐标和有限差分坐标(有限差分网格范围比坯料横断面几何大)示意图；

图2：本发明实施例的内部节点有限差分网格示意图；

图3：本发明实施例的右平直边界节点示意图；

图4：本发明实施例的右上凸角边界节点示意图；

图5：本发明实施例的右下凹角边界节点示意图；

图6：本发明实施例的某一轧制道次FEM和FDM耦合温度计算流程图；

图7：本发明实施例的重轨生产的工艺布置图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图7，为实施例的重轨生产的工艺布置图，其轧制工艺过程如下：坯料在加热炉加热到出炉温度--高压水除磷——BD1开坯轧制——BD2异形轧制——水雾冷却——万能轧制粗轧——水雾冷却——万能轧制精轧——冷床冷却——矫直。

该轧制过程可以快速计算某工艺条件下轧件断面温度，通过调整水雾冷却参数可以控制精轧后轧件断面温度的均匀程度，防止轧件轧后断面温度不均而在冷床上产生弯曲，从而无法进行矫直。

为了快速计算轧件断面温度，本实施例采用的技术方案包括以下步骤：

步骤1：对坯料横断面进行有限差分网格划分，建立有限差分分析模型；建立有限差分分析模型具体实现包括以下子步骤：

步骤1.1：建立有限差分分析模型如图1所示，在坯料横断面高度和宽度方向上采用均匀网格，网格大小根据轧件断面的复杂程度来调节；有限差分网格范围为矩形区域，整个有限差分网格覆盖整个坯料横断面，为了便于判断边界点，有限差分网格比坯料横断面大两个节点宽度和高度；

步骤1.2：构建内部节点的二维有限差分方程；

对于二维常物性控制方程(无热源)(热传导平衡方程)

$>\mathrm{\rho c}\frac{\partial T}{\partial \tau }=\lambda (\frac{{\partial }^{2}T}{{\partial x}^2}+\frac{{\partial }^{2}T}{{\partial y}^2})---\left(1\right);$>

图2为内部节点有限差分网格，计算内部节点(m,n)的温度涉及周围4个有限差分节点(m,n+1)、(m,n-1)、(m-1,n)、(m+1,n)。将式(1)中对时间的一阶微商用一阶向前差商代替，对坐标的二阶微商用二阶中心差商代替获得内部节点的显式差分方程为：

$>T_{m,n}^{k+1}=\frac{\mathrm{a\Delta \tau }}{{\left(\mathrm{\Delta x}\right)}^2}(T_{m+1,n}^k-2T_{m,n}^k+T_{m-1,n}^k)+\frac{\mathrm{a\Delta \tau }}{{\left(\mathrm{\Delta y}\right)}^2}(T_{m,n+1}^k-2T_{m,n}^k+T_{m,n-1}^k)+T_{m,n}^k---\left(2\right);$>

这里，Δτ为时间步长，如果是均匀离散Δx＝Δy，则有：

$>T_{m,n}^{k+1}=\mathrm{Fo}(T_{m+1,n}^k+T_{m-1,n}^k+T_{m,n+1}^k+T_{m,n-1}^k)+(1-4\mathrm{Fo})T_{m,n}^k---\left(3\right);$>

$>\mathrm{Fo}=\frac{\mathrm{a\Delta \tau }}{{\left(\mathrm{\Delta x}\right)}^2}=\frac{\lambda }{\mathrm{\rho c}}\frac{\mathrm{\Delta \tau }}{{\left(\mathrm{\Delta x}\right)}^2}---\left(4\right);$>

每个时间层k+1的节点(i，j)温度都表达成了上一时间层k的以该节点(i，j)为中心的5个节点温度的显函数，该格式的稳定性条件是：Fo≤1/4；

步骤1.3：构建边界节点的二维有限差分方程；

根据边界节点所处的位置不同，边界节点的二维有限差分方程分为12种；具体如下：

(1)平直边界的右边界：

当边界节点位于右平直边界时，参与温度计算的节点如图3所示，计算边界节点(m,n)的温度涉及周围3个有限差分节点(m,n-1)、(m-1,n)、(m+1,n)。对边界采用热平衡法(即传入节点的热流等于该节点所代表控制的内能变化率)导出二维有限差分边界节点T_m,n热平衡方程为：

$>\mathrm{\lambda \Delta y}\frac{T_{m,n-1}^{\left(k\right)}-T_{m,n}^{\left(k\right)}}{\mathrm{\Delta x}}+\lambda \frac{\mathrm{\Delta x}}{2}\frac{T_{m+1,n}^{\left(k\right)}-T_{m,n}^{\left(k\right)}}{\mathrm{\Delta y}}+\lambda \frac{\mathrm{\Delta x}}{2}\frac{T_{m-1,n}^{\left(k\right)}-T_{m,n}^{\left(k\right)}}{\mathrm{\Delta y}}-h(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})\mathrm{\Delta y}-ϵc_0[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]\mathrm{\Delta y}=\mathrm{\rho c}\frac{\mathrm{\Delta x}}{2}\mathrm{\Delta y}\frac{T_{m,n}^{(k+1)}-T_{m,n}^{\left(k\right)}}{\mathrm{\Delta \tau }}---\left(5\right);$>

取Δx＝Δy得

$>T_{m,n}^{k+1}=F_0(2T_{m,.n-1}^k+T_{m+1,n}^k+T_{m-1,n}^k)+T_{m,n}^k(1-4F_0)-\frac{2\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^k-T_{\mathrm{fluid}})-\frac{2ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^4-T_e^4)---\left(6\right);$>

(2)平直边界的左边界：

$>T_{m,n}^{k+1}=F_0(2T_{m,.n+1}^k+T_{m+1,n}^k+T_{m-1,n}^k)+T_{m,n}^k(1-4F_0)-\frac{2\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^k-T_{\mathrm{fluid}})-\frac{2ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^4-T_e^4)---\left(7\right);$>

(3)平直边界的上边界：

$>T_{m,n}^{k+1}=F_0(2T_{m-1,n}^k+T_{m,n-1}^k+T_{m,n+1}^k)+T_{m,n}^k(1-4F_0)-\frac{2\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^k-T_{\mathrm{fluid}})-\frac{2ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^4-T_e^4)---\left(8\right);$>

(4)平直边界的下边界：

$>T_{m,n}^{k+1}=F_0(2T_{m+1,n}^k+T_{m,n-1}^k+T_{m,n+1}^k)+T_{m,n}^k(1-4F_0)-\frac{2\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^k-T_{\mathrm{fluid}})-\frac{2ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^4-T_e^4)---\left(9\right);$>

(5)凸角边界的右上凸角：

当边界节点位于右上凸角边界时，参与温度计算的节点如图4所示，计算边界节点(m,n)的温度涉及周围2个有限差分节点(m,n-1)、(m-1,n)。

$>T_{m,n}^{(k+1)}=\frac{F_0}{2}(T_{m,n-1}^{\left(k\right)}+T_{m-1,n}^{\left(k\right)})+T_{m,n}^{\left(k\right)}(1-F_0)-\frac{\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]---\left(10\right);$>

(6)凸角边界的右下凸角：

$>T_{m,n}^{(k+1)}=\frac{F_0}{2}(T_{m,n-1}^{\left(k\right)}+T_{m+1,n}^{\left(k\right)})+T_{m,n}^{\left(k\right)}(1-F_0)-\frac{\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]---\left(11\right);$>

(7)凸角边界的左上凸角：

$>T_{m,n}^{(k+1)}=\frac{F_0}{2}(T_{m,n+1}^{\left(k\right)}+T_{m-1,n}^{\left(k\right)})+T_{m,n}^{\left(k\right)}(1-F_0)-\frac{\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]---\left(12\right);$>

(8)凸角边界的左下凸角：

$>T_{m,n}^{(k+1)}=\frac{F_0}{2}(T_{m,n+1}^{\left(k\right)}+T_{m+1,n}^{\left(k\right)})+T_{m,n}^{\left(k\right)}(1-F_0)-\frac{\mathrm{h\Delta \tau }}{\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{ϵc_0\mathrm{\Delta \tau }}{\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-{\left(T_e\right)}^4]---\left(13\right);$>

(9)凹角边界的右下角：

当边界节点位于右下凹角边界时，参与温度计算的节点如图5所示，计算边界节点(m,n)的温度涉及周围4个有限差分节点(m,n+1)、(m,n-1)、(m-1,n)、(m+1,n)。

$>T_{m,n}^{(k+1)}=\frac{4}{3}{F}_0(T_{m,n-1}^{\left(k\right)}+T_{m+1,n}^{\left(k\right)}+\frac{T_{m-1,n}^{\left(k\right)}}{2}+\frac{T_{m,n+1}^{\left(k\right)}}{2})+T_{m,n}^{\left(k\right)}(1-4F_0)-\frac{4\mathrm{h\Delta \tau }}{3\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{4ϵc_0\mathrm{\Delta \tau }}{3\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-T_e^4]---\left(14\right);$>

(10)凹角边界的右上角：

$>T_{m,n}^{(k+1)}=\frac{4}{3}{F}_0(T_{m,n-1}^{\left(k\right)}+T_{m-1,n}^{\left(k\right)}+\frac{T_{m+1,n}^{\left(k\right)}}{2}+\frac{T_{m,n+1}^{\left(k\right)}}{2})+T_{m,n}^{\left(k\right)}(1-4F_0)-\frac{4\mathrm{h\Delta \tau }}{3\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{4ϵc_0\mathrm{\Delta \tau }}{3\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-T_e^4]---\left(15\right);$>

(11)凹角边界的左上角：

$>T_{m,n}^{(k+1)}=\frac{4}{3}{F}_0(T_{m,n+1}^{\left(k\right)}+T_{m-1,n}^{\left(k\right)}+\frac{T_{m+1,n}^{\left(k\right)}}{2}+\frac{T_{m,n-1}^{\left(k\right)}}{2})+T_{m,n}^{\left(k\right)}(1-4F_0)-\frac{4\mathrm{h\Delta \tau }}{3\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{4ϵc_0\mathrm{\Delta \tau }}{3\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-T_e^4]---\left(16\right);$>

(12)凹角边界的左下角：

$>T_{m,n}^{(k+1)}=\frac{4}{3}{F}_0(T_{m,n+1}^{\left(k\right)}+T_{m+1,n}^{\left(k\right)}+\frac{T_{m-1,n}^{\left(k\right)}}{2}+\frac{T_{m,n-1}^{\left(k\right)}}{2})+T_{m,n}^{\left(k\right)}(1-4F_0)-\frac{4\mathrm{h\Delta \tau }}{3\mathrm{\rho c\Delta x}}(T_{m,n}^{\left(k\right)}-T_{\mathrm{fluid}})-\frac{4ϵc_0\mathrm{\Delta \tau }}{3\mathrm{\rho c\Delta x}}[{\left(T_{m,n}^{\left(k\right)}\right)}^4-T_e^4]---\left(17\right);$>

其中，ρ为密度,单位为kg/m³；c为比热容，单位为J/(kg·℃)；λ为导热率，单位为W/(m·K)；T为温度，单位为K；x、y为横纵坐标，单位为m；Δx、Δy为有限差分网格x、y方向的步距，单位为m；τ为时间，单位为s；Δτ为时间步长，单位为s；a为导温系数，单位为m²/s；为有限差分网格节点(m,n)在k+1时刻的温度，单位为K；Fo为傅里叶数；h为轧件与环境的换热系数，单位为W/(m²·K)；ε为黑度；c₀为黑体的辐射系数；T_e为辐射环境温度，单位为K；T_fluid对流环境温度或水温，单位为K。

步骤2：采用有限差分方法计算轧件从出炉、除鳞、及第一道次轧制前的断面温度；

步骤3：对第一道次轧制前的轧件进行有限元网格划分，采用弹塑性有限元理论建立三维有限元模型并进行计算，获得第一道次轧制前后轧件的有限元单元及节点数据；

利用弹塑性有限元理论建立各个轧制道次的三维有限元模型，异形道次轧件很难咬入下一个道次，基本上每个道次都需要单独模拟。在网格畸变严重的道次，对轧件进行网格重划。采用ANSYS APDL语言、几何模型更新法、LS-PREPOST和VB软件导出轧制前后某一截面的单元数据及节点坐标，并保存为单独的文件。

步骤4：建立有限元和有限差分耦合模型，请见图6，其具体实现包括以下子步骤：

步骤4.1：将第一道次轧制前轧件有限差分网格温度映射到第一道次轧制前轧件有限元网格上；

步骤4.2：将第一道次轧制前轧件有限元节点温度传递到第一道次轧制后轧件的有限元节点上，并将塑性变形热导致的轧件温度升高△T加到第一道次轧制后轧件的所有有限元节点上；

步骤4.3：对预先采用有限元方法计算好的第一道次轧制后的轧件断面进行有限差分网格划分，将第一道次轧制后轧件的有限元节点温度映射到第一道次轧制后轧件的有限差分网格上；

步骤5：采用第一道次轧制后轧件的有限差分网格进行轧制间隙的温度计算；

步骤6：以后各个道次轧制过程都采用步骤1对轧件横断面进行有限差分网格划分和步骤4的耦合模型进行计算，轧制间隙及水雾冷却的温度采用步骤5对轧件横断面进行有限差分温度计算，最终获得轧件各个时刻的断面温度。

其中步骤2中轧件从出炉、除鳞、及第一道次轧制前的断面温度计算，步骤3中第一道次有限元网格划分及计算，步骤4中所述的建立有限元和有限差分耦合模型进行轧制过程计算，步骤5中所述轧制间隙有限差分温度计算以及步骤6所述后续轧制道次、间隙和水雾冷却温度计算，计算所需要数据包括：坯料参数、物性参数、设备参数、轧制工艺参数、综合影响参数。

坯料参数包括坯料初始长度、坯料初始宽度、坯料初始厚度、坯料出炉温度、坯料圆角半径；所述的物性参数包括热传导系数、比热、密度、黑度；所述的设备参数包括各个设备的间距、除鳞及水雾冷却设备参数；所述的轧制工艺参数包括轧制各阶段的轧制时间和间隙时间；所述的综合影响参数包括单元步长、时间步长、环境温度、各阶段的换热系数。

坯料参数的确定由操作者给出，物性参数的确定由材料的成分和温度确定，设备参数由具体车间布置确定，轧制工艺参数由具体规格的轧制产品的轧制程序表确定。

各阶段的换热系数包括空冷换热模型、除鳞换热系数、水雾冷却换热系数、轧制换热模型。

空冷换热模型：轧件在放置或在输送辊道上运送时发生的热交换主要以辐射和对流的形式而产生的，对流系数为：

h＝η(T₀-T_fluid)^0.25  (18)；

轧件辐射换热热量为：

Φ_r＝εc₀[(T₀)⁴-(T_e)⁴]A  (19)；

轧件的辐射黑度ε为小于1的数，

$>ϵ=\lambda [1.1+\frac{T_0}{1000}(0.125\frac{T_0}{1000}-0.38)]---\left(20\right);$>

对于热轧件来说，要视其表面上的氧化铁皮的程度不同而取值也不同，当氧化铁皮较多时一般取0.8，而刚轧出的平滑表面一般取0.55～0.65。

除鳞换热系数的确定：

在高压水除鳞过程中，水流密度、水压和轧件表面温度对换热系数的影响比较大，其中表面换热方式主要为强迫对流，对流系数表达式为：

$>h_c=\gamma 107.2w^{0.663}\times {10}^{-0.00147T_0}\times 1.163---\left(21\right);$>

水雾冷却换热系数的确定：

轧件局部冷却时，为了防止水的溅射对不需要冷却部位的影响，采用冷却能力较小同时没有溅射的水雾冷却方式；水雾冷却的换热系数表达式为：

lgh_sw＝6.615-0.563lg(H/D)+0.200lgP-0.865lgT₀  (22)；

轧制换热模型：

轧制过程中，轧辊和轧件的接触面会产生摩擦热和热传递，摩擦热导致轧件温度升高，热传递导致轧件温度降低，由于轧件某一断面和轧辊接触时间很短，同时接触面很难确定，我们认为这两部分热量相互抵消，这里不予考虑。这里只考虑变形热。

轧制过程作为整体处理，轧制时轧件塑性变形所产生的热量造成温度升高ΔT：

ΔT＝0.183σlnμ  (23)；

根据式(23)可以算出整个轧制过程每个轧制道次的温升，将该温度均匀加到有限差分程序中。

其中，T₀为轧件表面温度，单位为K；η为修正系数；Φ_r为辐射交换热量，单位为W；A为辐射面积，单位为m²；λ为修正系数；h_c为除鳞换热系数，单位为W/(m²·K)；w为水流密度，单位为L/(min·m²)；γ为水压影响系数；h_sw为水雾冷却换热系数，单位为W/(m²·K)；H为喷嘴到钢板距离，单位为mm；D为喷嘴直径，单位为mm，P为水压，单位为MPa；σ为金属变形抗力，单位为MPa，ΔT为轧件在孔型中塑性变形引起的温升，单位为K；μ为延伸系数。

该轧制过程可以快速计算某工艺条件下轧件断面温度，通过调整水雾冷却参数可以控制精轧后轧件断面温度的均匀程度，防止轧件轧后断面温度不均而在冷床上产生弯曲，从而无法进行矫直。

本实施例借助编程软件VB将计算所参数作为可输入变量，将有限差分模型及方程，各个道次轧制有限元断面单元数据及节点坐标，有限差分节点和有限元节点耦合算法编制在程序中，实现快速自动的计算异形断面轧件断面温度场。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。