基于余能原理计算简单剪切变形

第8卷 第16期 2008年8月1671-1819(2008) 16-4445-05

科 学 技 术 与 工 程

Science T echno logy and Eng i neering

l 8 N o . 16 A ug . 2008 V o.

Z 2008 Sc. i T ech . Engng.

力 学

基于余能原理计算简单剪切变形

范志会 金 明 尚新春

(北京科技大学土木工程学院, 北京100083; 北京交通大学工程力学所, 北京100044)

摘 要 以高玉臣提出的基于基面力的弹性大变形余能原理为基础, 利用L agrange 乘子, 放松平衡条件和力边界条件对余能泛函的约束, 给出了基于基面力的广义余能泛函。利用极分解定理, 将表示单位质量的余应变能分为刚性转动和纯变形两部分, 考虑相邻单元间边界面力的平衡条件, 推导出了平面四边形单元的有限元公式, 并利用其进行了简单剪切变形的计算。数值计算结果与理论结果的对比表明, 当剪切角不超过35b 时, 二者是非常接近的。关键词 余能原理 有限元 简单剪切变形中图法分类号 O 34; 文献标志码

A

简单剪切变形作为一种简单的普适变形, 许多学者从理论上进行了研究

[1, 2]

的定义如下:I 1=tr (C) , I 2=

*

, 只要给出势函数, 就

可以确定应力分布。然而, 到目前为止, 还未发现用余能原理研究该变形的有关文献。本文首先利用线弹性本构关系给出了可压缩弹性体应力的理论解, 接下来以高玉臣提出的基于基面力的弹性大变形余能原理为基础, 利用Lagrange 乘子放松约束条件对余能泛函的约束, 给出了广义的余能泛函及相应的有限元公式, 最后将简单剪切变形简化为平面应变问题, 应用本文给出的平面四边形单元进行了有限元计算, 并与理论解进行了对比。

122*

[

I 1-tr (C ) ],I 3=det (C) 2

(2)

取初始的Lagrange 坐标系与空间坐标系为相同的直角坐标系, 令长方体初始构型的各面分别和直角坐标系的坐标面平行, 如图1所示, e 1, e 2, e 3分别为x, y , z 方向的单位矢量。

1 简单剪切变形的理论解

根据文献[1, 2],可压缩弹性体的本构方程为R =

C +

J 1

*

*

I 2

*

52*52*-1

I -I 3C *+I 3*

5I 25I 25I 3

(1)

(1) 式中I 1, I 2, I 3是Cauchy 变形张量C 的不变量, 2是弹性势, 它是不变量I 1, I , I 的函数。不变量

2008年4月23日收到

高等学校博士学科点专项科研基金

项目([1**********]) 资助

*

2

*3

图1 简单剪切变形

根据图1, 可知式(1) 的分量形式为

R

11

52=52+252+15I 25I 3

(3)

4446科 学 技 术 与 工 程

2

8卷

R

22

=(1+K )

5225252+(2+K ) *+*15I 25I 3

(4)

udV 0-i

S 0u

t u dS -Q (t - t ) udS Q

S 0R

[5]

(13)

(13) 式中T 为基面力

i

, u 为位移, W c 为单位质量

R

33

5225252

=+(2+K ) *+*

15I 25I 3

R

12

(5) (6) (7)

的余应变能, x (i =1, 2, 3) 为随体坐标, V P 为初始的物质坐标系下由协变基矢量组成的平行六面体的体积, Q 0, t 0分别为弹性体变形前的密度及单位面积上的边界面力, t 0为变形前的力边界上给定单位面积的面力, dV 0, dS 0分别是初始的物质坐标系下弹性体的体元和边界面元。

对任意的多面体单元, 假定作用在各面的应力矢量是均匀分布的, 那么单元应力写为

R =

A

T ªr A V e

[6]

=2K R

23

525+*15I 2

31

=R =0

式中K =tan C 。根据以上各式, 给定2的形式, 就能完全确定弹性体的应力状态。

对各向同性体, 仅保留应变分量的二次幂项, 其线弹性本构关系

2=

8

[3]

可以用不变量表示为

2

2K +2G I 1-K +G I 1-43*I 2+K +G 283

(8)

