求解温度场的非线性有限元方法

Ξ

求解温度场的非线性有限元方法

刘福来1,　杜瑞燕2

(1.东北大学信息科学与工程学院,辽宁沈阳　110004;2.河北青年干部管理学院教务处,河北石家庄　050031)

摘要:从G alerkin 有限元方法出发,对自由表面上的辐射换热的数学表达式不作线性化处理,而是把温

度场的求解问题转化为非线性代数方程组的求解问题,并且用Newton 迭代法计算了温度场.

关键词:温度场;有限元方法;Newton 迭代法

中图分类号:O 242.21 文献标识码:A 文章编号:100025854(2005)0120021204

由文献[1]知,求解二维待轧过程的温度场,就是要求下面微分方程和初值问题的解:

52T 5

x 2+52T 5y

2=1α5T

5t ;(1) -k 5T

5n

=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 n

5y 2

ωi d x d y =0,i =1,2,…,n.

(6)

对式(6)应用Green 公式,有

Ξ收稿日期:2004

0105;修回日期:20040420

作者简介:刘福来(1975),男,河北省唐山市人,东北大学博士研究生.

第29卷第1期2005年　1月河北师范大学学报(自然科学版)

Journal of Hebei Normal University (Natural Science Edition )Vol.29No.1Jan.2005

κ

Ω

ρc

5T n 5t ωi +k 5ωi 5x 5T n 5x +5ωi 5

y 5T n

5y

d x d y -

∫

S 3

k

5T n

n ωi

d s =0,i =1,2,…,n.(7)

取基函数ω1,ω2,…,ωn 为形函数,即求{T}e ,使得T e ={N }T {N }e 满足

∑

E

e =1κ

Ω

e

ρc 5T e

t

N i +k 5N i 5x 5T e 5x +

5N i 5y 5T e 5y

d x d y -

∫

S 3

k

5T n

n ωi

d s =0,i =1,2,…,n.(8)将边界条件代入式(8)且对5T

e 5t

作差分后有下面的等式:

∑E e =1κ

Ω

e

ρc T e -T e

t -Δt Δt N i +k 5N i 5x 5T e 5x +5N i 5y 5T e

5y 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)可以改写为

∑E

e =11

Δt κ

Ωe

ρc T e N i d x d y +k κ

Ω

e

5N 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)经整理后变为

∑E

e =1

1

Δ

t κ

Ω

e

ρc{N }

T

N i d x d y{T}e +k

κ

Ω

e

5N i 5x 5N 5x

T

+

5N i 5y 5N 5y

T

d x d y{T}

e -

1Δt κ

Ω

e

ρc{N }T N i d x d y{T}e

t -Δt -∫

S 3e

σεA T 4

∞N i d s +

∫

S 3e

σ

εA ({N }T {

T}e )4

N i d s =0,i =1,2,…,n.(11)

记式(11)为F (T )=

∑E

e =1

F 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 1

5T 4

………5f 45T 1…

5f 4

5T 4

=[K e 1]+[K e 2]+[K e

3].(12)

其中[K e 1],[K e

2]的元素分别为

K e 1ij =κ

Ωe

k 5N i 5x 5N j x +5N i y 5N j

5y d x d y ,K e 2ij =

κ

Ω

e

ρcN i N j d x d y.当单元位于上边边界时,[K e

3]=

∫

S 3

σεA

00

00

0N 2

2

0N 4N 2

00

00

N 2N 4

N 2

4

d s ;

2

2河北师范大学学报(自然科学版)第29卷