基于六面体体积坐标的新型8结点实体单元

清华大学学报(自然科学版) 2009年第49卷第11期

　

CN 11-2223/N J T singh ua Un iv (Sci &Tech ) , 2009, V o l. 49, N o. 11w 29

http://qhx bw. chinajo urnal. net. cn

　

基于六面体体积坐标的新型8结点实体单元

李宏光, 　岑　松, 　岑章志

1, 2

1, 3

1, 3

(1. 清华大学航天航空学院, 北京100084; 2. 空军航空大学, 长春130022; 3. 清华大学破坏力学教育部重点实验室, 北京100084)

摘　要:六面体体积坐标方法是构造高性能三维实体单元的新工具。基于六面体体积坐标方法, 构造了含有内参的8结点实体单元HV 3D 8。其基本位移场的形函数由点协调广义协调条件精确确定, 并按照Wilso n 非协调元的模式进一步建立了单元内部位移场, 这样使得该单元的位移场对整体坐标是二次完备的。数值算例表明:该单元在各种弯曲问题中不仅计算精度高, 而且抗网格畸变能力优于其他同类等参元, 显示了六面体体积坐标和广义协调理论相结合的特有优点。

关键词:有限元; 六面体体积坐标; 广义协调; 网格畸变; 8

结点三维实体单元

中图分类号:O 242. 21

文章编号:1000-0054(2009) 11-1856-05

文献标识码:A

20世纪60年代T aig 和Ir ons 提出的等参坐标方法现在已经被广泛地用于构造有限元模型。然而在一般情况下, 等参坐标与整体坐标之间并不

[3]

是线性关系, Lee 和Bathe 的研究表明:当单元发生畸变时, 一个用等参坐标表达的完备的高次插值函数仅仅等价于整体坐标表达的一阶完备的插值函数。因此当有限元网格发生畸变时, 数值模拟的结果通常和真实解相去甚远。

为了克服二维问题中等参坐标的上述缺点, 龙驭球和龙志飞等建立了第一类四边形面积坐标方法, 并成功地应用于对构造网格畸变不敏感的四边形平面膜元

[6-8]

[4-5]

[1][2]

、轴对称元、薄板元

[9][10]

等模型;

陈晓明和岑松等

[11]

又建立了更为简便的第二类面

8-node solid element based on the

hexahedral volume coordinate method

LI Hon gguang 1, 2, CE N S on g 1, 3, CEN Zhangzhi 1, 3(1. S chool of Aerospace , T s inghua University ,

Beij ing 100084, China ;

2. Aviation University of Air Force , Changchun 130022, China ; 3. Failure Mechanics Laboratory of the Ministry of Education ,

T s inghua University , Beij ing 100084, China )

Abstract :Th e h exahedral volum e coordinate method is a new tool for developin g h igh-perform ance 3D solid finite element m od els. An 8-node h exahedral elemen t HV3D8containing internal parameters w as s ucces sfully formulated us ing th e hex ahedral volume coordin ate meth od. T he sh ape function s for th e fun damental dis placemen t field of the elem ent are exactly determined by the n od al-version gen eralized confor ming conditions, wh ile the internal dis placemen t field is constructed u sing th e W ils on ’s n on conforming elemen t meth od. T hu s, the entire dis placem ent field for the elemen t poss es ses s econ d-or der completeness in the global coordinates. Numerical res ults indicate th at th e elemen t exh ibits better per formance than other con ventional is oparametr ic models for variou s b ending problems w ith m esh dis tortion. Th e res ults dem onstrate efficiency of th e HVCM and th e generalized conforming theory for developing simple, effective and reliab le finite elements. Key words :finite element; hexah edral volum e coordin ate;

ing; mes h dis 8-node 3D

s 积坐标方法, 目前也在构造高性能平面膜元方面有

[12]

成功应用。

六面体单元是在三维有限元分析中常见的单元, 传统的8结点六面体等参单元在弯曲等高阶问题中常常出现各类闭锁问题。为了构造出高性能的8结点六面体单元, 国内外许多研究者都进行了相关的研究, 并取得了一些进展。例如, 基于Wilson

[14]

[13]