(14)

(14) 式中V e 是该单元的体积, 指标A 表示单元体的边界面, T 是该面上的力矢量, r A 是当前构形上原点到该面几何中心的位置矢量。

利用极分解定理, 将余应变能分为表示纯变形和刚性转动的两部分

(9) (10) (11) (12)

[5]

A

(8) 式中K 为拉梅系数, G 为剪切模量。那么, 应力与变形的关系为

R

R

22

11

2

=K

+G 2

。根据各向同性材料的线弹性

本构关系及散度定理, 可写出单元的余能表达式(Wc ) =

e

=K

2

K 2

(1+K ) +(2+K )G

2

2

A B A M A 2

(TT ) r AB -(Tr A ) +T u r A

2EV e 1+M

(15)

A

R

R

12

33

=

K 2

K 2

(15) 式中u A 为刚性转动在T 作用处引起的位移,

r

=K

22

+(1+K )G 2

E 为杨氏模量, M 为泊松比, 并且r AB =r A r B 。

不考虑体积力的作用, 利用单元的平衡条件, 弹性体总的广义余能泛函为P c (T, u A , u 0, u A ) =

r

A

从以上各式可以看出, 应力的各分量值除了与弹性常数K 和G 有关, 还与表示变形的K 有关, 而与点的位置无关, 故对简单剪切变形来说, 它的变形是均匀的, 所有各点的应力状态是完全相同的。并且, 正应力分量是K 的偶函数, 而剪应力分量是K 的奇函数。

E

e

A

A B

(TT ) r A r B -2E V e

r

M A 2

(Tr A ) 1+M

S u

e

+T u A +

E

A

T u 0-0

A

u t d S -Q (t - t ) u d S Q

S R

e

(16)

2 基于广义余能原理的有限元模型

本节使用E i n ste i n 求和约定。根据高玉臣

[4]

(16) 式中T , u A 分别是单元体A 面上的力矢量及其作用处的位移矢量, u 0是该单元几何中心的位移。

上述泛函将单元的平衡条件和力边界条件进行了放松, 此外, 在相邻单元间的边界上还应满足边界面力的平衡条件, 当所有单元间的公共边界都被考虑时, 弹性体的广义余能泛函可写为

A

A

提出的弹性大变形余能原理, 利

用Lagrange 乘子, 放松平衡条件和力边界条件对余能泛函的约束, 可得到广义的余能泛函P c (T, u )=

i

V

Q 0W c P c (T, u A , u 0, u A ) r

16期范志会, 等:基于余能方法计算简单剪切变形4447

M A 2

(Tr A ) 1+M

A

+T u A +

r A

A

E

A

T u 0-(17)

A

T u A +T u A

A

A

(17) 式中T 为单元的力边界上给定的力矢量, 若A 面不是给定的力边界, 则T =0。假设发生的纯变形为小变形, 则变形后的位置矢量r A 可近似看作是变形前的位置矢量经过刚性转动得到的, 那么有r A r B =常数。泛函P c 取驻值时, 得

E

e

