第六章杆单元的有限变形理论及算法(徐春晖、李明瑞)

端的外力P平衡。可以把等式左边的SA0项称之为杆的Kirchhoff内力NK，有NK = S
A0。于是可得到杆的Cauchy内力Nσ与杆的Kirchhoff内力NK的转换关系，以及与
外力的平衡关系：
Nσ = NK( Lt/L0) = NK(1+ u, x)= P 。
(6.3)
再利用本构关系 S = E ε，可以得到简单拉伸时的有限变形理论解为：
及 P = NΣU −3 = eU −3 = 1 (1 − U −2 )U −3 。显然，选
E1A0 E1A0
2
取不同的应变度量，结果会有很大的不同。这是因为我们在选取 Lagrange 应变ε
的同时又假设了其共轭应力 PK2 应力 S 与ε之间满足线性本构关系 S=Eε。
类似，在选用 Almansi 应变 e 为应变度量时，也假定了线性本构关系 Σ = E1e 。
第六章 杆单元的有限变形理论及有限元算法
§1 杆的简单拉伸
取杆单元的应变为 Lagrange 应变，即
ε = u, x + (u, x)2/2 = (Lt2 – L02)/(2L02)。
(6.1)
其中，L0为杆的原长，Lt为变形后杆的现长。
取计算应力为与 Lagrange 应变共轭的第二 Piola-Kirchhoff 应力（PK2），S。并设
2 L2t
2 L2t 2
1
的 例 3 ， 可 知 其 共 轭 的 应 力 为 Σ = TU 2 = SU 4 = σU 2J = σU 3 At 。 所 以 令 A0
P=σ At ; NΣ = Σ A0 则有 NΣU −3 = P 。如果假设有本构关系 Σ = E1e 。类似于（6.3），
（6.4）我们有 NΣ = Σ A0;
u,x = u0/L0；v,x = v0/L0；w,x = w0/L0。 其他的位移偏导数为零。变形梯度矩阵为：
（6.5)
⎡1 + u, x 0 0⎤
[ ] F
=
⎢ ⎢
v, x
1 0⎥⎥ = Λx Λy Λz ；
(6.6)
⎢⎣ w, x 0 1⎥⎦
以下均采用Lagrange应变张量作为应变的度量。由（4.16），Lagrange应变张量的
At S11
A0 L0 2 At Lt 2
ΛxΛxTΛx
= NK
L0 2 Lt 2
ΛxΛxTΛx
=
NK Λx 。
（6.8)
这是因为有：Kirchhoff轴力NK = A0S11。 注意，Cauchy内力向量Tn 是沿着变形后的杆轴方向的，Tn在初始位形坐标系上 的投影分量为：
σ
=
1 J
Λx S11Λx T 。
(6.7)
由图（1）所示的变形前后杆的方向，可见杆变形后横截面上的法向量为
[ ] n = 1
Lt
L0 + u0
v0
w0
=
L0 Lt
Λx 。
由于有关系： Λx T Λx
=
Lt 2 L0 2
，所以该法向量是一个单位向量。
于是可得杆的 Cauchy 内力向量为
Tn
=
At σ ⋅ n =
对于同一种材料的的杆件而言，这两组本构关系是互相冲突不可能同时成立的。
选用不同应变度量，而采用类似的本构关系，表面上看都是选用了同一个材料常
数 E，本质上却反映了选用不同本构关系。所以有限变形问题的计算结果往往会 有有较大差异。实际上许多材料在弹性大变形时，载荷位移的实验曲线就是非线 性的，但是即使载荷位移的实验曲线完全是线性的，在选用一定的共轭应力应变 对后，它们间的本构关系材料系数也一定不是线性的。这方面的研究还很不够， 缺乏可资参考的资料。
分量可以写成 Eij
=
1 2
(
Λi
Λj
− 1) 。各个分量中只有E11
=
(ΛxTΛx
–
1)/2 非零外，其
他全为零。于是相应的PK2 应力分量也只有S11 = E E11 ≠ 0 , 其他全为零。
由 Cauchy 应力张量与 PK2 应力张量的转换关系式：
σ = 1 FSF T ， J
其中J=V/V0= AtLt/A0L0，得杆的Cauchy应力张量为
作为一个练习，读者可以试用对数应变来定义杆的应变，仿照上述方法可以
得出 P = U −1 lnU EA0 可能有的读者会提出这样一个问题，对于同一个问题，选用不同的应变度量，
是否计算结果一定不同。答案则是否定的。这个结论与上面的论述并不矛盾。关 键就在于如何确定本构关系。例如说，先选定应变度量 1 和本构关系 1 可以得到 第一组结果。如果另选定应变度量 2，也希望得到同样的计算结果，那么，应该 选择本构关系 2 与本构关系 1 满足某种张量变换关系。特别当本构关系 1 可以用 常数表达时，那么本构关系 2 就一定不是常数。在本章的第三节就将给出这样的 例子。 §2．杆件有限变形的一般情况
变形梯度矩阵也可简化为F =(1+ u,x)= Lt/L0=U。其中，A为杆的截面积。下标t与
0 分别表示变形后与变形前所在的时刻。于是应力转换关系简化为：
S =(AtLt)/(A0L0)( Lt/L0)-2σ。
或Hale Waihona Puke SA0( Lt/L0) = SA0 U =σAt。
(6.2)
注意到等式右端部分即为杆的当前内力，或杆的Cauchy内力Nσ，它将与施加在杆
本构关系为 ：S = E ε 。假设弹性模量 E 对初始位形不变。变形后，即当前时刻
的 Cauchy 应力张量σ与 PK2 应力张量 S 的转换关系为：
σ = FSFT/J。
杆是一维问题，所以张量记号均改换成标量记号。应力张量 6 个分量和应变张量
6 个分量中只有沿杆轴方向的一个分量不为零，其他全为零。且有：J=(AtLt)/(A0L0),
y
B’
n
Lt
v0
A(A’) L0B
B u0B
B
O
2
x
图 6.1 变形前后二维杆方向示意图
如图 6.1 所示设AB和 A’B’分别为变形前、后的杆件长度和位置，不失一般性，
可以假设A点不动，在三维空间中B点的位移为(u, v, w)。取变形前的杆轴方向为x
轴。杆变形前和后的横截面积分别为A0与A1。此时有：
P/EA0
=
NK EA 0
(1 +
u,X
)
=
S E
(1 +
u,X
)
= (1+
u,
x)ε
=
= 1U (U 2 − 1) 2
(1+ u, x)[ u, x + (u, x)2/2] （6.4)
这个结果与 Poisson 比无关。
如果选用 Almansi 应变，e = 1 ( L2t − L20 ) = 1 (1 − L20 ) = 1 (1 − U −2 ) ,由上一章第 7 节