超弹性平面应变问题的显式非线性有限元程序2D+Explicit+nonlinear+FEM+for+Rubber.docVIP

MACROBUTTON MTEditEquationSection2 SEQ MTEqn \r \h SEQ MTSec \r 1 \h SEQ MTChap \r 1 \h 非线性有限元 非线性有限元 课程作业 班级：力 学 学号：姓名：姜 琛 2013年1月13日 超弹性平面应变问题的显式非线性 有限元程序 力学 姜琛 程序概述 本程序是课程非线性有限元的期末作业，程序主要参考Ted Belystchko的著作《Nonlinear Finite Elements fro Continua and Structures》一书中关于几何非线性的完全拉格朗日格式（Total Lagrangian formulation, TL）的相关资料，以及改编Jack Chassa的三角形单元的显式TL格式大变形平面应力有限元Matlab程序得到。并参考了M.A Crisfield的著作《Non-linear Finite Element Analysis of Solids and Structures Vol.2》和Javier Bonet, Richard D. Wood的《Nolinear Continuum Mechanics For Finite Element Analysis》中，关于超弹性本构关系以及应力体积项和偏斜项的分解公式。与Jack Chassa的程序的主要区别在于，本程序在其基础上添加了超弹性材料的Mooney-Rivlin模型的相关代码，以及采用Q4单元的局部减缩积分来减少单元体积自锁。然后使用所编制程序，模拟了橡胶块均匀受压案例，并与相同案例的Abaqus结果对比，验证程序的正确性。 二．理论基础 根据非线性有限元课的要求，求解一个二维超弹性材料压缩问题，并且使用显式非线性有限元求解。该类问题是大变形小应变的问题，可采用如下的公式组合： 分析类型 描述 使用格式 应力应变度量 大变形大转动小应变 材料纤维的位移和转动很大，但是纤维的伸长和纤维间的角度很小； 应力应变关系可为线弹性或非线性的； Total Lagrangian （T.L.）； Updated Lagrangian （U.L.）； 第二Piola-Kirchoff应力张量，Green应变张量； Cauchy应力张量，Almansi应变张量； 一、超弹性材料本构方程(Hyper-elastic Constitutive Equations) 由于T.L.格式和U.L.格式本质上是相同的，只是应力和应变的参考构形分别为初始构形和当前构形，在此程序编制中采用了T.L.格式求解非线性问题。另外由于橡胶是超弹性材料，所以需要用应变能密度函数W求第二Piola-Kirchoff应力张量，其一般公式为（指标表示法）： 其中分别是变形三不变量；J是变形梯度张量的行列式值，在不可压缩条件下，J=1；是右Cauchy-Green变形张量，且。并且上式中各不变量的导数项可如下计算： 则式 GOTOBUTTON ZEqnNum977433 REF ZEqnNum977433 \* Charformat \! (1.1)也可以转换到如下格式： 在此使用超弹性材料所常用的Mooney-Rivlin模型，不可压条件下，其应变能密度公式如下： 此处的是描述材料剪切特性的参数，则PK2应力张量的表达式可以化为： 另由于本构关系采用应变能函数表示，而显式分析中需要条件稳定时间，其中弹性模量。 然而对于不可压缩材料，即泊松比为0.5时，将带来很严重的计算困难（除平面应力问题以外，和压力作为自变量的混合弱形式，而本文要求计算平面应变问题）。需要对所选用的超弹性应变能函数模型进行修正,以便于分解出体积项和偏斜项，并引入适当的可压缩性。其中采用的修正的Mooney-Rivlin模型如下[M.A. Crisfield]： 其中 且在初始条件下,当J=I3=1时，应变能为0。再导出PK2应力的形式为： 其中 并参照ABAQUS中Explicit的做法，加入一个微小可压缩性，泊松比，其中泊松比的方程为 则有，且，。因此可求得体积模量。 二、体积自锁 (Volumetric Locking ) 体积自锁问题，来源于单元无法表征等体积变形导致。而对于不可压缩材料，或近乎不可压缩材料都会伴随着较大的等体积变形，所以需要进行一些特殊处理。一般, 有应力的局部减缩积分方法，以及采用混合弱形式的控制方程，即杂交单元(Hybrid Elements)，增加一个压力的自由度，通过动量方程求出静水压。这里采用较简单的局部减缩积分(