（工程力学专业论文）基于Mooney—Rivlin本构模型Ⅰ型裂纹的有限元分析.pdf

北京交通大学硕士学位论文 y7 4 1 0 8 5 摘要 本文采用非线性有限元法对橡胶类材料i 型裂纹尖端区域进行了计算， 给出了不同大小均布拉应力作用时该区域的变形形貌，以验证扩张区与收 缩区存在的理论预测。 m o o n e y r i v l i n 本构模型是橡胶材料中具有代表性的本构模型，有限元 法是分析大变形问题的有效方法。本文运用大型通用有限元软件a n s y s ， 采用m o o n e y - r i v l i n 应变能函数，分析了平面应变i 型裂纹尖端的变形和应 力分布。数值结果表明：在裂尖附近不仅存在应力奇异性，而且变形存在 扩张区与收缩区；在裂尖附近较大的主应力是拉应力，其方向和裂纹方向 接近垂直。拉应力比另一个主应力高大约两个数量级，所以其对裂尖附近 的变形起主导作用；若以变形前状态为参考构形观察应力分南，则扩张区 与收缩区内的最大拉应力具有几乎相同的奇异性。 关键词：m o o n c y r i v l i n 应变能，裂尖场，有限元法，大变形 北京交通大学硕七学位论文 a b s t r a c t t h e r e g i o no fm o d e 1c r a c k t i pi nr u b b e r m a t e r i a li sc a l c u l a t e db yn o n l i n e a r f i n i t ee l e m e n tm e t h o d ( f e m ) a n dt h ed e f o r m a t i o nc o n f i g u r a t i o no f t h i sr e g i o ni s g i v e no u t u n d e rd i f f e r e n tu n i f o r mt e n s i l es t r e s s pt op r o v e t h et h e o r e t i c a l p r e d i c t i o na b o u t t h ee x i s t e n c eo f e x p a n d i n g a n ds h r i n k i n gs e c t o r s m o o n e y - r i v l i nc o n s t i t u t i v e f u n c t i o ni st h er e p r e s e n t a t i v eo n ei nr u b b e r m a t e r i a l sa n df e mi sae f f i c i e n tm e t h o dt o a n a l y z et h el a r g ed e f o r m a t i o n p r o b l e m t h em o d e 1c r a c kt i pf i e l du n d e r p l a n es t r a i nc o n d i t i o ni nh y p e r e l a s t i c m o o n e y r i v l i n m o d e lm a t e r i a li s i n v e s t i g a t e d a n u m e r i c a l a n a l y s i s i s p e r f o r m e dt oa n a l y z et h ed e f o r m a t i o na n dt h ed i s t r i b u t i o no fs t r e s sf i e l dw i t h g e n e r a lf i n i t ee l e m e n ts o f t w a r ea n s y s t h er e s u l t si n d i c a t et h a t ：( 1 ) n o to n l y t h es t r e s ss i n g u l a r i t yb u ta l s ot h ee x p a n d i n ga n d s h r i n k i n gs e c t o r se x i s tn e a rt h e c r a c kt i p ；( 2 ) i nt h et w o p r i n c i p a ls t r e s s e s ，t h el a r g e rs t r e s si sh i g h e ra b o u tt w o o r d e r so fm a g n i t u d et h a nt h eo t h e ro n e t h el a r g e rs t r e s si sat e n s i l es t r e s sw i t hi t s d i r e c t i o na p p r o x i m a t e l yv e r t i c a lt ot h ec r a c kt i pa n di t p l a y st h el e a d i n gr o l ei n t h ed e f o r m a t i o na tt h ec r a c k t i p ；( 3 ) i ft h e s t r e s s s i n g u l a r i t y i so b s e r v e di n u n d e f o r m e dc o o r d i n a t es y s t e m ，t h el a r g e s tt e n s i l es t r e s sh a sa l m o s tt h es a m e s i n g u l a r i t yi nt h ee x p a n d i n ga n ds h r i n k i n gs e c t o r s k e y w o r d s ：m o o n e y - r i v l i n s t r a i ne n e r g yf u n c t i o n ，c r a c kt i pf i e l d ，f i n i t e e l e m e n t m e t h o d ，l a r g ed e f o r m a t i o n j i 北京交通大学颂士学位论文 1 1 引言 第一章绪论 随着新材料的不断出现和计算精度要求的不断提高，线弹性理论在某些 领域已经不能满足要求。对于某类材料，受外力作用其应力应变不再保持 线性关系，但当载荷完全消除时，材料变形则完全恢复，橡皮、橡胶、高 分子聚合物等都属于这类材料。由于这类材料弹性非常好，且能产生很大 的变形，小变形理论已不能再满足要求。因此，人们建立了相应的非线性 理论。 通常，把非线性问题分为两大类，即几何非线性和材料非线性。所谓材 料非线性，是指非线性效应仅由应力应变关系的非线性引起，位移分量仍 假设为无限小量：几何非线性，其中包括两种情况：1 、大位移、大转角， 但是应变很小；2 、应变不再是无限小量。此外，还有一类非线性问题，即 边界条件非线性，即由于边界条件的性质随物体的运动发生变化而引起的 非线性响应。 对于橡胶材料的大变形问题，将涉及前两类非线性问题，而且在几何 非线性中，应变不是无限小量，即有限变形。 有限变形理论是弹性力学的一个分支，其在橡胶材料和高分子聚合材 料中有着广泛的应用前景，因而它的理论研究一直受到力学界和工程界的 重视。 早在1 6 9 1 年，b e r n o u l l i 通过实验发现弦的伸长和张力不满足线性的 h o o k e 定律。在此之后，通过l a n g r a n g e 和c a u e h y 等人的努力，提出了连 续介质力学的一些基本假设，有限变形理论的一些基本概念如应力、应变 等才建立起来。1 8 9 4 年f i n g e r 建立了超弹性体的有限变形理论，非线性理 论的基础得以建立。由于橡胶材料的双重非线性即几何非线性和物理非线 第一章绪论 性，使得有限变形弹性理论方程繁琐而复杂，特别是强烈的非线性，使得 当时在数学上对其进行一般性讨论没有多大的希望，因而这方面的研究难 以取得进展。但从1 9 4 0 年起，非线性弹性理论的研究得到了迅速的发展， 这一时期，出现了很多有价值的论文，并成功地解决了许多工程实际问题， 以文献l l 】代表了当时研究的最高成果。e r i n g e n l 2 1 、郭仲衡f 3 等学者在这方面 做出了较大贡献。在此之后，非线性研究进入了蓬勃发展的时期，一方面 非线性理论的逐步成熟，另一方面随着计算机的发展，大变形计算的研究 也得到了发展l “l ，有限变形弹性理论日臻成熟，工程中考虑大变形也逐渐 不可忽略。但非线性弹性力学固有的复杂性，使其无论在理论方面还是数 值求解上仍存在很多困难，有待于进一步研究。 l2 有限变形弹性体的本构关系 本构关系是反映物质性质的数学模型。一物体在外部因素作用下的响应， 必须知道描述构成物体的物质属性所特有的本构方程，因而寻求一种合理 而适用的本构方程一直受到人们的重视。 对于橡胶材料，其有如下丰要变形特点： 1 、大部分橡胶材料在受力状态下体积没有明显的变化，具有不可压缩 性； 2 、应力和应变关系是非线性的但加载和卸载是同条曲线； 3 、橡胶材料可承受大变形，因此不能用小变形理论描述其变形特征。 由于以上变形特点，使得其本构关系极为复杂。 从理论上讲，建立橡胶类材料即有限变形弹性体的本构关系，通常有 两种途径： 1 、c a u c h y 方法：从于弹性体的特征出发，即“一定的应力状态对应于 一定的应变状态”，通过分柝实验结果，假设一定的应力应变函数关系，再 通过实验确定其中的系数。直接由这种应力应变函数关系来描述的物体称 为柯西意义的弹性体。 2 、g r e e n 方法：即从势能函数出发来得到弹性体的本构方程，具有弹 性势的弹性体称为超弹性体或格林意义下的弹性体。 北京交通大学硕士学位论文 柯西弹性体是一个比超弹性体更为广泛的概念。但用c a u c h y 方法直接 建立有限变形材料的本构关系在数学和物理上都相当困难，人们倾向采用 g r e e n 方法来建立本构关系，且采用g r e e n 方法的关键，在于选择一个能正 确反映材料变形的应变能表达式。 对于各向同性材料，假设l l 、，2 、，3 为g r e e n 应变张量的三个不变量， 则其应变能彬可以表示为： w = w ( ，1 ，1 2 ，1 3 ) ( 1 2 - 1 ) 1 1 ；d ：e ，1 2 一d 2 ：e ，l = d 3 ：e ( 1 2 2 ) 11 l ：一专q ；- i z ) ，l ：= 之q ：3 i l f 2 + 2 i o ( 1 2 3 ) 其中d 为g r e e n 变形张量，e 为单位张量。 选取应变能函数矽不同，则描述物体的本构方程也不同，如何选择应 变能函数，反映材料变形的特征，是有限变形弹性理论的重要课题之一。 m o o n e y 【7 l 根据大量的实验结果，提出了橡胶材料的应变能形式为： w = c l ( j 1 3 ) + c 2 ( 一3 ) ( 1 2 4 ) 式中，c l 、c 2 为材料常数，这个函数可看作是在= 3 和f ；3 附近 的二重幂级数展开项的前两项。m o o n e y 形式的本构方程可以较好地拟合橡 胶材料中等应变范围的实验【聃l 。当c 1 ：0 是上式最简单的形式，r i v l i n 【l o 】称 这样的材料为n e o h o o k e a n 模型。 h a r t s m i t h 提出将i 矿取作两个函数暇“) 、( e ) 之后，即 矽= 瞩“) + ( e ) ( 1 2 5 ) 对于不可压缩材料，许多学者还提出了不同的应变能形式，较有代表 性的有： h a r t s m i t h “l ；s h i e l d 1 2 一1 3 l ；s t o t a k e r s l l 4 1 ： 等，很多都是( 1 2 4 ) 式在无变形状态处的级数展开，所不同的的是保留的级 第一章绪论 数项不同。 随着橡胶工业的发展，一些新型材料不断出现。泡沫聚胺脂橡胶是一 种典型的可压缩橡材料，因而一系列反映可压缩橡胶材料力学行为的应变 能函数也就相应而牛。由于要描述体积变形的影响，应变能表达式中必须 包含表示体积变化的不变量k 。b l a t z 和k o i 、b l a t z l l 翻、m u r n a g h a n l l 7 j 等学 者分别提出了可压缩材料的应变能函数形式。 o g d e n l l8 】通过引入改写的主伸长： 1 。= i ，3 ，j = 屯 ( 1 2 6 ) 将应变能函数记为： w = 矿( ，墨，j ) ( 1 2 7 ) k n o w l e s 和s t e r n b e r g l l 9 l 给出了应变能形式： t ( a 6 + 彤+ c i i j ) ” ( 1 2 8 ) 式中，a ，b ，c ，n 为材料常数，j 为变形前与变形后的体积比。 g a oy c 1 2 0 l 从分离形状变化和体积变化的角度出发，引入了表示形状 改变和体积改变的不变量： r = 1 1 k 1 ” k t 给出了一种应变能的具体形式，表示为： w a ( r 5 3 4 ) + 6 ( 鬣一1 ) 4 k 一9 式中，a , b ，坍，n ，q 为材料常数。 推导的本构方程为： tz 2 n a l ? 一1 k 亍札扣 tz 一1 3 2 ( d 一冬e ) + 2 b ( k 1 ) “k 1 一j 惭一q ) k + q 1 e ( 1 2 9 ) ( 1 2 1 0 ) ( 1 2 1 1 ) ( 1 。2 1 2 ) 北京交通大学硕士学位论文 式中，d 为c a u c h y 应变张量，第一项为偏应力项，第二项为静水应力 项，由此也说明了方程( 1 2 1 1 ) 中第一项表示形状改变，第二项表示体积 改变。 对于不可压缩材料，则有k ：1 ，方程( 1 2 1 1 ) 简化为： w = 日( ，? 一3 ”) ( 1 2 1 3 ) 应力应变关系表示为： t = 2 n a ( 1 ；。d + 口e ) ( 1 2 1 4 ) 式中，o 为静水应力，必须由平衡方程确定。 g a o y c 1 2 1 1 从分离材料的抗拉和抗压能力出发，引入了另一种描述橡 胶材料的应变能函数： 肌n 卜钏 2 ， 式中，n 和疗为材料常数。 推导出来的术构关系为： t m 。卜一钞叫 2 舶， d 。1 为c a u c h y 应变张量d 的逆变张量。 用这一模型，既可以描述材料的受拉变形，也可以描述材料的受压变 形。 综上所述，描述超弹性橡胶类材料的应变能函数形式多种多样，不同 学者得到了不同的应变能函数形式。 1 3 橡胶类材料i 型裂尖场的研究进展 当非线性弹性体本构关系确定以后，就可以利用它进行一些具体问题 的渐近分析。使用最多的还是在裂尖场方面，因为断裂现象的发生，大多 也是由裂纹扩展引起的，研究此处应力应变，对裂纹扩展开裂有着重要意 第一章绪论 义。因此，考虑断裂问题时，裂尖附近的应力奇异性便成为一个重要问题。 橡胶类材料由于具有良好的性能在工程上运用非常广泛，由以上有限变 形弹性体的本构关系的介绍可知，其本构关系十分复杂，并且受力时变形 伴随着大位移、大应变，特别是裂尖附近的应力应变场不易准确描述，必 须用非线性弹性理论分析其附近应力分布。 w o n 2 和s h i e l d | 2 2 1 最早研究了不可压缩n e o - h o o k e a n 材料的裂纹问题， 在无限远处受双向拉伸，采用渐近方法，得到了裂纹问题的平面应力渐近 解。 k n o w l e s 和s t e r n b e r g l 2 3 壤先采用完全非线性弹性理论，研究了均匀各 向同性超弹性体的平面应变i 型裂纹问题，给出了裂纹附近弹性场的大变形 渐近分析；其后引进高阶渐近分析1 2 4 1 ，对文献【2 3 】的结果进行了修正。但 k n o w l e s 和s t e r n b e r g l 2 3 刨对裂纹附近弹性场进行的大变形渐近分析及其高 阶分析并没有把这个问题完全弄清楚。 g a oy c 2 0 , 2 1 搬据裂尖区域变形剧烈，无法用统一映射函数描述其变形 的特点，提出了将裂纹尖端区域划分为扩张区与收缩区，采用不同的映射 函数分别描述其变形，利用式( 1 2 1 1 ) 、( 1 2 1 3 ) 两种不同的本构关系， 应用于裂尖场的研究都取得了理论解，最终得到了裂尖应力应变场的理论 分析结果，较简明的解决了l 型裂纹问题。l i u l 2 5 】利用g a o 的分区方法和 k n o w l e s 和s t e i n b e r g 提出的式( 1 2 8 ) 本构函数也得到了相同的结论。这 种分区方法应用于其他问题的研究也取得了系列的成果1 2 2 1 。 1 4 本文的研究背景 至今，各位学者提出了各种有关橡胶类材料的应变能函数及本构方程， 对于橡胶材料l 型裂纹尖端场，扩张区与收缩区的存在与否，需要理论、实 验和数值三方面的研究和验证。裂纹尖端区域有关扩张区与收缩区的数值 方面的研究，据作者掌握的资料，还未见相关报道。因此，本文利用有限 元法进行这方面的研究工作，希望能得出一些有意义的结果。 随着计算机的发展，有限元法在断裂力学中应用较多，但目前对裂尖 场的研究主要集中在裂纹参数计算方面，包括能量释放率和，积分的计算 北京交通大学硕士学位论文 1 3 3 - 3 5 。据作者掌握的资料，用有限元法分析基于m o o n e y - r i v l i n 本构关系的 i 型裂纹尖端的变形形貌，目前还没有这方面的研究报告。由于描述橡胶类 材料的应变能函数形式较多，采用不同的应变能函数，可能会得出不同的 结果，g a o l 2 0 , 2 1 1 提出的将裂尖划分为扩张区和收缩区进行分析的方法对于 m o o n e y - r i v l i n 本构模型是否仍然适用，是一个值得研究的问题。 本课题得到了国家自然科学基金( 9 0 2 0 5 0 0 7 ) 和北京交通大学校基金 ( t j 2 0 0 2 j 0 1 8 0 ) 的资助。 1 5 本文的工作 本文的目的就在于采用非线性有限元法对橡胶类材料l 型裂纹尖端区 域进行计算，以图谱方式直观展示裂纹尖端位移场和应力场。根据裂纹尖 端区域附近应力分布图谱对该区域变形形状的产生予以解释，以验证扩张 区与收缩区的存在的理论预测。 本文运用大型通用有限元软件a n s y s ，采用m o o n e y r i v l i n 应变能函 数，计算了平面应变i 型裂纹尖端的变形和应力分布，并对其裂尖变形场的 特点进行了分析。具体工作如下： 1 、裂纹尖端的有限元网格细分，原始数据的输入。 2 、从最小载荷到最大载荷的加载过程中，给出了五个中间载荷对应的有限 元计算结果。其中，最大载荷是最小载荷的2 0 倍。 3 、整理了相关数据，绘出了不同载荷作用时裂失区域的变形图及应力分布 图。通过选择不同大小范围的裂尖区域，对变形图谱的变化趋势进行对比 分析，从应力分布图谱对进行了变形场的产生进行了解释，并和相关理论 分析结果进行了比较。 4 、展望与本文有关问题的进步研究工作。 北京变通大学硕士论文 第二章大变形有限元法简介 本章介绍了大变形的几何描述方法、非线性方程式组的解法及收敛准 则等；重点介绍了非线性橡胶材料的有限元例式及其数值实现i 掘4 0 1 ，并对 通用有限元软件的特点作了简要介绍。 2 1 大变形的几何描述方法 在大变形理论中，对应力应变的描述有两种形式，l a g r a n g i a n 描述法和 e u l e r 描述法，前者是在初始构形下，用共轭应变对k i r c h h o f f 应力和g r e e n 应变描述，后者是在现实构形下，用共轭应力应变对c a u c h y 应力和e u l e r 应变描述。 在大变形有限元中，由于现实构形为一待求的未知构形，所以只可用 采用以已知构形为参考构形的描述方法，即l a g f a n g i a n 描述法。较为常用 的有两种算法即t l 法和u l 法，前者是以初始构形为参考构形的全 l a g r a n g i a n 法，即在整个分析过程中参考位形保持不变，用共轭应变对 k i r c h h o f f 应力和g r e e n 应变表示有限元例式；而后者是以现实已求解的构 形为参考构形的瞬时l a g r a n g i a 法，所有静力学和运动学的变量参考于每 一载荷或时间步长开始时的位形，即在分析过程中参考位形是不断被更新 的。 t l 法和u l 法在本质上没有区别，求解过程中初始构形始终不变；但 相邻构形却不断发生变化，因此两者在方法上又有所不同： 1 、在t l 格式中，单元刚度阵的积分是在初始构形上进行的，所以不 会因单元变形过大而引起奇异；在u l 格式中，单元刚度阵的积分是在t 时 刻的构形上进行的，所以在变形过大时有可能导致引起奇异。此外，u l 格 式必须在每次变形后修改节点坐标。 第二章大变形有限元法简介 2 、在t l 格式中，r c h h o f f 应力是在初始构形上描述的，应力全量可 由上一时刻的应力与应力增量直接相加得到；而在u l 格式中，c a u c h y 应 力不能由上一- 时刻的应力与应力增量直接相加得到。 3 、在t l 格式中，由于k i r c h h o f f 应力是在初始构形上描述的，所以应 力增率为客观应力增率；而在u l 格式中，由于c a u c h y 应力是在现实构形 上描述的，应力增率不是客观应力增率，在增率型本构方程中需扣除刚体 旋转的影响。 4 、对t l 格式和u l 格式的选取，直接与本构关系的描述密切相关。 如塑性问题，本构方程一般用c a u c h y 应力描述，所以选择u l 格式比较方 便；而粘弹材料这类与时间相关的材料，则因现实构型的不断变化，所以 本构方程一般用k i r c h h o f f 应力描述，所以选择t l 格式较为方便。 5 、在t l 格式中，保留了刚度阵中的所有线性项和非线性项：而在u l 格式中略去了高阶量。另外，在u l 格式中引入j a u m a n n 应力率时出现一与 体积相关的量，由于破坏了刚度阵的对称性，一般也被略去，对于能发牛 较大体积变形的材料，则不再适合。 6 、由于t l 格式和u l 格式分别是在初始构形和现时构形下描述的， 所以，在施加分布载荷时，二者在大小和方向上的含义都完全不同。 由于t l 格式无需略去高阶量及体积改变的相关量，所以有较高的精 度，而且无需修改坐标，应力全量也可以用直接与应力增量相加得到，方 法简便，计算量少。从各方面来看，以用t l 法最为合适。 2 2 非线性有限元法的基础 2 2 1 非线性方程组的解法 随着位移增长，一个有限单元已移动的坐标可以以多种方式改变结构 的刚度。一般来说这类问题总是是非线性的，需要进行迭代获得一个有效 的解。 非线性方程组一般可表示为： 妒( a ) z k ( a ) a + f = 0( 2 2 1 ) 北京变通大学硕+ 论文 由于k ( a ) 依赖于未知量a 本身则无法直接求解，因此必须用迭代法求 解。迭代法有多种算法，最为常见的有五种：直接迭代法、n r 法、m n r 法、切线刚度法和初始刚度法。 1 、直接迭代法 假设n ：a o ，则可以得到迭代h 次的近似解 矿= _ k ( 矿1 ) 】- 1 f ( 2 2 2 ) 若迭代收敛，当九趋于无穷大时，矿趋于真解。 使用直接迭代法时，每次迭代都需要重新计算和形成新的系数矩阵 k ( a ”1 ) ，并对它进行求逆计算。缺陷在于不能保证迭代的收敛性。 2 、n e w t o n - r a 【p h s o n 方法( 简称n r 法) 如果式( 2 2 1 ) 的第月次近似解a “已经得到，为了得到进一步的近似 解a ”1 ，将1 ；f ，( a ”1 ) 表示成在a “附近的仅保留线性的t a y l o r 展开式，并经进 一步计算，可以得到 a a “- ( k ；) 。1 妒” ( 2 2 3 ) 式中； k ；k ，( a “) ( 2 2 4 ) 石d * p k r 。) ( 2 2 5 ) 一般情况下，它具有良好的收敛性，但每次迭代也需要重新形成和求 逆一个新的切线矩阵k ；。 3 、修正的n e w t o n r a p h s o n 方法( 简称m n r 法) 迭代方程可以修正为 a a “一( k 0 ) 。v “ ( 2 2 6 ) 第二章大变形有限元法简介 每次迭代只需求解一相同方程组，计算比较经济。虽然付出的代价式 收敛速度较低，但总体上还是可能合算的。如和加速收敛的方法相结合， 计算效率还可进一步改进。 4 、切线刚度法 用增量法求解非线性方程组，将载荷分成多个增量步，并求得增量型 代数方程 a ”1 = _ 【k ( a 5 ) 】_ 1 f( 2 2 7 ) 切线刚度法的优点是收敛速度较快，但在每次迭代过程中，都需要重 新计算刚度矩阵，所以每步的计算量都比较大。 5 、初始刚度法 与切线刚度法基本相同，也是基于增量代数方程组，只是仅在每增 量步的第一次计算刚度矩阵，而以后每次则无需重新计算刚度矩阵 a a ”1 k 1 ( a o ) r 1 a f( 2 2 8 ) 初始刚度法，每次迭代计算量小，但收敛速度快。 2 2 2 收敛准则 对于大变形闯题，由于研究的橡胶材料具有物理和几何非线性，从而 使得用有限元计算时，会出现许多非线性方程组，如平衡方程组： 妒( u ) 一k ( u ) u r = 0 ( 2 2 9 ) 必须经过多次迭代才能求得精确解。对于相同的问题，如选择不同的 收敛准则，则会宣接影响求解的计算清度和经济性，因此，如何选择合适 的收敛准则也是一个重要问题。常用的三种收敛判断准则为不平衡力的收 敛准则、位移收敛准则及能量收敛准则。对于不平衡力的收敛准则： f o l e r - 1 1 ：1 忙6 ( 2 2 - 1 0 ) 式中： 怕l ：表示向量的范数 m ：表示载荷分级数 北京交通大学硕士论文 n ：表示该载荷分级的迭代次数 6 ：表示收敛精度 2 3 非线性橡胶材料的有限元例式 采用变形前的构形为参考状态的l a g r a n 西a n 法，以该状态的g r e e n 应 变张量的物理分量和k i r c h h o f f 应变张量的物理分量建立起适用于大位移和 大应变的增量形式的有限元计算例式，考虑了不平衡力的收敛精度和橡胶 材料的不可压缩性。有限元法考虑橡胶材料不可压缩的影响，通常可以采 用罚有限元法1 4 “2 1 、混合元法【”。4 6 】及杂交元法【4 7 - 4 9 1 等。混合元法和杂交元 法在有限元分析中引入静水压力作为未知量，罚有限元法把静水压力项通 过罚函数代替。本文以罚有限元法为例给出有限元计算例式。 2 3 1 变形描述 设任一变形体初始时刻岛的构形为r ，x 表示尺上的空间坐标系；t 时 刻的即时构形用r 表示，石表示r 上的空间坐标系，表示位移，妒表示函 数，则： z 一妒( x ，f ) 一x + “( x ，f ) ( 2 3 1 ) 取初始构形作为参考构形，应力、应变分别取第二p i o l a k i t c h h o f f 应 力和柯西一格林应变，分刺用和表示，则内力所做虚功u “为 u “。f r s o e o d r ， ( 2 3 2 ) 式中， 。i o w ( 2 3 3 ) & = 三喝一如) ( 2 3 t4 ) g ! f t 矗毛 ( 2 3 5 ) 第二章大变形有限元法简介 。熹 眩。e ， 式中，一应变能密度函数，6 , i - 黼，d “- - k r o n e c k e r 符号， b 变形梯度。 2 3 2 平衡方程的线性化 根据虚功原理，平衡方程为 6 u 耐= 6 u ( 2 3 7 ) 瓦甲， u “表不外力在虚位移上所作的虚功，5 表示对位移的变化。 由式( 2 3 2 ) 取变分，位移用增量形式表示，则 ( s u j c t t ) = j f 咒可0 6 u j 洲面o a u k 崛+ 矗辟c 矗等+ 等，h 三c f q ko 可o u q + 乃等，d r ( 2 3 8 ) 1 弹2 6 口 d 。一! 竺 侧喊o e u ，刁知分别为初应力张量和材料张量。 外力虚功u “的变分为 5 u “一c r “6 “d r j r ， 5 将式( 2 3 8 ) 一( 2 3 1 1 ) 代入( 2 3 7 ) ， ( 2 3 9 ) ( 2 3 1 1 ) 最终得到增量位移形式的 北京交通大学硕士论文 有限兀万栏： k l 血 。= r 7 一 ，r ( 2 3 ，1 2 ) 式中，h 为载荷的分级数。 戊= 厶 丑f s ) 戤 ( 2 3 1 3 ) k ，k l “k k ( 2 3 1 4 ) 陋l = 厶 岛+ 口j d 8 0 + 吼】媳 ( 2 3 1 5 ) k l 2 厶 b 驴p 】g 戤 ( 2 3 1 6 ) 矩阵【& 】，【吼】， b o a r 的具体形式都可以由式( 2 3 7 ) 导出。 2 3 3 本构关系 第一章已经提到，尽管橡胶材料的应变能函数表达形式不同，一般说 来，总可以用其变形张量的三个不变量来表示，即 w = “，? 2 ，1 3 )( 2 3 1 7 ) 如果材料是不可压缩的，则 ，3 = 1 ( 2 3 1 8 ) 对于不可压缩橡胶材料，在工程上被广泛应用的应变能函数为 ；矽( ，l ，2 ) + j 1 口( ，3 1 ) 2 ( 2 3 1 9 ) 式中，a 为罚因子，缈“，：) 的解析形式用、1 2 的双重傅利叶级数表 示： w ( 6 ，2 ) 。k “一3 ) ( ，。一3 ) “ ( 2 3 2 0 ) 式中，c 为材料常数。( ，1 ，：) 取级数的一阶项，则有： “，2 ) = c l ( 一3 ) + c 2 ( ，2 3 )( 2 3 2 1 ) 第二章大变形有限元法简介 由( 2 3 3 ) 得： s 一里：坐堕+ 塑堕+ 业堕 ( 2 3 2 2 ) 3 一一一十一一+ 一一 、二- 。a e 。|a 1 1 娩j a j 2a e qa 1 3a e b 式中， 参观如 眨3 ， 面0 , 1 2 4 2 ( ，l 岛一呜) ( 2 - 3 2 4 ) 酱以吲 眨，z s ， 式中，唏1 为q 的逆张量。 将( 2 3 2 3 ) ( 2 3 2 5 ) 代入式( 2 3 2 2 ) ，得 s 。2 c 1 6 i j + c 2 ( 1 6 口一g ) + a g 一1 ) i s g ( 2 3 2 6 ) 柯西应力与s 的关系为 5 击s “弓 ( 2 3 2 ，) 由式( 2 3 1 0 ) 可知，求得增量应力砸与增量应变a 的关系，得： = 4 c 2 【6 “6 f 一靠6 + 孙f 2 j 2 ，2 ，- 1 - 1 + 2 1 3 ( 1 9 1 。- 1 瓯1 一嚷1 蛎1 ) 】 ( 2 3 2 8 ) 式中，k 为四阶张量，共8 1 个分量。考虑对称问题，则可转换成四 阶矩阵，便于有限元积分计算，并且含有罚因子吐项，采用降阶积分， 可以避免有限元方程的“自锁”现象。 北京交通大学硕士论文 2 3 4 数值计算步骤 非线性有限元方程( 2 ，3 1 2 ) 计算步骤： ( 1 ) 将整个求解时域分成若干步长 t o = o ，t i ，t 2 ，l 或者凡= 0 ，凡l ，r 2 ， ( 2 ) 时闯步长( f 一一f ，f 。= f + & 】，开始的变量，已，& ，已知， 计算r “，“。 ( 3 ) 建立剐度阵【k f 的各项，并求解缸1 ： 血1 = 【屉 ? “r “一 ， “) ( 2 3 2 9 ) ( 4 ) 计算失衡力，进行平衡迭代： 矿= r ) “一k 降+ ：】1 ( 氏+ a s “溉 ( 2 3 3 0 ) 跏”【k r l 妒4 ( 2 3 3 1 ) 缸一a + & k ” ( 2 。3 3 2 ) 当妒“小于给定的收敛精度示，迭代终止。 ( 5 ) 计算k “= t + a t 时刻的各变量+ 1 ，已+ l ，& 】， ( 6 ) 重复计算( 2 ) ( 5 ) 步，计算下一个时间步长。 2 4 非线性有限元的应用现状分析 在许多重要的实际问题中，例如橡胶类材料受载荷产生大变形，此时， 基于小变形假设的线弹性有限元法已不再适用。因而非线性有限元日益发 展起来。 随着橡胶材料以及高分子材料的广泛应用，非线性有限元取得了很大 进展，但在橡胶材料有限元分析中仍存在许多困难，主要是因为： 第二章大变形有限元法简介 1 、有限变形弹性材料在承受载荷时，同时发牛大位移、大转动和大应变， 材料具有非常复杂的非线性。 2 、通常橡胶材料的体积在变形过程中不变或接近不变，即具有不可压缩 性。对于不可压缩材料，应力张量不能由变形唯一确定，当计算应力张量 时，盛须考虑静水压力。因此给位移型有限元模型分析非线性弹性问题带 来很大困难。 自1 9 6 5 年，h e r r m a n n 发表了关于不可压缩弹性体的变分原理文献，许 多计算工作者给予高度重视，主要是围绕不可压缩性提出不同的分析方法， 包括上文提到的罚有限元、混合元及杂交元等。其后，o d e n 5 0 规! t a y l o r 蚓， s i m o 5 4 1 ，b a t h e l 5 5 】等这些计算力学家的出色工作，推动了非线性有限元的发 展。 随着计算机技术的高度发展，运算速成度不断提出高，用有限元求解 不可压缩弹塑性问题的进展很大，它的数学基础日趋完善。有限元法经过 几十年的发展获得了巨大的成功，这不仅表现在有限元理论已经逐步完善， 还表现在已经开发了大量用之有效的有限元软件，许多大型商用有限元程 序中已经包含了一些非线性弹塑性材料有限元分析，但是能广泛处理橡胶 材料大变形问题的有限元并不多见。 在目前较流行的有限元软件中，使用较广泛的如美国的s a s i 公司的 a n s y s 软件，荷兰t n o 公司的d i a n a 软件，澳大利距的g + d c o m p n t i n g 公司的s t r a n d 软件，a l g o r 公司a l g o r 软件，s t r u c t u r a lr e s e a r c h & a n a l y s i s 公司的c o s m o s m 软件等f 5 “”。这些大型通用有限元分析系统一般有以下特 点： 1 、优秀的人机晃面，完全的c a d 式操作环境和前后处理功能。用户 可以直观生成所要分析的模型，全自动的复杂模型有限元网格生成系统能 按用户的要求来生成有限元离散网格，强大的后处理能按用户的需要产生 各种各样的后期动静态彩色效果等。 2 、强大的分析功能，非富的单元库和求解器。现化的有限元分析系统 几乎都包含了静动态，线性非线性，优化设计，力学，热力学，电磁学及 及耦合问题等分析和求解功能，丰富的单元种类库和材料种类库可以满足 北京交通大学硕士论文 各种要求。 3 、优秀的数据与图形输入输出交互功能。人机对话式的几何与物理信 息输入输出，以及与国际上广泛使用的众多的c a d 系统数据与图形交换功 能，使用厂t 更方便地构造分析模型。 由丁通用有限元程序的广泛使用，力学计算工作逐渐从力学工作者手 里解放出来，工程技术人员对构件进行数值分析变得比较容易。通用有限 元程序的计算结果已经成为各类工业新产品结构设计的可靠性依据。要对 裂纹尖端的位移场和应变场进行准确数值分析，要求数据结果具有足够的 精度和权威性，本文的主要计算工作采用a n s y s 软件进行。 北京交通大学硕士学位论文 第三章裂纹尖端的有限元模型 3 1 扩张区与收缩区旧3 考虑个平面应变裂尖试件，其变形前和变形后的状态分别如图3 1 ( a ) 、( b ) 所示，设俾，e ) 为变形前极坐标，( r ，日) 为变形后极坐标，两者的坐 标原点都在裂尖点。由于裂纹尖端变形十分剧烈，无法用统一的映射函数 来描述其变形，g a 0 1 6 2 】提出将裂尖区域划分为个扩张区和两个收缩区。在 变形之前，裂尖附近的扩张区域e x 非常狭窄，而收缩区域s h 和s h l 贝0 几 乎占据整个裂纹尖端区域；在变形之后，扩张区e x 占据了几乎整个角区， 而收缩区s h 和跚贝1 3 变得非常狭窄。通过分区方法，采用不同的映射函数 描述扩张区与收缩区域的变形，最终得到了理论分析结果。 描述橡胶类材料的应变能函数形式有多种，m o o n e y r i v l i n 型应变能函 数是一种典型的描述橡胶材料本构关系的理论模型，其合理性已被大量实 验证实。若以m o o n e y r i v l i n 型本构关系为基础，分析裂尖场，是否存在扩 张区与收缩区? 裂纹尖端的变形场和应力场怎样? 这是本文要回答的问 题。 第三章裂纹尖端的有限元模型 ff?f 。 l 援p i s h ： i iii 、 ( a ) 变形前( b ) 变形后 图3 1 裂尖尖端的扩张区与收缩区示意图 3 2 有限元模型的建立及求解 应用通用有限元软件a n s y s 建立模型及求解的一般步骤为 1 、根据要求对所提出问题进行建模： 2 、对模型进行有限元网格划分； 3 、根据载荷条件进行加载求解： 3 2 1 i 型裂纹的有限元模型 m o o n c y r i v l i n 应变能函数在橡胶类材料中具有代表性，本文用此类材 料建立有限元分析模型。 图3 2 为平面应变l 型裂纹问题。 i l 、，2 、1 3 为应变张量的三个独立 不变量，假定材料不可压缩，l j 1 3 ；1 。m o o n e y r i v l i n 应变能函数为： w t j ( ，l 一3 ) + c 2 ( t 一3 )( 3 2 1 ) 式中，w 为单位体积的应变能，c ，、c ：为材料常数。 北京交通大学硕+ 学位论文 a n s y s 中m o o n e y r i v l i n 本构关系中两参数材料常数取i d i n g 等人 6 3 】 曾给出材料常数， c i = 0 1 4 8 8 6 2 m p a m2 ，c 2 = 0 1 0 8 4 9 6 m p a m 2 。 模型的几何尺寸见图3 2 ，图中村。0 4 m ，厚度为l m ，在中心有一长 2 a = 0 0 4 8 m 的对称裂纹。 p 图3 2 i 型裂纹模型 3 2 2 有限元网格划分 一z 对橡胶类材料裂纹尖端场进行有限元网格划分，要注意以下要点： 1 、裂纹尖端区域没有一个明确的概念，一般分析区域都非常小。分析 区域必须有足够精细的网格划分，但网格过小，其余周围区域就存在过渡 困难。因此，区域网格密度变得特别地重要，要求得一个精确的分析，必 须提供足够分析应力的网格密度，同时考虑区域之间的连接。本文采用了 等比单元划分。 第三章裂纹尖端的有限元模型 2 、为保证计算结果的可靠性，应尽量避免低劣的单元形状，例如防 止过大的纵横比，过度的顶角等。同时，必须注意防止网格单元的原始形 状和变形后形状发生过度扭曲，必要时对某些部位迸行的人工网格划分以 产牛合理的计算结果。 考虑对称性，以降低计算时间，取图3 2 中的四分之一计算，有限元网 格划分见图3 3 ( a ) ，a h 、o f 为对称边，o a 为裂纹边缘。图中，a h = a f = 么， o a = a 。 考虑到裂尖附近可能存在应力奇异性，以裂尖点为圆心，r t o 0 2 4 m 为 半径划出一个半圆形区域，加密剖分，见图3 3 ( b ) 。先将半圆沿环向等分 为1 8 个区域，每个区域的角度为1 0 0 ，然后进步细分0 0 一6 0 0 区域，见图 3 3 ( c ) 。裂尖处的网格划分采用了等比划分，公比的大小可以由计算需要 的精度来确定，适当采用公比，可以达到加密尖端网格，又能控制单元数 多少，以适合计算量大小的需要，降低求解的时间成本。考虑以上因素， 沿半圆的径向布置了6 0 个节点，由密到疏等比划分，公比为1 ：1 1 7 ，共 6 0 圈有限元网格。裂尖处的第1 圈单元由于两个节点合二为一，即采用三 角形单元。单元沿径向的长度为3 4 6 9 1 x 1 0 m ，其它区域使用了四边形单元。 采用a n s y s 中超弹性h y p e r 7 4 单元，在本文模型的有限元网格划分中共 有3 9 7 5 个单元，总节点数为1 2 1 9 0 。 北京交通大学硕士学位论文 h g ( a ) 结构 i 瓣瑟 ( b ) 裂纹模型局部 ，2 3 第三章裂纹尖端的有限元模型 一一 溢熬 ( c ) 裂纹尖端附近 图3 3 模型有限元网格剖分 3 2 3 载荷加载 大应变处理对一个单元经历的总旋度或应变没有理论限制。然而，应 限制应变增量以保持糙度，因此使用小的载荷增量非常重要。总载荷应当 被分成几个较小的步，缓慢加载，通过逐渐地施加载荷来提高保守系统的 收敛特性，从而使所要求的n e w t o n r a p h s o n