和T ay lor 等提出的非协调元方法的模型, 基于

[15]

Hu -Washizu 变分原理的多变量模型, 基于增强假设应变(EA S) 方法的模型

[17]

[16]

, 基于杂交元方法的

模型等, 都在一定程度上提高了单元抗网格畸变能力。然而, 由于以上这些单元都是采用等参坐标方法, 因此高阶问题中的单元网格畸变敏感问题并没有得到根治。

收稿日期:2008-10-06

基金项目:国家自然科学基金资助项目(10872108) ;

高等学校全国优秀博士论文作者专项基金资助项目(200242) ;

教育部新世纪优秀人才支持计划项目(NCET -07-0477)

作者简介:李宏光(1973—) , 男(汉) , 吉林, 博士研究生。:

李宏光, 等:　基于六面体体积坐标的新型8结点实体单元

1857

借鉴四边形面积坐标方法在二维问题中的成功

[18]

应用, 李宏光和岑松等提出了针对三维单元的六面体体积坐标方法。该方法把六面体内任一点与单元表面形成的四棱锥体积与单元总体积的比值作为该点的坐标分量。由于这种体积坐标与整体坐标之间的线性关系, 使得构造对网格畸变不敏感的新型六面体单元成为可能。同时, 这种坐标也可以用等参坐标来表示, 因此也可以应用于各类复杂形状。本文

[19]

将六面体体积坐标与广义协调理论相结合, 构造出一种新型的六面体单元HV3D8。

　　对应的结点1—8的体积坐标分别为:

(0, g 12, 0, g 14, g 15, 0) , (g 21, 0, 0, g 24, g 25, 0) ,

(g 31, 0, g 33, 0, g 35, 0) , (0, g 42, g 43, 0, g 45, 0) , (0, g 52, 0, g 54, 0, g 56) , (g 61, 0, 0, g 64, 0, g 66) , (g 71, 0, g 73, 0, 0, g 76) , (0, g 82, g 83, 0, 0, g 86). (1) 其中g ij (i =1, 2, …, 8; j =1, 2, …, 6) 是单元的特征参数, 其定义由文[18]给出。

假设基本位移场采用六面体体积坐标表达为

u =A 1+A 2(L 2-L 1) +A 3(L 4-L 3) +A 4(L 6-L 5) +A 5(L 2-L 1) (L 4-L 3) +

A 6(L 4-L 3) (L 6-L 5) +A 7(L 2-L 1) (L 6-L 5) +A 8(L 2-L 1) (L 4-L 3) (L 6-L 5) .

1　基本位移场

每个表面都是平面四边形的六面体单元如图1所示。

(2)

式中:u 表示u 方向的基本位移分量, 因为u 、v 、w 具有相同的形式, 这里只列出u ; L 1～L 6是该点的六面体体积坐标分量, 其定义由文[18]给出。为了确定式中的待定常数A 1～A 8, 引入广义协调的点协调条件可以得到

(u -~u ) j =0, 　(j =1, 2, …, 8) . (3)

其中~u j 表示单元结点的位移。将式(1) 和(2) 代入式(3) 得到

图1　六面体单元的结点编号

A A =~u .

其中:

-g 12g 14-g 21g 24-g 31g 33-g 42g 43-g 52g 54-g 61g 64-g 71g 73-g 82g 83

-g 14g 15-g 24g 25-g 33g 35-g 43g 45-g 54g 56-g 64g 66-g 73g 76-g 83g 86

-g 15g 12-g 25g 21-g 35g 31-g 45g 42-g 56g 52-g 66g 61-g 76g 71-g 86g 82

T

(4)

1-g 121-g 211-g 31

A =

1-g 421-g 521-g 611-g 711-g 82

-g 14-g 15-g 24-g 25-g 33-g 35-g 43-g 45-g 54-g 56-g 64-g 66-g 73-g 76-g 83-g 86

-g 12g 14g 15-g 21g 24g 25-g 31g 33g 35-g 42g 43g 45-g 52g 54g 56-g 61g 64g 66-g 71g 73g 76-g 82g 83g ,

A =[A 1　A 2　A 3　A 4　A 5　A 6　A 7　A 8], T T

