边界面力的一致力为(s为面力矢量)：
nfTi??NinTsnd? (4.216)
n?式中, ?为受面力作用的边界与单元边界重合的部分。
利用Gauss-Legendre乘积方法，依据单元形函数，可以对上述各式进行数值计算。 设材料的粘弹性本构方程???e??ve??vp中各部分如下： 1) 弹性应变
?e?其中D为弹性矩阵，设为常数矩阵。
2) 粘弹性应变增量
粘弹性应变率可以表示为
n?ve??(1? (4.217) D1nE2 (4.218) A?n??ve)E2?1即
n?ve??E1?E2nE2 (4.219) A?n??veE2?1?1在tn?1时刻为
n?1?ve??E1?E2n?1E2 (4.220) A?n?1??veE2?1?1n?1nn???????? (4.221) ?n?1nn????????且有
采用如下插分格式
n?1nn?1?ve?ve??ve??tn[(1??)????] (4.222)
其中，?为插分系数，如0，0.5，1等。则粘弹性应变增量与应力增量之间的关系为
n??ve?1???tnE2n?ve???E2A?tn??E2??1???tn??1??1????n (4.223)
?1?tn 3)粘塑性应变增量
当满足一定的屈服准则时，材料产生粘塑性变形。采用类似的插分方程
n?1nn?1?vp?vp??vp??tn[(1??)????] (4.224)
?vp用Taylor级数来近似 其中?n?1n?1n?vp?vp????Hn??n (4.225)
n其中H由屈服准则和粘塑性流动准则确定
?aTd???TnH?()??{????aa} (4.226)
????dfnn?vp?? 采用Ducker-Prager准则：
??k? (4.227) ?J1?J2?为应力偏量的第二不变量，式中?,k?由材料粘聚力和内摩擦角定；J1为应力第一不变量；J2当取关联的流动法则时有
a??n?J2?Q?f?J???(1)n?() (4.228) ????????这样，由式(4.224)和式(4.225)得到
nn?vp??vp???tn??Hn??ntn (4.229)
4)总应变增量
由应变分解可知
nn (4.230) ??n???en???ve???vp其中总应变增量可以由结点位移得到
??n?Bn?dn (4.231)
n式中, B为位移应变矩阵；?d为单元结点位移增量。
n由式(4.223)、式(4.229)和式 (4.230)可得应力增量为：
?1n?nnn????n?Dvevp?????t???t?venvpn? (4.232)
1??E2?tn/?1??n是粘弹粘塑性模量 Dvevp???E2A?tnnDvevp??D?1???Hntn? (4.233)
(1??E2?tn/?1)?1??5) 动态有限元刚度矩阵
在任一tn时刻，可以证明，该时刻的切线刚度矩阵为
nnKT??BnTDvevpBndV (4.234)
V?1 3. 显式时间积分法
任一tn时刻的动平衡方程可以写成如下矩阵形式：
??n?Cd?n?Pn?fn (4.235) Md式中M为总质量矩阵，C为总阻尼矩阵，P为内力的总矢量，fn为作用的体力、面力的一致结点力矢量。
为简便起见，用中心差分公式对上述动平衡方程进行离散化。加速度、速度分别为：
n??n?d1nnn?1{d?2d?d} (4.236) 2(?t)?n?1{dn?1?dn?1} (4.237) d2?t通过差分处理，在tn?1时刻的位移可用tn?1,tn时刻的位移显式地表过出来。
根据上述思路，可以对这类非线性材料进行动力有限元计算，得到相应的动应力、动应变响应等结果。
9.4 应用举例
4.6.4 粘弹粘塑性动力有限元分析举例
编制粘弹粘塑性动力有限元程序，对受到振动压实物料的位移及应力应变响应特性进行了如下研究。
假设物料是均匀的，形成粘弹粘塑性材料结构，如图4-8所示为该结构某断面的平面四边形单元网格划分。在上边界的某一点受到周期性集中载荷的作用，载荷具有不对称性。两侧边界无约束。与基础层结合处的底边结点均受约束。假设材料满足等向强化的Drucker-Prager屈服条件。材料参数如表4-3。
表4-3 材料参数取值 Poision比 弹性 模量 弹性 模量 粘滞系数 粘滞系数 粘聚力 内摩擦角 强化系数 密度 流动性系数 ? 0.3 E1 6.8804 E2 2.164 ?1 ?2 (N/mm2) (N/mm2) (sec-1) (sec-1) (N/mm2) 7.5e-5 1.9e-7 7.07 c ? ? ? (kg/mm3) 3.0e-6 ? 1.0e-2 62.73 1.0e-3 YZX
图4-8 受集中载荷的物料层的网格划分
设载荷P(t)为
2?vt;当P(t)?0?b0sin?P(t)??当P(t)?0 (4-174) ?0;?其中，v?60Hz，b0??1.0?105N。取时间步长为dt=3.333X10-6sec。在不同的时刻，材料断面的应力场分布不同。某结点垂直方向(y方向)上的应力时间历程如图4-9（a）所示。由于作用力的不对称性，应力应变关系呈现近似的不对称滞回特性，如图4-9(b)所示。
(a) (b)
图4-9 材料断面某结点垂直方向的应力历程及其应力应变关系
习题