(C A B #T +u A +u 0-u A ) =0

r

B

(18) (19)

图2 平面四边形单元

E

e

(T

? A

-T )

A

S

R

=0

R =

co s H si n -sin H co s 0

(25)

各个单元还应满足

根据以上各式, 显然有

E T

A A

A

=0(20) (21)

[6]

x A =x A cos H +y A sin H y A =-x A sin H +y A cos H r A #r B =r A #r B =x A x B +y A y B r A ªr B =x A x B e 1ªe 1+x A y B e 1ªe 2+

y A x B e 2ªe 1+y A y B e 2ªe 2

矩阵C A B

(26) (27) (28) (29)

T #D u A =0

r

(18) 式中C AB 为单元内任意两点之间的柔度矩阵

C A B =

r A B I -r A ªr B

E V e 1+M

(22)

(22) 式中I 为单位二阶张量。

根据式(22), 还可得到该单元任意两点之间的柔度

00

1 00M x A x B x A y B

C A B =(xA x B +y A y B ) -EA 0 1+y A x B y A y 3 平面四边形单元

本节取消脚标的求和约定。

对平面问题, 取初始的Lag range 坐标系和平面坐标系为相同的直角坐标系, e 1, e 2分别为x , y 方向的单位向量。

假设作用在各边界上的应力是均匀分布的, 则可认为各边上的力矢量均作用在其中点位置, 将其取为节点进行编号, 其单元节点编号顺序如图2所示。节点在初始时刻的坐标用(x A , y A ) (A =i , j , k, l ) 表示, 当单元绕e =e 3逆时针刚性转过H 后, 相应的位置坐标变为(x A , y A ) , 那么, 该单元上各边中点在初始时刻和刚性转动后的位置矢量分别为

(30)

(30) 式的每一个分量都是位置坐标和转角H 的函数, 式中A 为平面四边形单元的面积。

由刚性转动在各边中点引起的位移为u A =r A -r A =[x A (co s H -1) +y A si n H ]e 1+

r

[-x A sin H +y A (cos H -1) ]e 2

00

(31)

r

显然, 在转轴已知的情况下, 对给定的物质点, u A 只是刚性转角H 的函数。

若作用在节点上的力分别用它在x ,

y 方向的分量X A 和Y A 表示, 根据图2所示, 显然有

T =p A X A e 1+p A Y A e 2

A

(32)

r A =x A e 1+y A e 2r A

=x A e 1+y A e 2=R r A

式中R

为转动张量, 在二维情况下

(23) (24)

(32) 式中

p A =

1, -1,

A =j , k A =i , l

(33)

4448科 学 技 术 与 工 程

A

8卷

T D u A =

r

E

A

[pA X A (-x A sin H +y A cos H ) +

00

定铰链支座, 左端为滑动铰链支座, 令其发生简单剪切变形。

(34) (35) (36)

根据结构特点, 将该问题简化为平面应变问题。由于其变形是均匀的, 故只取一个单元进行研究, 其节点的总体编号如图3所示。表1给出了单

元编号与总体编号之间的关系。

p A X A (-x A cos H -y A sin H ) ]D H

u 0=u 01e 1+u 02e 2u A =u A 1e 1+u A 2e 2

的有限元公式

(34) 式即为单元在刚性转动后位置的力矩表达式。令

根据式(18) ) 式(21), 可以得到该平面四边形单元C ii (1, 1) p i X i +C ii (1, 2)p i Y i +C ij (1, 1)p j X j +C ij (1, 2) p j Y j +C ik (1, 1) p k X k +C ik (1, 2) p k Y k +C il (1, 1) p l X l +C il (1, 2)p l Y l +x i (cos H -1) +

y i si n H +u 01-u i 1]=0

(37)

图3 横截面为矩形的结构表1 单元编号与总体编号的关系

单元编号

节点i 1

节点j 3

节点k 4

节点l 2

C ii (2, 1) p i X i +C ii (2, 2)p i Y i +C ij (2, 1)p j X j +C ij (2, 2) p j Y j +C ik (2, 1) p k X k +C ik (2, 2) p k Y k +

C il (2, 1) p l X l +C il (2, 2) p l Y l -x i sin H +

y i (co s H -1) +u 02-u i 2]=0

s

-X i +X j +X k -X l =0-Y i +Y j +Y k -Y l =0

A =i , j , k, l

(38) (39) (40)

总体编号

采用增量法, 利用式(14) 确定应力增量。取结构的长为1m, 宽为0. 5m, E =210M Pa , M =0. 3, 当K 取不同数值时, 所得各应力分量的值与理论值的对比见图4。

(41)

E

[p A X A (-x A sin H +y A co s H ) +

00

p A Y A (-x A cos H -y A sin H ) ]=0

对弹性体而言, 式(37) 和式(38) 需要考虑全部单元进行整合, 式(39) ) 式(41) 反映的是单元的平衡方程, 每一单元均需分别给出。显然, 上述各式是关于T , u A , u 0和H 的非线性方程组。

若将X A 和Y A 看作是由于载荷增量引起的内力增量, 令

$T

A 1A

=p A ?X A , $T

A 2

=p A Y A (A =i , j , k, l )

(42)

将其代入式(18) ) 式(21), 就可以得到形如式(37) ) 式(41) 的平面四边形单元增量形式的有限元公式。

图4 简单剪切变形的结果对比

