瞬态热传导问题的加权最小二乘无网格法_署恒木

收稿日期:2006-05-08;修改稿收到日期:2008-05-271作者简介:署恒木*
(1957-),男,教授
(E -mail:sh uhm@hd ).
第25卷第6期2008年12月
计算力学学报
C hinese Journal of C omputational Mechanics
V ol.25,N o.6December 2008
文章编号:1007-4708(2008)06-0904-05
瞬态热传导问题的加权最小二乘无网格法
署恒木*1, 黄朝琴2, 李翠伟3
(1.中国石油大学(华东)储运与建筑工程学院,青岛266555;
2.中国石油大学(华东)石油工程学院,青岛266555;
3.上汽通用五菱汽车股份有限公司青岛分公司技术工程部,青岛266555)
摘 要:加权最小二乘无网格法是一种基于节点信息的纯无网格法,该方法使用最小二乘法建立系统的变分原理,通过移动最小二乘法构造近似函数,控制方程在节点处的残量使用最小二乘法予以消除,边界条件通过罚函数法引入。

本文推导了瞬态热传导问题的加权最小二乘无网格计算格式,编制了相应的计算程序,算例结果表明,该方法具有精度高、前后处理简单的优点,是一种高效的的新型无网格法。

关键词:无网格法;移动最小二乘法;权函数;瞬态热传导中图分类号:T K124;O241182 文献标识码:A
1 引言
目前在求解瞬态热传导问题时,有限体积、有限差分和有限单元等网格方法为主要的数值方法,近年来兴起的无网格法[1],因不需要划分网格而受到关注。

无网格法研究始于20世纪70年代,历史最悠久的是光滑质点流体动力学SPH (Smoothed Particle H ydrodynam ics)方法
[2]
;1992年Nayr oles
等提出了漫散元DEM (Diffuse Elem ent M ethod)方法,并首次使用了移动最小二乘法M LS(Mo v -ing Least -Square M ethod)
[3]
;Belytsc -hkohe 改进
了这种方法,并命名为无网格伽辽金方法EFGM (Element -Free Galerkin M ethod)
[4]。

无网格法在我国也得到了不同程度的应用和发展,文献[5-7]将其应用于稳态温度场问题的求解,文献[8]应用EFGM 对瞬态热传导问题进行了求解,但该方法需要背景网格进行积分,并非纯无网格。

张雄等[9~11]提出了加权最小二乘无网格法MWLS(M eshless Weighted Leas-t Square Method),是一种纯无网格法,并应用于稳态热传导问题的求解
[12]。

本文推导了瞬态热传导问题的MWLS 计算
格式,编制了相应的计算程序。

算例结果表明该方法具有精度高、前后处理简单的优点。

2 MWLS 基本理论
M WLS 应用M LS 方法来构造近似函数(或称为试探函数),通过加权最小二乘法来离散方程。

211 MLS 基本原理
设在求解域8内有场函数u(x )及一组随机分布的离散节点x I (I =1,2,,N ),用8i 表示节点x I 的影响域,在二维问题中,影响域常取为圆形或矩形。

使用M LS 进行全局近似,对任意的x I 8有 u h
(x,x -)=
E m
i=1
p i (x -)a i (x )=
p T (x -
)a(x )
(1)
式中 x -=[x ,y ,z ]T 是计算点x 的邻域8x 内各点的空间坐标,u h (x,x -
)或u h (x )是函数u(x -
)的近似表达式;p i (x -)为基函数,m 为基函数的项数,a i (x )相应的系数。

对于二维问题常用的是线性基和二次基。

线性基:p T
(x -)=[1,x ,y ],m =3
(2)
二次基:p T
(x -)=[1,x ,y ,xy ,x 2
,y 2],m =6(3)系数a(x)的选取使得近似函数u h
(x,x -)在计算点x 的邻域8x 内,是待求函数u(x )在某种最小二乘意义下的最佳近似。

设计算点x 的邻域8x 包括n 个节点,近似函数在x -=x I 的误差的加权平方和为 J =
E N I
X I (x )[u h (x,x I )-u(x I )
]2=E
N
I X I (x )
E
m
i=1
p i (x I )a i (x )-u I
2
(4)
式中X(x)为权函数,权函数X I(x)=X(x-x I),只在节点x I周围的一个区域8I中值大于零,而在该区域外为零,即该函数是紧支的。