~u =[~u 1　~u 2　~u 3　~u 4　~u 5　~u 6　~u 7　~u 8]=[u 1　u 2　u 3　u 4　u 5　u 6　u 7　u 8].

因此待定参数可以由下式确定:

-1

A =A ~u . 　　再将式(5) 代入式(2) 可得

8

i

a 7i (L 6-L 5) (L 2-L 1) +

(5)

a 8i (L 2-L 1) (L 4-L 3) (L 6-L 5) ,

(i =1, 2, …, 8).

-1

(7)

u =

0i

∑N

i =1

u i . (6)

其中a 1i ～a 8i 是矩阵A 的元素。这里的形函数是具有精确封闭解的。

　　其中形函数N 可表示为

N i =a 1i +a 2i (L 2-L 1) +a 3i (L 4-L 3) +a 4i (L 6-L 5) +a 5i (L 2-L 1) (L 4-L 3) +

6i L 32　内参位移场

为使位移场达到二次完备, 在基本位移场的基

1858

清华大学学报(自然科学版) 2009, 49(11)

u K =K 1L 1L 2+K 2L 3L 4+K 3L 5L 6. (8)

B K =[B K 1B K 29N K i

00

B K i =

9N K i 09N K i e

B K 3]; (17)

其中K 1、K 2、K 3为单元内部自由度。显然, 内参位移场满足点协调条件, 且内参形函数为:

N K 1=L 1L 2, 　N K 2=L 3L 4, 　N K 3=L 5L 6. (9) 当单元退化为平行六面体时, 它们退化为三维Wil-son 单元的内参形函数。

因为六面体体积坐标与整体坐标之间是线性关系, 所以单元的位移模式具有整体坐标的二次完备性。

09N K i 09N K i 9N K i 0

009N K i 9z 09N K i 9N K i (18)

, 　(i =1, 2, 3) .

3　单元刚度矩阵

基于上述基本位移场和内参位移场的假设, 可以得到单元的总体位移场

u v =w

e

u v w

u K

+

v K =N q q +N K K . w e

e

e

对内参K 进行凝聚, 最后得到单元刚度矩阵

K =K qq -K q K K K K K K q .

其中:

K qq =K K K

e

-1

00

(10)

(19)

其中:q 是由结点位移表示的单元自由度矢量, K 是由单元内部自由度矢量,

N q =

N 100

0N 10

00N 1

N 200

0N 20

00

…N 8…

00

0N 80

00, N (11)

K K q

∫=∫B DB d V , =∫B DB d V .

T V e V

e

B q DB q d V ,

T

K T K

K

(20)

V e

q

N 2…

这里D 是三维弹性系数矩阵。通过以上过程得到的单元记为HV3D8。单元列表见表1。

表1　单元列表

单元名称C 3D 8C3D8R C 3D 8I HV 3D 8

简要说明

A BAQ U S 中的8结点三线性六面体单元[20]A BAQ U S 中的采用减缩积分和沙漏控制的8结点六面体单元[20]

A BAQ U S 中的8结点非协调六面体单元[20]本文的8结点六面体单元

N K =

N K 100N K 200…N K 8000N K 100N K 20…0N K 80. 00N K 100N K …00N K 2(12)

　　对位移场进行微分, 可以得到单元应变

e e e E =B q q +B K K . (13)

其中应变和应变矩阵分别为:

e T E =[E (14) x 　E y 　E z 　2C xy 　2C yz 　2C z x ];

B q =[B 19N

9x 009N i 0i

000i

4　数值算例

例1　分片检验

由于HV 3D 8单元采用的协调条件是点协调, 因此对于图2所示的有限尺寸常应力分片检验不能

B 20

…B 8];

(15)

9N 09N i 09N i 0

00i

09N i 09N i 0i 00

B =

0i

, 　(i =1, 2, …, 8) ;

李宏光, 等:　基于六面体体积坐标的新型8结点实体单元

1859

表3　梁承受横向载荷P 时端部挠度

挠度

网格图3a

/(°) H [**************]75

C 3D 80. 1000. 0730. 0500. 030—0. 0540. 0340. 022—