图中各点是应用广义余能原理计算所得的数值结果, 各条线是根据式(9), 式(10) 和式(12) 画出的理论曲线, 其中实线、短划线和点线分别是R 11, R 22和R 12的理论曲线。

从图4可以看出, 当K 达到0. 7时, 此时剪切角b , 4 基于广义余能原理的简单剪切变形有限元计算

如图3所示横截面为矩形的结构, 其右端为固

16期范志会, 等:基于余能方法计算简单剪切变形4449

好的吻合。当剪切角继续增大, 数值计算结果与理论曲线之间就会逐渐产生偏差, 这可能是由以下几个方面的原因造成的:

(1) 在推导增量形式的有限元公式时做了近似处理, 即有限元公式本身存在误差。其中任意两点间的柔度矩阵是在刚性转动后的构形上给出的, 而实际应在变形后的构形上进行计算。

(2) 若给定外载荷, 通常认为在实际加载过程中载荷的大小和方向不变, 但是构件发生简单剪切变形时, 其各面上不仅有正应力的存在, 还有剪应力, 这些应力分量随着构件的变形, 其大小和方向是不断发生变化的, 故加载过程存在误差。若给定位移条件, 则不能保证构件精确满足简单剪切变形的位移模式, 同样也存在误差。

(3) 数值计算过程中存在误差。

极分解定理, 把表示单位质量的余应变能分为刚性转动和纯变形两部分, 得到了基于广义余能原理的有限元公式。利用平面四边形单元进行了简单剪切变形的数值计算, 结果表明, 该有限元模型简单, 计算简便, 当剪切角不超过35b 时, 数值结果与理论解是非常接近的, 故应用余能方法进行该变形的计算是可行的。

致谢 北京交通大学工程力学所高玉臣院士生前对作者给予了悉心指导和大力帮助, 在此深表感谢。

参 考 文 献

1 郭仲衡. 非线性弹性理论. 北京:科学出版社. 19802 黄筑平. 连续介质力学基础. 北京:高等教育出版社. 20033 吕和祥, 蔡和洋. 非线性有限元. 北京:化学工业出版社. 19924 G ao Y C. Co m p l e m en tary energy princi p l e for l arge elasti c def or m a -tion , Science i n Ch i na :series G Phys i cs . M echan ics &A stron o my ,

5 结论

根据高玉臣提出的基于基面力的弹性大变形余能原理, 推导出了无条件的广义余能原理。利用

2006, 49(3); 341) 356

5 高玉臣. 固体力学基础. 北京:中国铁道出版社, 1999

6 Gao Y C . A ne w descri p tion of t he stress state at a poi n tw it h app lica -tion . Archive of Appli ed M echan i cs , 2003, 73, 171) 183

Calculation of Sim ple Shear D eformati on Based on the

Co mplem entary Energy Pri nci ple

FAN Zh-i hu, i JIN M i n g , S HANG X i n -chun

(Univers i ty of S ci en ce and Technol ogy Beiji ng , C i vil&Env i ronm en t Engi n eeri ng Schoo, l Beijing 100083, P 1R 1Ch i na ;

Instit u t e of Eng i neeri ng M echan ics , B eiji ng Jiaotong Un i versity 1, Be iji ng 100044, P 1R 1Ch i na)

1

[Abstract] U si n g t h e co mp le m entary energy pri n ciple for large elastic ity proposed by G ao Yuchen based on t h e base forces and Lagrange mu lti p liers , the genera lized co mp le m entary ener gy f u nctional based on base forces w as ob -tained by re leasi n g constraint cond iti o ns o f equ ili b ri u m equations and force boundary cond itions . According to po lar deco m position, the co mp le m entary energy o f un itm ass could be deco m posed i n to t w o parts :ri g id rotation and pure defor m ation parts . The finite ele m en t for m u l a e for plane quadrilatera l e le m entw ere ga i n ed w ith the consi d eration of equ ilibri u m conditions bet w een adjacent ele m ents , w hich w as used to the co m putation for si m ple shear defor m a ti o n . Co m parison sho w s t h at the num erical resu lts are very closed to theoretical solutions when the shear ang le is not ex -ceed 35b .

[ co le m energy i n p le le m si p le defo r