?(??NVeVeT????(?cNTNdV)u?NdV)uVeTT???NpdV???Nqd????B?dVT?eVe? (7.267)
上式可简写为
???Cu??P?F (7.268) MuT式中，M是一致质量矩阵，M????NNdV；
VeTC?cN? C是阻尼矩阵，?NdV；
VeTTP?NpdV?N P是结点外力向量，????qd?；
Ve?eTF?B???dV。 F是结点内力向量，
Ve 通常，动力显式积分算法采用集中质量矩阵，既M是一对角矩阵，并取C=aM。 则(7.268)式表示的联立方程组变成了（结点数?结点自由度数）个相互独立的方程
iiiii&&&mu+cu=P-Fi (7.269) 这里对i不求和，i＝1～总自由度数。 (2) 临界时间步长的计算
由于中心差分格式的算法是条件稳定的，为了保证系统计算的稳定性，对时间增量步长?t的大小必须加以限制。稳定性条件通常由系统的最高频率?max决定，满足稳定性条件的时间增量步长为
?t?2?max(1??2??) (7.270)
其中?是具有最高频率的模型的临界阻尼比。与工程直觉相反，阻尼的引入实际上降低了系统的临界稳定性条件。系统的最高频率由网格中最大的单元膨胀模式决定。
满足稳定性条件的时间增量步长可以由膨胀波沿网格中任意单元的最小穿越时间近似得到
Le?tn??n (7.271)
c其中??0.6~0.8，c为膨胀波在材料中的传播速度
r1/2c=() （7.272）
E，Len为第n状态（tn时刻）单元e的名义长度。
稳定性条件可以保证在一个时间增量步内，扰动只传播网格中的一个单元。如果系统只包括一种材料，则满足稳定性条件的时间增量步长与网格中最小的单元尺寸成正比；如果系统划分的单元网格尺寸比较均匀，但包括多种不同材料，则具有最高膨胀波速的材料中网格尺寸最小的单元决定系统的稳定时间步长。
对于一个简单的桁架单元，在团聚质量矩阵的情况下，稳定性准则给出一个临界时间步长：?t?l，式中，c为材料声速，l为单元长度，?t表c示膨胀波穿越长度为l的单元所需要的时间。这就是所谓的Courant-Friedrichs-Lewy(CFL)稳定性条件。
对于三角形单元和四边形板单元来说，临界时间步长的选取依赖于单
元名义长度的确定，一般按照图8.7的原则来确定单元的名义长度lcrit，对于高阶单元来说，临界时间步长远比低阶单元的临界时间步长来的小， 这一事实使得高阶单元对于显式积分算法相当不合适。 l2 l1 d1
d2 s1 s1 s1 lcrit?min(l1,l2)
lcrit?A/max(d1,d2)
图8.7 单元名义长度
lcrit?2A/smax
虽然上面给出的稳定性准则严格来说是对线性系统而言的，但对于非线性问题也给出了有用的稳定性估计。对线性问题时间步长的80%~90%缩小，对于大多数非线性问题保持其系统稳定性是足够的。然而，十分重要的是：在整个计算过程中，要不断的检查能量的平衡问题。任何总能量的增加或损失（5%或更多）都将导致失稳。Belytschko指出：常增量时间步不能保持解的稳定性，即使系统的最高频率3. 接触问题的动力平衡方程
对于考虑阻尼的接触问题的动力平衡方程又可表达为
?max不断减小。 &&&fMu+Cu+=f+l (7.273)
其中，M是质量矩阵；C是阻尼矩阵；fint和fext分别是结点内力和外力；??&&&是摩擦接触反力；u,u和u分别是结点的位移、速度和加速度向量。采用
中心差分法求解，参照（7.256）式得:
?uin?122???tn?1i(1??)?tn?1extiintii?n?1??u?((f)?(f)??nnn)i2????tn?12(2????tn?1)m (7.274)
4材料本构列式和应力计算
基于动力显式算法的大变形弹塑性问题可以采用速率型的本构方程。采用具有客观性的Cauchy应力?ij的Jaumann导数
?ij??ik?kj??kj?ki （7.275） ?ij????ij是Cauchy应力的物质导数，?ij是旋转率张量。 其中?大变形弹塑性速率型（或增量型）本构方程与小变形弹塑性增量本构方程在形式上类似，只要把其中的应力和应变改用Cauchy应力的Jaumann导数和变形率，应力偏量和等效应力分别改用Cauchy应力偏量和等效的Cauchy应力即可。本构方程的张量分量形式为
ep?ij?Dijkldkl （7.276）
ep式中Dijkl是弹塑性本构矩阵，dkl是变形率张量。
?
n?1t??t时刻(n+1状态)构形下的Cauchy应力?ij可以表达为（参考图
8.6）
?n?1ij?????nijn?12ij?t （7.277）
实际计算时引入了下面的差分格式，令
??n?1ij???(??nijnikn?12kj???njkn?12ki)?t （7.278）
式中?n?12ijn?2n?1??v?vi21?j???2?xi?xj?1???为旋转率张量。 ?则得到t??t时刻构形下的Cauchy应力
?n?1ij???n?1ij?(D)depnijkln?12kl?t （7.279）