C 3D 8R 0. 9930. 8330. 7860. 797—0. 2670. 2180. 199—

C 3D 8I 0. 9930. 8320. 7850. 794—0. 2640. 2140. 194—

HV 3D 80. 9920. 9920. 9930. 9951. 0060. 9930. 9930. 9961. 009

给出精确解, 但是当网格细分时该单元可以收敛到精确解, 因此可以通过弱式分片检验。

例2　含畸变网格的MacNeal 细长梁如图3所示, 在不同网格划分情况下悬臂梁分别承受弯矩M 和面内横向载荷P , 计算结果用以检验单元的抗网格畸变性能。弹性模量E =107M Pa, Poisson 比L =0. 3, 倾角为H

。

图3b

图3c

　　注:精确解为1. 000b , b 标准值为0. 1081。

5　结　论

本文基于六面体体积坐标方法和广义协调理论, 构造了一个8结点六面体单元。数值算例表明该单元收敛可靠, 计算精度高, 特别是具有很强的抗网

格畸变能力, 从而证实了六面体体积坐标与广义协调理论相结合是一种构造高性能单元的有效手段。

参考文献　(References )

[1][2]

T aig I C. Structu ral analysis by th e matrix dis placement m ethod [R]. English Electric Aviation Report, 1961:S 017. Irons B M. Engineering application of nu merical integration in s tiffness method [J]. J A merican Institute of Aer onautics and A str onautics , 1966, 14:2035-2037. [3]

Lee N S, Bathe K J. Effects of elem ent distortion on th e performan ce of isoparametric elem ents [J ]. I nt J Nume ric al M e thods in Eng ineer ing , 1993, 36:3553-3576. [4]

龙驭球, 李聚轩, 龙志飞, 等. 四边形单元面积坐标理论[J ]. 工程力学, 1997, 14(3) :1-11. LONG

Yuqiu,

LI

Juxu an,

LONG

Zhifei,

et

al. [J].

Area-coordinate theor y for quadrilateral elemen ts [5]

图3　MacNeal 梁弯曲表2　梁承受弯矩M 时端部挠度

挠度

网格图3a

H /(°)

C3D 8

[1**********]

图3c

456075

0. 1000. 0710. 0470. 026—0. 0480. 0280. 016—

C3D 8R 1. 0000. 8750. 8390. 849—0. 2070. 1590. 138—

C3D8I 1. 0000. 8740. 8380. 847—0. 2050. 1570. 136—

HV 3D81. 0001. 0001. 0001. 0001. 0001. 0001. 0001. 0001. 000

图3b

E ngineering M ec hanics , 1997, 14(3) :1-11. (in Chinese) 龙志飞, 李聚轩, 岑松, 等. 四边形单元面积坐标的微分和积分公式[J]. 工程力学, 1997, 14(3) :12-20.

LONG Zhifei, LI Juxuan, CEN Song, et al. Differential and in tegral formulas for area coordinates in quadrilateral elemen ts [J]. Eng inee ring M echanic s , 1997, 14(3) :12-20. (in Chinese) [6]

龙志飞, 陈晓明. 采用面积坐标的四结点四边形膜元[J]. 中国矿业大学学报, 2000, 29(3) :310-313. LONG Zhifei,

CHEN Xiaoming. Four-node quadrilateral

plane elemen t usin g area coor dinates [J]. J Ch ina Univ ersity of M ining &T echnology , 2000, 29(3) :310-313. (in Ch ines e) [7]

