第五章：大变形问题的有限单元法

第五章 大变形问题的有限单元法
1. 弹性大变形问题的有限元法
2. 弹性分支点稳定问题有限元分析
3. 物质描述大变形增量问题的T.L 、 U.L法
大变形时平衡方程和虚位移原理
变形体初始和现时位形 如图所示，以欧拉应力表 述平衡时
V : ij , x j Fi 0 S : i ij n j
e Ve 0
1. 弹性大变形问题的有限元法
弹性大变形问题，需要考虑变形的非线性项 和变形对平衡的影响。 若以初始自然平衡状态作初始位形，则物质 描述的格林应变为 ui uk uk 1 u j L E ij ( ) ij ij 2 X i X j X i X j 式中 ui 1 u j 1 uk uk L ij ( ) ij 2 X i X j 2 X i X j 线性部分 非线性部分
x i ui ij X j X j
x i V: ( S kj ) Fi 0 0 X j X k
对比小变形情况，可见大变形时变形对平衡的 影响，是通过变形或位移梯度表现出来的。
ui V: [ S kj ( ik )] Fi 0 0 X j X k
为便于计算机编程，将张量转换为矩阵： 格林应变矩阵和张量的分量间有如下关系 T E E11 E22 E33 2E12 2E23 2E31 T 11 22 33 2 12 2 23 2 31
2 2 2 对应的克希荷夫应力矩阵和张量分量间关系为 T S S11 S22 S33 S12 S23 S31 引入两个算子矩阵 T d,X1 0 0 d,X 2 0 d,X 3 A 0 d,X 2 0 d,X1 d,X 3 0 0 0 d,X 3 0 d,X 2 d,X1
S E W ij ij
因此克希荷夫应力和格林应变在能量上共轭。
由于变形率和欧拉应力张量是对称的，因此
再利用欧拉和拉格朗日应力间关系，可得
v i w ijVij ij x j
ki
ij J
1
x j X k
W ij
v j X i
ij
i i i i ij ij V S V
柯西应变
为建立物质描述虚功方程，先讨论能量共轭 关系。 空间描述中 变形率张量是 w ij ij ijVij 相对现时位形 定义的柯西应 单位体积变形功率 变的速率 又因为 x i x j x k x l 1 ij J S kl E ij Vkl X k X l X i X j
L L 11 L 22 L 33 L 12 L 23 L T 31



