Ξ求解温度场的非线性有限元方法刘福来1,　杜瑞燕2(1.东北大学信息科学与工程学院,辽宁沈阳　110004;2.河北青年干部管理学院教务处,河北石家庄　050031)摘要:从G alerkin 有限元方法出发,对自由表面上的辐射换热的数学表达式不作线性化处理,而是把温度场的求解问题转化为非线性代数方程组的求解问题,并且用Newton 迭代法计算了温度场.关键词:温度场;有限元方法;Newton 迭代法中图分类号:O 242.21 文献标识码:A 文章编号:100025854(2005)0120021204由文献[1]知,求解二维待轧过程的温度场,就是要求下面微分方程和初值问题的解: 52T 5x 2+52T 5y2=1α5T5t ;(1) -k 5T5n=0,(x ,y )∈S 2;(2) -k 5T 5n=σεA (T 4-T 4∞),(x ,y )∈S 3;(3) T (x ,y ,0)=T 0(x ,y ).(4)其中:α=λρc称为导温系数,λ,ρ和c 分别为热导系数、密度和比热;S 2为给出热流强度Q 的边界面;T ∞为环境温度;S 3为给出热损失的边界面.对轧制问题的温度场,常常考虑的几种边界面[1]是:对称面、自由表面和轧件与轧辊的接触面.在辐射面上,边界条件的数学表达式为σεA (T 4-T 4∞)(其中:σ为Stefan 2Boltzmann 常数,ε为物体表面黑度,A 为辐射面积,T ∞为环境温度)是温度T 的4次幂,具有强烈的非线性.以往在实际计算中有2种处理方法[2],一种是简化问题的物理模型,有时将表达式看成常数,有时将边界条件转化成h r A (T -T ∞)(其中h r =σε(T 2+T 2∞)(T +T ∞)),在轧制问题中求解温度场时文献[1,3]都采用了这一方法;另一种是处理问题的数学方法,即用近似方法求解非线性的偏微分方程问题.例如,用数值分析的方法,文献[4]中利用了差分方法.本文中,笔者从G alerkin 有限元法出发,对自由表面上辐射换热的数学表达式不作线性处理,而是直接对非线性代数方程组用Newton 迭代法计算温度场,以二维待轧过程温度场的有限元解析进行讨论.1　G alerkin 有限元方法简介将待求解区域Ω剖分为E 个单元,每个单元4个节点.设N i 是形函数(i =1,2,3,4),用4节点线性等参单元,则单元内的温度为 T e =N 1T 1+N 2T 2+N 3T 3+N 4T 4={N }T {T}e ,(5)其中:{N }=(N 1,N 2,N 3,N 4)T ;{T}e =(T 1,T 2,T 3,T 4)T .设ω1,ω2,…,ωn 是一组基函数,用G alerkin 方法求方程(1)～(4)的解,实际上是求c 1,c 2,…,c n ,使T n =c 1ω1+c 2ω2+…+c nωn 满足 κΩρc 5T n 5t -k 52T n 5x 2+52T n5y 2ωi d x d y =0,i =1,2,…,n.(6)对式(6)应用Green 公式,有Ξ收稿日期:20040105;修回日期:20040420作者简介:刘福来(1975),男,河北省唐山市人,东北大学博士研究生.第29卷第1期2005年　1月河北师范大学学报(自然科学版)Journal of Hebei Normal University (Natural Science Edition )Vol.29No.1Jan.2005κΩρc5T n 5t ωi +k 5ωi 5x 5T n 5x +5ωi 5y 5T n5yd x d y -∫S 3k5T nn ωid s =0,i =1,2,…,n.(7)取基函数ω1,ω2,…,ωn 为形函数,即求{T}e ,使得T e ={N }T {N }e 满足 ∑Ee =1κΩeρc 5T etN i +k 5N i 5x 5T e 5x +5N i 5y 5T e 5yd x d y -∫S 3k5T nn ωid s =0,i =1,2,…,n.(8)将边界条件代入式(8)且对5Te 5t作差分后有下面的等式: ∑E e =1κΩeρc T e -T et -Δt Δt N i +k 5N i 5x 5T e 5x +5N i 5y 5T e5y d x d y + ∫S 3eσεA [({N }T {T}e )4-T 4∞]N i d s =0,i =1,2,…,n.(9)2　非线性有限元求解公式这里从式(9)出发,构造出Newton 迭代法(参见文献[5])求解方程组.式(9)可以改写为 ∑Ee =11Δt κΩeρc T e N i d x d y +k κΩe5N i 5x 5T e 5x +5N i 5y 5T e 5y d x d y + ∫S 3eσεA ({N }T {T}e )4N i d s -∫S 3eσεA T 4∞N i d s -1Δt κΩeρcN i d x d y{T}e t -Δt =0.(10)利用式(5),(10)经整理后变为 ∑Ee =11Δt κΩeρc{N }TN i d x d y{T}e +kκΩe5N i 5x 5N 5xT+5N i 5y 5N 5yTd x d y{T}e - 1Δt κΩeρc{N }T N i d x d y{T}et -Δt -∫S 3eσεA T 4∞N i d s + ∫S 3eσεA ({N }T {T}e )4N i d s =0,i =1,2,…,n.(11)记式(11)为F (T )=∑Ee =1F e(T )=0(其中:F e (T )=[f 1(T ),f 2(T ),f 3(T ),f 4(T )]T;T =[T 1,T 2,T 3,T 4]T ;0=[0,0,0,0]T),则F e (T )的Jacobi 矩阵为 F ′e (T )=5f 15T 1…5f 15T 4………5f 45T 1…5f 45T 4=[K e 1]+[K e 2]+[K e3].(12)其中[K e 1],[K e2]的元素分别为 K e 1ij =κΩek 5N i 5x 5N j x +5N i y 5N j5y d x d y ,K e 2ij =κΩeρcN i N j d x d y.当单元位于上边边界时,[K e3]=∫S 3σεA00000N 220N 4N 20000N 2N 4N 24d s ;22河北师范大学学报(自然科学版)第29卷当单元位于右边边界时,[K e3]=∫S 3σεA 000000000N 23N 4N 30N 4N 3N 24d s.这样F (T )的Jacobi 矩阵为 F ′(T )=∑Ee =1F ′e(T ).利用文献[5]采用下列形式 T k +1=T k +ΔT k ,F ′(T k )ΔT k +F (T k )=0,k =0,1,…,n ,(13)即每一步要解1个n 阶线性方程组.这样利用初始条件(4)认为t -Δt 时刻温度场{T}t -Δt 已知,由式(13)进行迭代求解,给出终止条件‖ΔT k ‖≤ε1‖T k‖(范数‖‖取分量的绝对值的最大值,ε1为事先给定的精度).当‖ΔT k ‖≤ε1‖T k ‖时,取Tk +1=T k +ΔT k作为解.3　计算机模拟利用上述非线性有限元理论,计算板坯从出炉到出炉后30s 时的温度场,并与线性有限元方法计算结果进行比较.计算时取板坯的对称中心为坐标原点,坐标系参见图1.图1　笛卡尔坐标系的选取 由板坯的对称性只需计算第一象限内板坯的温度场,计算时取表1中的参数.板坯右上角的一点(图1中a 点)温度变化最为剧烈,以这点的温度变化为代表对这2种方法的计算结果进行分析.图2～图4中曲线B 和曲线C 分别为a 点用线性方法和非线性方法计算的每间隔5,3,2s 温度变化情况;图5是板坯上边端出炉后第30s 的温度分布,从图2到图5可以看出以下几点:表1　相关的参数Boltzmann 常数/(nW ・m -2・K -2)56.75环境温度/K 298　出炉温度/K 1553　板坯上表面空气温度/K 323　板坯侧面空气温度/K323　板坯高度/mm150　黑度0.8密度/(kg ・m -3)7800　比热容/(J ・kg -1・K -1)670　导热系数/(W ・m -1・K -1)30　板坯宽度/mm 80　计算时间/s 30　预定精度0.1 图2　间隔5s 的温度曲线 图3　间隔3s 的温度曲线 1)非线性有限元方法计算的温度值稍高于用线性方法计算的温度值,并且温度差随时间的间隔减小而减小,实际计算中经常使用线性化的有限元方法,非线性有限元方法在理论上有严格的数学推导,32第1期刘福来等:求解温度场的非线性有限元方法42河北师范大学学报(自然科学版)第29卷 图4　间隔2s的温度曲线 图5　板坯的温度场计算结果与实测温度值较好的吻合;2)板坯表面温度下降幅度大,初始温度差最大达到340K,这是由于板坯尺寸小放热剧烈的结果;3)板坯表面温度分布均匀,只是在边界上与中心点温度差稍大一些,尤其右上角一点与中心点的温度差达到25K,主要是因为边界通过辐射与外界热交换的结果.表2　Newton迭代法计算结果迭代次数kΔT k T k实测温度/K015531213.31-16413892-112.91276.13-54.51221.6表2给出了板坯右上角(a点)在出炉后30s时用Newton迭代法的计算情况,其中ε1=0.1.从表2可知Newton迭代法收敛速度非常快,只需要3次迭代就可以求解,并且求解的温度和实测温度非常吻合.4　结　论本文中,笔者从G alekin有限元方法出发,对自由表面上辐射换热的数学表达式不作线性处理,而对非线性代数方程组,用Newton迭代法计算温度场,得到了令人满意的结果.非线性有限元方法在理论上有严格的数学推导,在计算中构造非线性方程组时有一定的难度.线性化的有限元方法在理论上虽然不太严格,但求解过程比非线性有限元方法简单,而且计算结果与非线性方法的计算结果差别不大,这2种方法计算的轧件的温度曲线基本一致.因此这里进一步验证了线性化的有限元方法求解轧件温度场的实用性.参考文献:[1]　刘相华.刚塑性有限元及其在轧制中的应用[M].北京:冶金工业出版社,1994.[2]　俞昌铭.热传导及其数值分析[M].北京:清华大学出版社,1981.[3]　李长生.棒线连轧轧件温度场有限元解析[D].沈阳:东北大学,1997.[4]　HILL R.New horizons in the machine of solids[J].Mech Ph ys S olids,1956,4(5):66267.[5]　李庆扬,莫孜中,祁中群.非线性代数方程组的迭代解法[M].北京:科学出版社,1997.Nonlinear FEM for C alculating T emperature FieldL IU Fu2lai1,　DU Rui2yan2(1.School of Information Science and Engineering,Northeastern University,Liaoning Shenyang　110004,China;2.Teaching Department,Hebei Y outh Administrative Cadres College,Hebei Shijiazhuang　050031,China)Abstract:Instead of the linearization of expression of thermal radiation on free surface,start at G alerkin finite element method,convert the calculation problem of the temperature field into a calculation of nonlinear equations,and then calculate the temperature by Newton iteration method.K ey w ords:temperature field;finite element method;Newton iteration method(责任编辑　白占立)。