13-2.静力隐式-显式和动力显式有限元列式

弹塑性大变形有限元方法
胡平
§大变形弹塑性本构方程
最简单的材料大变形弹塑性本构方程是不考虑加载历史的形变理论的全量本构方程。

这类本构方程基于加载历史是比例加载的基本假定。

对于许多金属冲压成形问题，比例加载的基本假定是近似正确的，因此，基于这类理论的所谓one step inverse algorithm，由于其高效高速的计算效率，近年来被广泛应用于汽车概念设计阶段的冲压工艺性快速校核。

关于全量理论的One-step inverse algorithm的理论知识将另行介绍（参见理论手册）。

然而，更精细和准确的成形问题数值分析应该是能够考虑加载历史的增量本构理论。

为了反映与加载历史的相关性，大变形弹塑性问题需要采用速率型(增量型)的本构方程。

首先，需要在大变形条件下，按照符合客观性的要求，建立增量型本构方程。

由(7.26)式可知，变形率张量[d]是客观张量，但由(7.36)式可知Cauchy应力张量的物质导数[]σ 不是客观张量，所以若用变形率张量和Cauchy应力张量来建立本构方程，则将不满足本构方程的客观性条件(7.42)式。

为此，需要另外定义Cauchy应力张量的导数。

1．Cauchy应力张量的Jaumann导数（率）
而Cauchy应力的Jaumann导数定义为
01lim tj
dt dt
(1)
推导后得到
ij
ij ik kj kj ki σσσωσω∇=--
（2）
其矩阵形式为
2. 第一Piola 应力（名义应力）张量的本构导数
（3）
用Cauchy 应力Jaumann 导数表示的第一Piola 应力的本构导数为
()ij ij ik kj kj ki ik jk ij kk t t d d l l σσσσσ∇
=--++ （4）
其矩阵形式为
()][t t σ∇
= （5）
3. 第二Piola 应力张量的本构导数
仿照第一Piola 应力本构导数的推导，有
()ij ij ij kk kj ik ik jk T t l l l σσσσ=+--
（6）
把（7.179）式代入（7.182）式，便得用Cauchy 应力Jaumann 导数表示的第二Piola 应力的本构导数为
()ij ij ik kj kj ki ij kk T t d d l σσσσ•
∇=--+
（7）
其矩阵形式为
（8）
4. 大变形弹塑形本构方程
三维大变形弹塑形问题属于大变形问题。

对于大应变问题，
Pr Re andtl uss -流动理论关于应变增量等于弹性与塑性应变增量之和的
假定一般是不成立的。

但对于金属材料的大应变问题，弹性应变比塑性应变小得多，所以仍然可以近似的用Pr Re andtl uss -增量流动理论。

大变形弹塑性加——卸载准则与小变形弹塑性加——卸载准则在形式上相同，只要把其中的应力认为是Cauchy 应力即可。

大变形弹塑形速率型（或增量型）本构方程也与小变形弹塑性增量本构方程在形式上类似，只要把其中的应力和应变分别改用Cauchy 应力的Jaumann 导数和变形率，
应力偏量和等效应力分别改用Cauchy 应力偏量和等效的Cauchy 应力便可。

本构方程的张量分量形式为
ep ij ijkl kl D d σ∇=
（9）
若令
{}11
22331223
31
T σσ
σσσσσ∇
∇∇
∇∇∇∇⎡⎤= ⎣⎦
（10）
{}[]
112233122331222T
d d d d d d d =
（11）
则对各向同性材料有矩阵形式，即
{}{}ep D d σ∇⎡⎤=⎣⎦
（12）
考虑到弹性变形微小，塑性变形时体积不变，（由（4）式便得
()ep
ij ijkl kl ik kj kj ki ik jk t t D d d d l σσσ=--+
（13）
()ep
ij ijkl kl ik kj kj ki T t D d d d σσ=--
（14）
大变形弹塑性静力隐式和静力显式有限元法
对于弹塑性大变形有限元而言，也象非线性弹性大变形问题一样，可以有TL 和UL 两种算法。

这取决于本构方程用何种应力应变去描述。

一般而言，能够用第二Poila 应力和Green 应变表示的本构方程，可以用TL 或UL 算法建立有限元方程。

1. 第一种作法（静力隐式UL 列式理论基础）
（1） 增量形式本构方程
把（8）式和（12）式联立可得
()()()()e p
ij pq kj ki i jpq
ik kj T t D
E t E t E t σσ•
•••
=--
（15）
由于其中e p i j pq D 和ij σ与时间t 无关，所以在上式两端乘以dt ，并写成有限增量形式为
e p ij i j
pq
pq ik kj kj ki T D E E E σσ ∆=∆-∆-∆
（16）
若把上式展开，并考虑到ij σ、ij E ∆的对称性和材料的各向同性，则可写成矩阵形式，即
{}[]()
{}ep d T D E σ⎡⎤∆=-∆⎣⎦ （17）
其中
[]11121233200000020d σσσσσσσσ312223=
0 0 2 ()()2331121211223123
233122331231312311102221112221102σσσσσσσσσσσσσσσσσ23 + 0 + ()1211331
22σσσ⎛⎫ ⎪ ⎪ ⎪ ⎪
⎪ ⎪ ⎪
⎪
⎪ ⎪ +⎪⎝⎭
（18） (2) 单元增量刚度方程
由于t t +∆时刻的单元应变和单元应力为
{}{}E E =∆
（19） {}{}{}T T σ=+∆ （20）
所以把上述各式式虚功方程，有
{}{}[]()
{}{}
(
){}{}{}{}{}T
ep d e
T
e e
T
T
e
a E D E dv u f u p dv u p da
σ
δσσδδδ⎡⎤⎡⎤∆+-∆⎣⎦⎣
⎦
=∆+∆+∆⎰⎰⎰ （21）
于是仿照几何非线性三维大变形有限元法的U.L 法，便得单元增量刚度方程，即
[][][](){}
{}{}{}0
e
e
e
e
e e e
l
k k k u f q r σ++∆=+- （22）
其中小位移弹塑性刚度矩阵
[]
[][]0
e
T
L ep L e
k B D B dv ⎡⎤=⎣⎦⎰
（23）
初应力刚度矩阵
[][][][][][][](
)
e
T T
L d L e
k G G B B dv σσσ=-⎰
（24）
大位移弹塑性刚度矩阵
[]
[][]()[]()
[][]()
****
T T L ep d N N ep d L e
l T e N ep d N B D B B D B k dv B D B σσσ⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤-+-⎣⎦⎣⎦⎣⎦⎣⎦⎢⎥=⎢⎥⎡⎤⎡⎤⎡⎤+-⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦
⎰ （25） 初应力节点力
{}[]{}T
e
L e r B dv σ=⎰ （26）
2. 第二种作法（静力显式有限元理论基础）
（1）持续平衡方程
把虚功率方程（5.233）式用于t 时刻和t dt τ=+时刻两个充分临近的构形，有
()()()
T
T
j
Ij I i I i I A j
I j I j I I i I I i I A v t d P v d P v dA X v t dt d P dP v d P d P v dA X δδδδδδΩΩ Ω
Ω∂⎛⎫
Ω=Ω+ ⎪∂⎝⎭∂⎛⎫
+Ω=+Ω++
⎪∂⎝⎭⎰⎰⎰⎰⎰⎰ 上两式相减后可得
T
j Ij I i I i I A v dt d dP v d d P v dA X δδδΩ
Ω∂⎛⎫
Ω=Ω+ ⎪∂⎝⎭⎰⎰⎰ 若把其中的Ij dt 、I dP 和I d P 表示成物质导数，则得持续平衡方程，即
T
j ij i i I I I A v t d p v d p v dA X δδδ••
•Ω
Ω∂⎛⎫
Ω=Ω+ ⎪∂⎝⎭⎰⎰⎰ （2）单元持续平衡方程
由于把t 时刻构形作为参考构形，仿照（7.175）式的推导，t τ→时，有
j ji I
v l X ∂→∂，()()I i p p t τ••
→，()()I i p p t τ••
→，所以把上面持续平衡方程
用于单元e ，并考虑t 时刻单元节点力的时间变化率()
()()1,2,...,8k i f t k •=，可得单元持续平衡方程，即
()()
()()
()()1
ij ji
e
a
k k i
i
i i i i i e
a t t l
dv
f t v p t v dv p t v da
σ
δδδδ==++⎰∑⎰⎰ （27）
再把（4）式代入上式，并考虑到体积不变，略去含kk l 的项，便得
()()
()()
1
ij ik kj kj ki ik jk ji e
k a
k i
i i
i
i i
k e
a d d l l dv f t v p v dv p v da
σ
σσσσδδδδ∇
•
••
=--+=++
⎰∑⎰⎰ （28）
由于
,,ij ji ij ji ij ji d d σσσσ∇∇
===
适当更换哑标，有
ij ji ij
ij l d σδσδ∇
∇= （29）
()()ik
kj kj ki ji ij ik kj d d l d d σ
σδσδ+=
（30）
()1
2
ik jk ji ij kj ki l l l l σδσδ=
（31）
于是，把上述各式代入（28）式，便得
()()()
8
1
122ij ij ij ik kj kj ki e k k i i i i i k e
a d d d l l dv f v p v dv p v da
σ
σδσδδδδ∇•
=⎡⎤
--⎢⎥⎣
⎦=++⎰∑⎰⎰ （32）
（3）单元变形率与速度梯度
若以节点速度为基本未知量，并令
{}()()()()()()111888123123T
e
v v v v v v v ⎡⎤= ⋯ ⎣
⎦ （33） {}[]111213212223313233T
l l l l l l l l l l =
（34）
单元变形率和单元速度梯度为
{}[]{}
e
L d B v =
（35） {}[]{}e
v l B v =
（36）
其中
[][][]v v B L N =
（37）
[]123123
123000000000000000000v x x x L x x x x x x ⎛⎫∂∂∂ ⎪∂∂∂ ⎪ ⎪∂∂∂
= ⎪∂∂∂ ⎪ ⎪∂∂∂ ∂∂∂⎝
⎭T
⎪⎪
（38）
而[]L B 与线性弹性有限元法一样。

（4）单元应力本构导数 单元应力本构导数为
{}[]{}e
ep
L
D B v σ∇
⎡⎤=⎣⎦
（39）
（5）单元速率刚度方程
由于
{}{}{}
(
)
[][]{}T
ij ij e
e
T
T
e e L ep L e d dv d dv
v B D B dv v σδδσδ∇∇
=⎛⎫⎡⎤ = ⎪⎣⎦⎝⎭
⎰⎰⎰ （40） (){}[]{}{}
(
)
[][][]{}T
ij ik ki d e
e
T
T
e e L d L e d d dv d d dv
v B B dv v σδδσδσ=⎛⎫ = ⎪⎝⎭
⎰⎰⎰ （41）
(){}[]{}{}
(
)
[][][]{}12T
ij kj ki v e
e T
T
e e v v v e l l dv l l dv v B B dv v σδδσδσ=⎛⎫ = ⎪⎝⎭
⎰⎰⎰ （42） 其中
[][][][][][][][][][]000000v σσσσ⎡⎤
⎢⎥= ⎢⎥⎢⎥ ⎣⎦
（43）
[]111231122223312333σσσσσσσσσσ ⎡⎤⎢⎥= ⎢⎥⎢⎥ ⎣⎦
（44）
()()
{}(
){}1e
w
T
e k k i i
k f v v f δδ•
•
==∑
（45）
{}
()()()()()()111888123123T
e
f
f f f f f f ••••••
•
⎡⎤= ⋯ ⎢⎥⎣⎦
（46）
{}()[]
{}
T
T
e i
i
e
e
p v dv v N p dv δδ•
•=⎰⎰ （47）
{}()[]
T
T
e i
i
a a p v da v N p da σ
σ
δδ•
•⎧⎫
=⎨⎬⎩⎭
⎰⎰ （48）
各式代入（40）式，并考虑到{}(
)T
e v δ的任意性，便得单元速率刚度方程，
即
[][](){}
{}{}
e
e
e
e
e
e
k k v f
p p σ
•
•
•
⎧⎫
+=++⎨⎬⎩⎭
（49）
其中
[][]
[]0e T
L ep L e
k B D B dv ⎡⎤=⎣⎦⎰
（50）
[][][][][][][](
)
e
T
T
v v v L d L e
k B B B B dv σσσ=-⎰
（51） {}
[]
{}
T
e
p N p dv ••
=⎰
（52）
[]e
T
a p N p da σ
••⎧⎫⎧⎫
=⎨⎬⎨⎬⎩⎭⎩⎭
⎰ （53）
（6）总体速率刚度方程
仿照线性弹性有限元法，可由单元速率刚度方程组集成总体速率刚度方程，即
[][](){}{}
K K V p σ•
+=
（54）
其中
{}1
111
2
3
1
2
3
...T n n
n V v v v v
v v ⎡⎤= ⎣⎦
（55） [][]()[][]
01
m
T
e e
e
n
e K A k A ==
∑
（56） [][]
()
[][]1
m T
e e
e
e K A k A σσ==∑
（57） {}{}{}
p F Q •
•
•
=+
（58）
{}
[]
()
{}
1
e e m
T
e e Q A p p ••
•=⎛⎫
⎧⎫ ⎪=+⎨⎬ ⎪⎩⎭⎝⎭
∑ （59）
由于刚度矩阵[]0K 、[]K σ与时间t 无关，所以求解时可把（54）式写成以节点位移增量为基本未知量的增量形式，即
[][](){}{}0K K U P σ+∆=∆
（60）
显然，需用增量法求解，而把第i-1增量步所得[]v σ和[]d σ作为第i 增量步的初应力，按（7.236）式求得该增量步的节点位移增量{}U ∆，把i-1步所得单元形状作为计算单元刚度矩阵和单元体、面力节点载荷增量的体积和面积。

这样求解，就是把第i 增量步起始时刻（t 时刻）的构形作为参考构形，把第i 增量步终了时刻()t t +∆时刻的构形作为现时构形。

所以，看起来
似乎用的是Lagrange 描述方法，只不过是逐步修改的Lagrange 方法，即U.L 法。

大变形弹塑性动力显式有限元法
动力显式有限元算法的理论基础是连续介质振动理论中的动力学方程和中心差分格式算法。

首先，我们先简单地回顾一下最基本的单自由度动力学方程及其解法。

1 单自由度系统
对于只有一个自由度的线性弹簧系统，如图7-7.1所示，该系统的平衡方程为：
n n n n P kx cv ma =++ (7.237)
其中，n 代表时间增量步的第n 瞬时。

图 7－7.1 线性弹簧系统
P, x, v, a
m
(1)．无阻尼系统的自由振动
如果(7.237)式中的阻尼系数等于零，而且右端载荷项也为零，则(7.237)式可以写成无阻尼的自由振动方程，
0=+n n kx ma (7.238)
令m
k =
2
ω 方程(7.238)改写为
02=+n n x x
ω (7.239) 其解为
t A t A x n ωωsin cos 21+= (7.240)
式中，A 1和A 2为待定常数，由系统的初始条件确定。

设t=0时，
00,x x x x ==，则由(7.240)式得ω0
201,x A x A =
= 式(7.240)还可以变换
为相角形式的简谐函数，
)cos(ϕω-=t A x n (7.241)
式中，1
2
2
221,A A arctg
A A A =+=
ϕ 简谐函数的性质可用图(7-7.2)来表示。

图7-7.2 简谐振荡
(2)．有阻尼系统的自由振动
如果(7.237)式中仅右端载荷项为零，则(7.237)式可以写成有阻尼的自由振动方程，
0=++n n n kx cv ma (7.242)
引入临界阻尼系数C c ，
(7.243)
(1) 欠阻尼（C<C c ）（图7-7.3） 方程（7.242）的解为
)1cos(2ϕωςςω--=-t Ae x t n (7.244)
(2) 临界阻尼（C=C c ）（图7-7.4） 方程（7.242）的解为
t n e t x x
x x ςωω-++=))((000 (7.245) (3) 过阻尼（C>C c ）（图7-7.5） 方程（7.242）的解为
t s t s n e A e A x 2121+= (7.246)
其中，1
2)1(,1
2)1(2
02022
0201--+--=
--++=
ςωωςςςωωςςx x
A x x
A
ωςς)1(22,1-±-=S
图7-7.3 欠阻尼情况
图7-7.4 临界阻尼情况
]
图7-7.5 过阻尼情况
2、动力学平衡方程的中心差分求解格式
为了求解动力平衡方程（7.237），
n n n n P kx cv ma =++
将时间历程分成无数个离散的时间点，每两个时间点之间称为时间步长，这样，通过中心差分法就可以计算出每个时间点处的位移、速度、加速度、应力、应变…等物理量。

中心差分方法的具体推导过程如下（见图8.6）：
设t 时刻的状态为n ，t 及t 时刻之前的物理量已知，定义：
设：
1-∆-n t t ：n-1状态；
121-∆-n t t ：21
-n 状态； n t t ∆+21：21+n 状态；
n t t ∆+：n+1状态。

再设：1
-∆∆=
n n
t t β，并将（7.237）式 n n n n P kx cv ma =++
中的速度n v 和加速度n a 分别用差分格式改写为：
212
1
111-++++=
n n n v v v β
β
β
（7.247）
)()1(2
21211
-+--∆+=n n n n v v t a β （7.248）
n t t ∆+（n+1状态）时刻的总位移可以累加得到
n
n n n t v x x ∆+=++2
11 （7.249）
将（7.247）、（8.248）两式带入（7.237）式，并令m c ⋅=α（比例阻
尼）得
图7-7.6 离散的时间点
(n ) (21+n ) (n +1)
(21-n )
(n -1)
n
t t ∆+21n
t t ∆+1
21-∆-n t t 1
-∆-n t t
t
1
-∆n t n
t ∆
n n n n n n n kx P v v m v v t m -=++++-∆+-+-+-)11
1()()1(2212121211β
ββαβ
（7.250）
由上式解出：
)()2()1(2211
11212
1n n n n n n n n kx P m
t t v t t v -⋅∆+∆++∆+∆-=
-----+αββαβα （7.251） 对于多自由度系统，只要将（7.247）~（7.251）式中的v, a, x, m, k, P 分别用矢量或矩阵表示，即
{}{}{}212
1
11
1-+
+++=
n n n
v v v β
β
β
（7.252） {}{}{})()1(2
2
1211
-+--∆+=n n n n
v v t a β （7.253） {}{}{}n n n n t v x x ∆+=+
+2
11 （7.254）
[]{}{}[]{}{}{}[]{}n
n n n n n n x k P v v m v v t m -=++++-∆+-+-+-)11
1()()1(221212
1211βββαβ
（7.255）
如果[m]可以写成对角矩阵，则上式又可以写成分量的形式，从而得到速度
的显式表达格式，
)()2()1(221111212
1
i
n i n i
n n i n n n i
n kx P m
t t v t t v -⋅∆+∆++∆+∆-=-----+αββαβα （7.256） 其中，i 表示第i 个自由度。

为了保证系统计算的稳定性，中心差分法对允许使用的时间步长有一个限制（参见振动分析的数值方法），即，
cr t t ∆≤∆ （7.257）
式中，cr t ∆表示离散系统的临界时间步长。

对于单自由度系统而言，临界时间步长按下式计算
k m
t cr 22
0==∆ω （7.258）
对于多自由度系统，临界时间步长取决于系统的最高频率分量，即
m cr t ω2
=∆ （7.259）
其中，m ω表示系统的最高频率分量。

3. 动力显式算法
(1) 动力显式有限元方程
取出一个微小的正平行六面体，它在x 和y 方向的尺寸分别为dx,dy,它在z 方向的尺寸取单位长度。

在x 轴方向，有
dy
y
y y ∂+
σ
dx x x
∂∂σ
01)(11)(11)(=⨯⨯--+⨯-⨯∂∂++⨯-⨯∂∂+dxdy u c u
P dx dx dy y
dy dy dx x x x x yx yx yx x x
x ρτττσσσ等式两端同除以dxdy 得
0=--+∂∂+∂∂x x x yx
x u c u
P y
x ρτσ 同理，沿y 轴方向有，
0=--+∂∂+∂∂y y y xy y u c u
P x
y
ρτσ 采用动力显式算法，综合上面两个方程，得板料的运动学微分方程为，
(7.260) 式中，
ρ是材料的质量密度，c 是阻尼系数，i u
和i u 分别是材料内任一点的速度和加速度，p i 是作用在该点上的体积力，ij σ是该点处的Cauchy 应力。

根据散度定理
ij
i j
u dV
x i i ij ij
V
q u d
dV
以及边界条件，由(7.260)式，可以得到系统的虚功方程为
(7.261)
式中，i u
δ是虚速度，ij εδ 是对应于Cauchy 应力ij σ的虚应变速率。

把物体离散成m 个单元，对于任一单元，有α个结点，取其形函数为
αN ，单元内任意点的位移分量i u 、速度分量i u
和加速度分量i u 分别为
⎪⎩⎪⎨⎧===ααααααi i
i i
i i u N u u N u u N u (7.262)
由几何方程可得
ααεi j ij u B = (7.263)
式中，α
i u 、α
i u
和α
i u 分别为结点α的位移分量、速度分量和加速度分量，α
j B 为第五章给出的应变张量（具体求解时，采用线性化处理，可以按照小变形的应变矩阵考虑，但要进行构形更新）。

将(7.262)式和(7.263)式代入(7.261)式，可得：
⎰
⎰⎰⎰⎰
-Γ+=
+Γe e
e
e
e
V i j ij i i V i i V V i i i i dV u B d u N q dV u N p dV u N u cN dV u N u N
βββββββ
βααββααδσδδδδρ (7.264)
式中βδi u
是结点β的虚速度，上式还可以写成 ⎰⎰⎰⎰⎰
-Γ+=
+Γe e
e
e
e
V j ij i V i V V i i dV
B d N q dV N p u dV N cN u dV N N
β
ββα
βααβασρ (7.265)
写成矩阵形式为
(7.266)
将单元方程集合，即得整体有限元方程
∑⎰∑⎰∑⎰∑⎰∑⎰-Γ+=+Γe e
e
e
e V T T V T V T V T
dV
d dV dV c dV
σB q N p N u N N u
N N
)()(ρ (7.267)
上式可简写为
F P u C u
M -=+ (7.268) 式中，M 是一致质量矩阵，∑⎰=e
V T
dV N N M ρ；
C 是阻尼矩阵，
∑⎰=e
V T
dV c N N C ； P 是结点外力向量，∑⎰∑⎰ΓΓ+=e
e
d dV T
V T q N p N P ；
F 是结点内力向量，
∑⎰=e
V T
dV σB F 。

通常，动力显式积分算法采用集中质量矩阵，既M 是一对角矩阵，并取C =aM 。

则(7.268)式表示的联立方程组变成了（结点数⨯结点自由度数）个相互独立的方程
i i
i i
i
i
u
c u
P
F
这里对i 不求和，i ＝1～总自由度数。

(2) 临界时间步长的计算
由于中心差分格式的算法是条件稳定的，为了保证系统计算的稳定性，对时间增量步长t ∆的大小必须加以限制。

稳定性条件通常由系统的最高频率max ω决定，满足稳定性条件的时间增量步长为
)1(2
2max
ξξω-+≤
∆t (7.270)
其中ξ是具有最高频率的模型的临界阻尼比。

与工程直觉相反，阻尼的引入实际上降低了系统的临界稳定性条件。

系统的最高频率由网格中最大的单元膨胀模式决定。

满足稳定性条件的时间增量步长可以由膨胀波沿网格中任意单元的最小穿越时间近似得到
(7.271)
其中8.0~6.0=γ，c 为膨胀波在材料中的传播速度
()E
（7.272） ，e n L 为第n 状态（t n 时刻）单元e 的名义长度。

稳定性条件可以保证在一个时间增量步内，扰动只传播网格中的一个单元。

如果系统只包括一种材料，则满足稳定性条件的时间增量步长与网格中最小的单元尺寸成正比；如果系统划分的单元网格尺寸比较均匀，但包括多种不同材料，则具有最高膨胀波速的材料中网格尺寸最小的单元决定系统的稳定时间步长。

对于一个简单的桁架单元，在团聚质量矩阵的情况下，稳定性准则给出一个临界时间步长：c
l
t ≤
∆，式中，c 为材料声速，l 为单元长度，t ∆表示膨胀波穿越长度为l 的单元所需要的时间。

这就是所谓的Courant-Friedrichs-Lewy(CFL)稳定性条件。

对于三角形单元和四边形板单元来说，临界时间步长的选取依赖于单
元名义长度的确定，一般按照图8.7的原则来确定单元的名义长度lcrit ，对于高阶单元来说，临界时间步长远比低阶单元的临界时间步长来的小， 这
3. 接触问题的动力平衡方程
对于考虑阻尼的接触问题的动力平衡方程又可表达为
(7.273)
其中，M 是质量矩阵；C
是阻尼矩阵；f int 和f ext 分别是结点内力和外力； λ是摩擦接触反力；u,u 和u 分别是结点的位移、速度和加速度向量。

采用中心差分法求解，参照（7.256）式得:
d 1
d 2
),m in(21l l l crit =
),m ax (/21d d A l crit =
max /2s A l crit =
图8.7 单元名义长度
))()(()2()1(22int 1111212
1
i
n i n i ext n i n n i n n n i n f f m
t t u t t u
λαββαβα+-⋅∆+∆++∆+∆-=-----+ (7.274)
4材料本构列式和应力计算
基于动力显式算法的大变形弹塑性问题可以采用速率型的本构方程。

采用具有客观性的Cauchy 应力ij σ的Jaumann 导数
ki kj kj ik ij ij ωσωσσ
σ--=∇
（7.275） 其中ij σ
是Cauchy 应力的物质导数，ij ω是旋转率张量。

大变形弹塑性速率型（或增量型）本构方程与小变形弹塑性增量本构方程在形式上类似，只要把其中的应力和应变改用Cauchy 应力的Jaumann 导数和变形率，应力偏量和等效应力分别改用Cauchy 应力偏量和等效的Cauchy 应力即可。

本构方程的张量分量形式为
kl ep ijkl ij d D =∇
σ （7.276）
式中ep
ijkl D 是弹塑性本构矩阵，kl d 是变形率张量。

t t ∆+时刻(n+1状态)构形下的Cauchy 应力1
+n ij σ可以表达为（参考图
8.6）
t
n ij
n ij
n ij
∆+=++21
1
σ
σσ
（7.277）
实际计算时引入了下面的差分格式，令
t n ki
n jk
n kj
n ik
n ij n ij
∆++=+++)(ˆ21
21
1ω
σω
σσσ
（7.278）
式中⎪⎪
⎭
⎫ ⎝
⎛∂∂-
∂∂=+++j n i i n j
n ij
x v x v 2121
21
21ω
为旋转率张量。

则得到t t ∆+时刻构形下的Cauchy 应力
t d
D n kl
n ep ijkl
n ij
n ij
∆+=+++21
)(ˆ11σ
σ
（7.279）。