A
L
u ,X 1 0 0
0 u,X 2 0
0 0 u,X 3
u ,X 2 u,X 1 0
0 u,X 3 u,X 2
u,X 3 0 u,X 1
T
式中
u T d,X i u,X i u u1 u2 u3 Xi Xi 3阶单位矩阵 再引入位移梯度向量的记号 T ui u I3 I3 I3 X j X1 X2 X3
由(AL)可见，格林应变-位移关系是非线性的。 为用虚位移原理建立单元刚度方程，还需有 应变增量和位移增量间的关系。对线性部分 δ ε Bδ e 非线性部分因为 E Au ( 1 A H )u L L 2 A θ 2A θ
同理可验证，对任意j恒有
X k (J )0 X k x j
ij X k ji V :J Fi 0 Fi 0 0 S : ji N j i 0 X k x j X j
利用克希荷夫应力和拉格朗日应力的关系，可 X j 将平衡条件改为 S
空间描述的虚位移原理
设现时位形微小虚 位移在V内单值连续、 在位移边界上为零。 则外力总虚功为
w e Fi ui dV
V
i ui dS
S
考虑到位移边界虚位移为零和应力边界条件
we Fi ui dV ij n j ui dS
V S
将平衡方程引入，考虑到虚位移微小，则
G e 所以 1 L L L L ε ( A θ ) A θ A Gδ e 2 1 L 综上所述，格林应变增量为 L B A G L 2 E ( B A G)δe ˆ ˆ 如果记 B ( B 2B L ) ，则 E Bδe 。 ~ ~ L E B e B e B e B B B L
u j X i
因此拉格朗日应力和初始位形位移梯度在能量上共轭。
由此即可得到物质描述的虚功方程为
Fi 0ui dV0 i 0ui dS0 S ijE ij dV0
V0 S 0 V0
V0
Fi 0ui dV0 i 0ui dS 0 ij ( X )dV0
we
V
ij x j
ui dV ij n j ui dS
S
利用格林公式，可得
虚 we ij ij dV w变 wd 功 V 率 也即虚位移原理的虚功方程为 方 Fi ui dV i ui dS ij ij dV 程 V S V 若虚位移用虚速度、虚应变用虚应变率代替， 则 F v dV v dS dV
对弹性问题，在物质描述下本构关系(克希荷 夫应力和格林应变)为 o Skl Dklpq E pq X p X q X k X l o Dklpq J Dmnij x m x n x i x j
在上述基础上，由虚位移原理可得 T T ˆ e e T S Ed V0 S Bδ e d V0 ( F j FE ) δ e 式中 F 为单元结点力矩阵，F 为单元等效结 点荷载矩阵 FEe N T F0 d V0 N Tφ d S 0
证明
x1 ,1 J x 2 ,1 x 3 ,1 x1 , 2 x2 ,2 x3 ,2 x1 , 3 x 2 , 3 e ijk x i , 1 x j , 2 x k , 3 x3 ,3
X k 1 xi xl J ekmn e jil x j 2 X m X n
以。 在外荷为保守力系 时
Fi 0 dV0 Fi dV i 0dS0 i dS
由复合函数求导数，平衡方程改写为
上式乘以J，由于 平衡方程成为
V : ij , X k X k , x j Fi 0
J dV dV 0
X i ij J lj xl
e j
e E
V0
V0e
e S
按集成规则集装后可得 ˆ T Sd V P P R B 0 d E
e Ve 0
ˆ E Bδe
Pd和PE分别为直接和等效结点荷载矩阵，R为 综合等效结点荷载矩阵。上式也可写为 ˆ (U ) B T Sd V0 R 0
X k 1 x 3 x 2 x 2 x 3 j 1 J e kmn ( ) x1 2 X m X n X m X n
X k x2 x3 x3 x2 (J ) ( ) X k x1 X 1 X 2 X 3 X 2 X 3 x2 x3 x3 x2 x2 x3 x3 x2 ( ) ( )0 X 2 X 3 X 1 X 3 X 1 X 3 X 1 X 2 X 1 X 2
H u
u N e
G e
式中AL为如下6×9的矩阵
A
L
1T 0 0 0 2T 0 0 T 0 3 T 2 1T 0 T T 0 3 2 T 0 T 1 3
u i Xi
式中G为如下9×3m的矩阵
N1 I3 X1 N1 G HN I3 X2 N1 I3 X 3
Nm N2 I3 I3 X1 X1 Nm N2 I3 I3 X2 X2 Nm N2 I3 I3 X3 X3
ij X k X k V :J Fi 0 ( J ij ) X k x j X k x j 对任意j恒有 X k ij (J ) Fi 0 0 X k X k x j (J )0 经推证得 X k x j
ij X k ji V :J Fi 0 Fi 0 0 S : ji N j i 0 X k x j X j
H
X 1 H u 0 0
0 X 1 0
0 0 X 1
X 2
0 X 2 0 X 2
X 3
X 3
0
u X 3
T
在上述符号基础上，格林应变有位移表为 1 L L E Au ( A H )u 2 设单元位移场为 u N e 则单元格林应变为 ~ ~ L L E B e B e B e B B B 其中线性应变矩阵B和线性分析一样 B AN 非线性部分“应变矩阵”为 1 L 1 L L B A HN A G G HN 2 2
x i V: ( S kj ) Fi 0 0 X j X k x i S : S jk N k i0 X j
ij
x k
ik
上述拉格朗日和克希荷夫应力表示的平衡条件 都是以初始位形作参考的物质描述。
如果考虑到变形梯度和位移梯度间的关系，
则克希荷夫应力表达的平衡方程
可改为
S 0 V0 i
u j
由虚位移原理可得
V0
ˆ δ d V ( F e F e ) T δ e S Ed V0 e S B e 0 j E e
T T
式中 F 为单元结点力矩阵，F 为单元等效结 点荷载矩阵e T T FE N F0 d V0 N φ d S 0
L L
θ 1T 0 0 L A θ T θ 2 0 θ T 3
0
T θ 2
0
θ 1T T θ 3
0
T T θ1 θ1 θ1 θ1 0 T T θ 2 θ 2 θ 2 θ 2 0 θ1 T T T θ3 θ3 θ3 θ3 θ 3 2 2 AL θ θ2 T T T T 0 θ 2 θ1 θ1 θ 2 θ 2 θ1 θ1 θ 2 θ T T T T T 3 θ 2 θ3 θ 2 θ 2 θ3 θ3 θ 2 θ 2 θ3 θ T θ θ T θ θ Tθ θ Tθ θ 1T 1 3 3 1 3 1 1 3