为了求得系数a(x),令J(x)取最小值,可得
a(x)=A-1(x)B(x)u(5)式中A和B矩阵如下:
A(x)=E N I=1X I(x)p(x I)p T(x I)(6) B(x)=[X1(x)p(x1)X2(x)p(x2),X n(x)p(x n)]
(7)将式(5)代入式(1)中得
u h(x,x-)=N(x,x-)u=N I u I(8)其中形函数:
N(x)=p T(x)A-1(x)B(x)(9)令r=A-1p(10)并对式(9)求导,可得到形函数的一阶和二阶导数为
N k,i=r T,i B+r T B,i(11)
N k,ij=r T,ij B+r T,i B,j+r T,j B,i+r T B,ij(12)
权函数是M LS近似中的重要组成部分,其选择目前还没有理论上的具体规则,带有某种任意性,本文采用C2(8)连续的Gauss权函数:
X(r)=e-r2B2-e-B2
1-e-B2
,r[1
0,r>1
(13)
式中r=+x-x I+/d mI。

影响域对计算的精度和计算效率有直接的关系,本文的影响域半径取为d mI=scale@s[k],其中s[k]为节点x I与距其最近的第k个节点之间的距离,scale是大于或等于1的乘子。

212MWLS基本原理
加权最小二乘法是用来求解具偏微分方程的一种加权残量法,其权函数(或称为检验函数)就是残量本身。

令残量R和R-平方和加权积分:
J R=Q8A R T R d8+Q#B R-T R d#(14) A和B为权系数取极小值,可得
Q8A5R5u j T R d8+Q#B5R-5u j T R d#=0
(j=1,2,,,n)(15)为避免积分,采用如下离散格式:
E r1
k=1
A k5R(x k)
5u j
T
R-(x k)+
E r2
k=1
B k5R
-
(x k)
5u j
T
R-(x k)=0(16)式中A k和B k为权系数,x k为所取节点,r1为计算微分方程残量的配点数,r2为计算边界条件残量的配点数。

若残量的近似函数使用M LS构造,便得到了M WLS方法。

3瞬态热传导MWLS计算格式
311瞬态热传导数学模型
控制微分方程:
k¨2T+q#=Q cT
#
(17)
初始条件:T
t=0
=T0(18)边界条件:T=T-,在#1上(16) n#k¨T=q-,在#2上(17)
n#k¨T=h(T a-T),在#3上(18)式中T
#
=5T/5t,k为热传导系数,q#为热源项,Q 为材料密度,c为材料比热,n为边界的外法线单位矢量,T-为#1边界上的给定温度,q-为#2边界上给定热流值(沿n方向),h为热对流系数,T a为环境温度。

312MWLS计算格式
按照式(14)建立系统方程残量如下:
0=Q8R2(x)d8+Q#1K1R21(x)d#+
Q#2K2R22(x)d#+Q#3K3R23(x)d#=0(19)式中K1,K2和K3为罚函数用来施加边界条件。

对上式求变分,为了避免积分,采用其离散形式,并写成矩阵形式如下:
[K t]T#+[K]T=f(20)式中
[K t]=-E N
s=1
Q ck[¨2N I(x s)]#[N J(x s)](21) [K]=E N s=1[k¨2N I(x s)]#[k¨2N J(x s)]+
E N1
s=1
K1N I(x s)N J(x s)+
905
第6期署恒木,等:瞬态热传导问题的加权最小二乘无网格法
E N
2
s=1K 2
[n #k ¨N I
(x s
)]#[n #k ¨N
J
(x s )]+
E N
3
s=1
{K 3
[n #k ¨N
J
(x s )+h N J (x s )]#
[n #k ¨N I (x s )+h N I (x s )]}
(22)
f =-E
N
s=1
k ¨2
N I (x s )q #
+
E
N
1
s=1K 1N I (x s )T -
+
E N 2
s=1K 2
n #k ¨N I
(x s
)q -
+
E N
3
s=1
K 3
[n #k ¨N
J
(x s )+h N J (x s )]hT a (23)T =[T 1
T 2
,T N ]
T
(24)
为了保证式(20)中各项的数量级相当,罚函
数可取为
K 2=105
,K 1=K 2k 2
,K 3=K
2min 1,k hL 2
L 为问题的特征长度。

对于稳态问题:
[K]T =f
(26)
对于瞬态问题,进一步求得时间向后差分格式:[K ]+1$t [K t ]T (t n +1)=f (t n+1)+1$t
[K t ]T (t n )
(27)
上式即为瞬态热传导问题的M WLS 计算格式。

4 算例分析
411 稳态算例
如图1所示,中心开孔环形板:外圆半径为R b
=100cm ,内孔半径为R a =011cm ,外圆给定恒温T b =20e
,内圆孔给定温度T a =8e 。

由对称性可首先考虑任意一条半径上的布点,在此条半径上
图1 圆环平板示意图
Fig.1 Annulu s plate s cheme
采用对数等分方法布置节点,第i 个节点位置为
r(i)=exp ln R a +
ln R b -ln R a n -1
i
(28)
式中n 为节点数,取n =40(i =1,2,,,n)对于圆盘,在H I [0,2P ]上进行20等分。

本算例理论解为式(29)。

图2所示为沿半径r 方向计算结果的温度曲线分布图。

节点数和影响域半径对计算精度影响较大,本例中B =210,s cale =310,计算误差范数为0122%,基本与理论解一
致,说明了MWLS 方法的可靠性。

T (r ,H )=T b -T b -T a ln R b R a
ln R b
r (29)
误差范数定义为
L T =
E N
s =1
(T
h
s
-T s )
2
E N s=1
(T
s
)2
@100%(30)
式中T h s 和T s 分别为节点x s 的近似值和精确值。

412 瞬态算例
考虑某110m @110m @011m 混凝土块,已知混凝土的热传导系数k =101836kJ/(m #h #K),容积比热Q c =252kJ/(m 3#K),平板厚度为011m 。

初始条件:T |
t=0
=400x (1-x )
边界条件:T |x =0=T |x=1=0e
采用均匀规则布点,如图3所示,应用M WLS 进行求解,权函数参数B =210,d m I =317@s
[4],时间步长$t =014h(hour)。

本文求解了[0,214h ]时间区间上各时刻的温度分布。

图2 温度沿半径分布图
Fig.2 T emperatu re distribution along radius
906
计算力学学报
第25卷
图3 混凝土块及其节点布置
Fig.3 Con crete plate and its n od es
distribution
图5 误差范数比较
Fig.5 The comparison of er ror norm
本算例在x 方向上的理论解见式(31)。

图4为各时刻沿x 方向的温度分布图。

图5为M WLS 与有限元(FEM )的计算误差范数比较,其中,有限元采用三角形二次单元进行求解,网格剖分如图6所示。

求解软件为ANSYS,在求解时间上两者相当。

由此可知M WLS 的计算精度和效率较有限元高。

T(x ,t)=
E
]
n=1-800
(n P )3
(n P sin n P +2cos n P -2)#e -(n P a)2t
sin n P x
(31)
式中a =k/(Q c)。

5 结 论
本文应用M WLS 对热传导问题进行了求解,该方法不需要划分网格,只需要节点信息。

算例结果表明M WLS 具有前后处理简单、计算量小、精度高的优点,
为进一步研究复杂温度场问题提供了基础。

M WLS 是一种高效的纯无网格方法。

图4 各时刻温度沿x 轴分布曲线
Fig.4 T he tem
perature distrib ution alongxaxis
图6 FEM 网格剖分Fig.6 FEM g rids
参考文献(Referneces):
[1] BELY T SCHK O T ,K RO N GA U Z Y ,O RGA N D,et
al.M eshless methods:A n o verv iew and recent deve-l opments[J].Comp ut M ethods A p p l M ech Engr g ,1996,139:3-47.
[2] L U CY L B.A numerical approach t o t he testing of
the fissio n hypothesis[J].T he A s tr on ,1977,8(12):1013-1024.
[3] Nay ro les,T o uzot,V illo n.G ener alizing the f inite ele -ment met ho d:diffuse approx imat ion and diffuse ele -ments[J].Comp M ech ,1992,10:307-318.
[4] BEL Y T SCHK O T ,L IU Y Y,G U L.Element -free
G aler kin methods[J].I nter national J our nal f or N u -mer ical M ethods in Engineer ing ,1994,37:229-256.[5] 陈莘莘,李庆华.用无单元法求解稳态热传导问题
[J].株洲工学院学报,2003,17(5):86-88.
[6] 苑素玲,葛永庆,王章奇.无网格法在温度场中的应
用[J].华北电力大学学报,2003,30(2):82-86.[7] 李庆华,陈莘莘.用加权最小二乘无网格法求解稳态
热传导问题[J].株洲工学院学报,2005,19(4):71-
907
第6期
署恒木,等:瞬态热传导问题的加权最小二乘无网格法
73.
[8]高志华,曾辉辉,等.无单元伽辽金法及其在瞬态温度
场中的应用研[J].冰川冻土,27(4):557-562.
[9]张雄,胡炜,潘小飞,等.加权最小二乘无网格法
[J].力学学报,2003,35(4):425-431.
[10]张雄,刘岩.无网格法[M].清华大学出版社,
2005.[11]ZH A NG X,SO N G K Z,LU M W.M eshless meth-
od based on collo cat ion with radial basis functio ns [J].Comp ut M ech,2000,26(4):333-343.
[12]L IU Y,ZH A N G X,L U M W.M eshless least-
squar es method for so lv ing the steady-state heat con-
ductio n equat ion[J].T s ing hua Science and T echnol-og y,2005,10(1):61-66.
Meshless weighted least-square method for transient
heat conduction problems
SH U H eng-mu*1,H UANG Zhao-qin2,LI Cu-i w ei3
(1.Co llege of Sto rag e&T r anspor tatio n and Ar chitect ur al Engineer ing,China U niver sity o f Petr oleum,
Q ingdao266555,China;
2.Colleg e of P et ro leum Eng ineer ing,China U niv ersity of P et roleum,Q ingdao266555,China;
3.SA IC G M Wuling A utomobile Co.,L td Q ing dao Branch Co mpany T echnolog y&Eng ineer ing D epar tment,
Q ingdao266555,China)
Abstract:M eshless Weig hted Least-Square M ethod(M WLS),w hich is only based on no rma-tion.In MWLS the variatio n principle of system is constructed w ith w eighted least-square residual meth-o d,and the approx im ation function is built by mo ving least-square metho d.T he boundary conditio ns ar e enforced w ith a penalty function.Appling M WLS to so lve transient heat conduction pro blems,the nu-m er ical results show that M WLS has sever al advantag es such as hig h accuracy and easy for pre-process and post-process.
Key words:meshless method;mov ing least-square method;w eig ht functio n;transient heat co nductio n
908计算力学学报第25卷。