, . [J]. 清华大学学() (　　注:精确解为1. 000a , a 标准值为0. 0054。

　　从表2和表3的结果可以看出, 对于畸变网格, 其他单元的计算精度都发生不同程度的下降, 甚至不能完成计算。而本文单元HV 3D 8却显示了很好

1860

清华大学学报(自然科学版) 2009, 49(11)

LONG Yuqiu , CHE N Xiaoming. T wo robus t qu adrilateral membrane elem ents [J]. J T sing hua Univ (Sc i &Te ch ) , 2003, 43(10) :1380-1385. (in Chinese) [8]

陈晓明, 岑松, 宋德坡. 采用面积坐标的抗畸变四边形曲边膜元[J ]. 清华大学学报(自然科学版) , 2007, 47(2) :248-255, 259.

CHEN Xiaoming , CE N Song, SONG Depo. Curved-side membrane elements ins ens itive to mes h distortion usin g quad rilateral area coordinates [J]. J T sing hua Univ (S ci &T ech ) , 2007, 47(2) :248-255, 259. (in Chin es e) [9]

管楠祥, 岑松, 陈晓明. 基于四边形面积坐标的广义协调轴对称单元[J ]. 工程力学, 2007, 24(Su p II) :161-167. GUAN Nanxiang,

CEN S ong,

CHEN

Xiaoming.

New

gen eralized confor ming axis ymmetr ic elemen ts b ased on the quad rilateral ar ea coordinate m ethod [J]. E ng ineer ing M echanics , 2007, 24(Sup II) :161-167. (in Ch ines e) [10]陈晓明, 岑松, 龙驭球. 采用面积坐标和基于假设转角的薄

板元[J ]. 工程力学, 2005, 22(4) :1-5, 30.

CHEN Xiaoming, CEN Song, LONG Yuqiu. Tw o th in plate elements

developed b y ass uming

rotations

and

usin g

quad rilateral area coor dinates [J]. Eng ine ering M echanics , 2005, 22(4) :1-5, 30. (in Chinese)

[11]陈晓明, 岑松, 龙驭球, 等. 含两个分量的四边形单元面积

坐标理论[J ]. 工程力学, 2007, 24(Su p I) :32-35. CHEN Xiaom ing, CE N Song , LONG Yuqiu,

et al.

A

two-componen t area coor dinate meth od for qu adrilateral elements [J]. E ngineering M ec hanics , 2007, 24(Sup I) :32-35. (in Chines e)

[12]陈晓明, 岑松. 基于四边形面积坐标法的平面单元解析试函

数法[J ]. 清华大学学报(自然科学版) , 2008, 48(2) :289-293.

CHEN Xiaoming,

CEN S ong.

Analytical trial function

[13]W ils on E L, T ayler R L, Doh erty W P, et al. Incom patible

dis placemen t models [M ]//Fen ves S J, et al. Numerical and Com putational M ethods in S tr uctural M ech an ics. New York:Academic Press, 1973:43-57.

[14]T aylor R L, Beresford P J, W ilson E L. A non -conforming

elemen t for stres s analysis [J ]. I nt J N umerical M e thods in E ngineering , 1976, 10:1211-1219.

[15]Cao Y P, Hu N , Lu J, et al. A 3D b rick elemen t b as ed on

Hu -Washiz u variational principle for mesh distortion [J]. I nt J N umerical M ethods in Eng ineer ing , 2002, 53:2529-2548.

[16]S imo J C, Rifai M S. A class of mixed assu med strain

m ethods and the m ethod of incompatible m od es [J]. I nt J N umer ical M ethods in E ngineering , 1990, 29:1595-1638. [17]Ch eun g Y K, Chen W J. Isoparam etric h ybrid hexah edral

elemen ts for three dimensional s tres s analysis [J ]. I nt J

N umer ical M ethods in Eng ine ering , 1988, 26:677-693. [18]李宏光, 岑松, 龙驭球, 等. 六面体单元体积坐标方法[J].

工程力学, 2008, 25(10) :12-18.

LI Hongguang, CEN S on g, LONG Yuqiu, et al. Volum e coor dinate m ethod for hexahedr al elem ents [J]. Eng ineering M e chanics , 2008, 25(10) :12-18. (in Chinese) [19]Lon g Y Q,

Xu Y.

Generalized conformin g triangular

m embr ane element with vertex r igid rotational freedom [J]. F inite Elements in A naly sis and Desig n , 1994, 17:259-271.

[20]ABAQU S Inc. ABAQUS Docu mentation Version 6. 5[R].

ABAQU S Inc, Rawtucket, Rhode Island, 2004.

meth od for plan e elem ents bas ed on quadrilateral ar ea coordinate theory [J]. J T sing hua Univ (S ci &Te ch ) , 2008, 48(2) :289-293. (in C hinese)