（光学工程专业论文）marc软件辅助表面接触特性的研究及在汽车结构分析中的应用.pdf

、 ， iifm t h ed e g r e eo fm a s t e ro fp h i l o s o p h y f o r m a r cs o f t w e a ra i d e ds u r f a c ec o n t a c t r e s e a r c ha n di t sa p p l i c a t i o ni n a u t o m o b i l es t r u c t u r e a n a l y s i s s c h o o l d e p a r t m e n t ：s c h o o lo fa u t o m o b i l e d i s c i p l i n e ：m e c h a n i c a le n g i n e e r i n g 一 一一 m a j o r ：v e h i c l ee n g i n e e r i n g 一 c a n d i d a t e ：c h e n 物b a o s u p e r v i s o r ：p r o f c h a oz h a n g m a y , 2 0 0 6 学位论文版权使用授权书 了解同济大学关于收集、保存、使用学位论文的规定， 同意如下各项内容：按照学校要求提交学位论文的印刷本和电子版 本；学校有权保存学位论文的印刷本和电子版，并采用影印、缩印、 扫描、数字化或其它手段保存论文；学校有权提供目录检索以及提供 本学位论文全文或者部分的阅览服务；学校有权按有关规定向国家有 关部门或者机构送交论文的复印件和电子版；在不以赢利为目的的前 提下，学校可以适当复制论文的部分或全部内容用于学术活动。 学位论文作者签名：白皂臣专 6 年6 月f 复e t 经指导教师同意，本学位论文属于保密，在县# 年解密后适用 本授权书。 指导教师签名：丐茛争 学位论文作者签名：鱼睦芦号 勿奶年月1 2e t 2 6 年6 月f 文日 指导下，进行 研究工作所取得的成果。除文中已经注明引用的内容外，本学位论文 的研究成果不包含任何他人创作的、已公开发表或者没有公开发表的 作品的内容。对本论文所涉及的研究工作做出贡献的其他个人和集 体，均已在文中以明确方式标明。本学位论文原创性声明的法律责任 由本人承担。 签名：毙p 寄 o 6 年6 月d 日 - 研究将有利于改 占有极其重要的 作用。有限元法由于其对复杂几何构形适应性和适合计算机实现的高效性等特 点，成为研究表面接触特性的常用工具之一。 本文前半部分利用有限元软件m a r c 对粗糙表面接触特性进行了研究。建 立了考虑真实表面形貌，随温度变化的材料特性以及采用完全热机耦合的接触 模型，并首次对于不同的滑压模式进行了分析。在比较了各工况接触状态，温 度场以及应力应变场的分布等计算结果后，得出接触过程受滑压模式影响的结 论。文章的后半部分利用多体动力学分析软件e x c i t e 对汽车关键零部件之一 轴承的润滑和接触特性进行了分析，并进行了疲劳寿命计算，得到了与实验相 一致的结果。 最后，关于进一步工作的方向进行了简要的讨论。 关键词：接触，有限元，热机耦合，疲劳寿命 a u t o m o b i l e ，t h es t u d yo ft h e i r ep e r f o r m a n c eo ft h ev e h i c l e p a r t i c u l a r l y , t h es t u d yo fs u r f a c ec o n t a c tp r o b l e mp l a y sa l li m p o r t a n tr o l ei nt h e t r i b o l o g yr e s e a r c h f i n i t ee l e m e n tm e t h o db e c o m em o r ea n dm o r ep o p u l a ri ns u c h r e s e a r c hf i e l dd u et oi t sa d a p t a b i l i t yt oc o m p l e xg e o m e t r ya n de f f i c i e n c yw h i c hc a n b ea c h i e v e db yc o m p u t e re a s i l y t h i st h e s i sa r ed i v i d e di n t ot w o p a r t s i nt h ef i r s tp a r tw e u s et h ef e ms o f t w e a r m s c m a r ct os o l v et h ec o n t a c tp r o b l e mb e t w e e nar i g i df i a ts u r f a c ea n da d e f o r m a b l er o u g hs u r f a c e t h ec o n t a c tm o d e lc o n s i d e r e dr e a ls u r f a c et o p o l o g y , t a k e t h em a t e r i a l p r o p e r t y a s t e m p e r a t u r e d e p e n d e n t ，a n d u s e sa f u l l yc o u p l e d t h e r m o m e c h a n i c a lf i n i t ee l e m e n tm o d e l t h em o d e lc o m p a r et h et w od i f f e r e n ti n i t i a l a l g o r i t h mf o rt h ef i r s tt i m e r e s u l t ss h o wt h a tt h es t e a d yc o n t a c ts t a t u sd e p e n d s s i g n i f i c a n t l ya n dc o m p r e h e n s i v e l y o ni n i t i a la l g o r i t h m ，m a t e r i a lp r o p e r t ya n d o p e r a t i o nc o n d i t i o n s i nt h es e c o n dp a r tw ea n a l s y st h el u b r i c a t i o na n dc o n t a c t p r o p e r t yo ft h eb e a r i n gw h i c hi so n eo ft h ek e yp a r t si na na u t o m o b i l eb yu s i n g m u l t i b o d yd y n a m i c ss o f t w e a re x c i t e t h ef a t i g u el i f ea n a l s y sa l s op e r f o m e da n d r e s u l t ss h o w st h ec o n s i s t e n c yw i t ht h ee x p e r i m e n t f i n a l l y , t h ep r o b l e m sr e q u i r i n gf u r t h e r s t u d i e sa r ed i s c u s s e d k e yw o r d s ：c o n t a c t ，f e m ，t h e r m o - m e c h a n i c a lc o u p l e d ，f a t i g u el i f e 1 2 表面接触问题的发展1 1 2 1 光滑弹性接触2 1 2 2 考虑粗糙度的表面接触。2 1 2 3 摩擦引起温度变化的接触情况3 1 2 4 同时考虑粗糙度和摩擦生热时的表面接触情况4 1 3 本文的研究目的和意义4 第2 章有限单元法7 2 1 弹性力学问题的有限元方法7 2 1 1 弹性力学的基本方程7 2 1 2 弹性力学的有限元格式1 l 2 2 接触问题的有限元法1 4 2 2 1 接触问题求解的一般过程1 4 2 2 2 接触界面的定解条件和校核条件1 5 2 4 3 引入接触约束的方法。1 7 2 3 热机耦合问题的有限元法1 9 第3 章接触模型的建立2 3 3 1 几何模型的建立2 3 3 2 有限元模型的建立2 7 3 2 1 有限元网格的划分。2 7 3 2 2 材料参数的设定2 8 3 2 3 边界条件的设定3 4 i i i 目录 3 2 4 接触定义3 6 3 2 5 载荷工况的设定及任务提交3 9 3 2 6 模型并行计算4 0 3 3 磨损模型建立4 1 3 3 1 子程序格式4 2 3 3 2 子程序思路4 3 3 3 3 子程序的调用4 4 第4 章计算结果分析4 7 4 1 接触程度控制模式4 7 4 1 1 接触状态结果分析。4 7 4 1 2 温度结果分析5 3 4 1 3 应力结果分析6 1 4 1 4 磨损模型结果。7 0 4 1 5 结果分析总结：。7 1 4 2 接触载荷控制模式7 l 4 2 1 接触状态结果分析7 2 4 2 2 温度结果分析7 4 4 2 3 应力结果分析8 0 4 1 4 磨损模型结果8 4 4 2 5 结果分析总结：。8 6 第5 章轴承润滑分析8 7 5 1 发动机轴承润滑分析8 7 5 1 1 发动机运行工况简述。8 7 5 1 2e h d 计算模拟8 8 5 1 3 计算结果分析9 3 5 2 试验台轴承润滑计算9 6 5 2 1 蓝宝石试验台简述9 6 i v 目录 9 8 1 ( ) ( ) 1 0 5 1 0 5 1 0 ：； 6 1 2 载荷谱的分解1 0 6 6 1 3 线性累积损伤理论1 0 6 6 2 疲劳分析1 0 7 6 2 1 节点力转化1 0 7 6 2 2 有限元计算1 0 8 6 2 3 疲劳寿命计算1 1 0 第7 章结论与展望1 1 3 7 1 结论儿3 7 2 进一步工作的方向1 1 3 致谢儿5 参考文献1 1 7 个人简历在读期间发表的学术论文与研究成果1 1 9 v 摩擦学是有关摩擦、磨损与润滑科学的总称。它是研究在摩擦与磨损过程 中两个相对运动表面之间的相互作用、变化及其有关的理论与实践的一门学科。 由摩擦引起的磨损、润滑、材料与能源消耗等一系列摩擦学问题普遍存在并对 社会、经济的发展产生了巨大影响。由于摩擦学科学所涉及的问题与节约能源、 节约材料、减少磨损、提高资源利用率和保护环境等密切相关，一直受到科技 界的高度重视。 摩擦学问题中各种因素往往错综复杂，涉及到多门学科，与流体力学、固 体力学、流变学、热物理、应用数学、材料科学、物理化学，以及化学和物理 学等学科的内容有关。因此摩擦学是一门多学科交叉的边缘科学，时至今日， 这门科学仍然处在发展之中，随着理论研究的日益深入和实验技术的日益先进， 摩擦学研究方法的发展由宏观进入微观；由定性进入定量；由静态进入动态； 以及由单一学科角度的分析进入多学科的综合研究。 在摩擦学研究中表面接触的研究分析占有极其重要的地位。通过在不同的 运动条件及材料性质情况下对接触状况、接触应力等的分析可以得到一些相关 规律，其为分析零件接触面的磨损、腐蚀、疲劳失效提供了依据，也为零件的 进一步优化设计提供了平台。近几年计算机辅助分析受到日益关注，有限元分 析软件由于其计算模块、计算方法日趋完善，在接触分析中往往能用较短的时 间得到相当准确的分析结果，正在为越来越多的科研人员、工程师们所使用。 1 2 表面接触问题的发展 自从h e r t z 在1 8 8 2 年解决了理想光滑弹性体的接触力学问题后，接触力学 有了很大的发展。比较突出的是1 9 8 5 年，j o h n s o n 阐明了半无限弹性体变形特 点，并获得了压力分布和次表面层的应力分布的解析解。这些结果很有价值， 为设计某些机械零件，如滚动轴承、齿轮、密封件和机器关节等提供了理论依 第1 章引言 据 1 2 1 光滑弹性接触 h e r t z 接触i 司题的假设： 1 互相接触物体是光滑和均质的。 删 2 在接触区域仅有弹性变形发生。 3 接触力垂直于接合面。 。牡 ， 图1 1h e r t z 接触问题的描述 早期的接触力学研究重点基本都放在理想光滑表面的微接触变形行为上。 第一个关于圆柱接触的弹性分析是p o r t i s k y 在1 9 5 0 年发表的，得到了在受到正 压力和切向牵引力作用时的接触应力和变形。在1 9 6 6 年h a m i l t o n 和g o o d m a n 推导出了弹性圆柱接触在受到半球状正压力和按比例分布的剪切牵引力作用时 的表面和次表面层的应力。在处理弹性固体接触的难点是这些弹性固体的形状 是多样的，并不是h e r z 接触理论中所用的椭圆体，但是使用了数值技术后这一 难点就被克服了；一般的方法是把复杂的解析表达式转化为线性代数的表达式， 很容易用数字计算机对其进行求解。这样，在1 9 7 9 和1 9 8 0 年，h a r t n e r 采用一 种数值迭代方法，获得了线接触问题的接触压力分布。1 9 8 6 年，s a c k f i e l d 和h i l l s 给出了不同形状的压头对压力分布和次表面层应力的影响，得到在比圆柱体平 展压头情况下的相当均匀的压力分布和非h e r t z 接触次表面层的应力场和应力松 弛。 1 2 2 考虑粗糙度的表面接触 在九十年代以前，对粗糙接触的研究相对于光滑接触要少很多，这主要是 由于粗糙接触的相对复杂的数学模型造成的。g r e e n w o o d 和w i l l i a m s o n 【1 】在 1 9 6 6 年提出了一种方法来估计表面粗糙度，基于对许多工程粗糙表面的观察和 2 概率密度为 率半径r 。 对于某一微凸峰，设在z 。一d 处第i 个微凸峰发生接触，按两球面接触理论， 其接触区圆半径为口；一( r 6 ；) 1 他，接触面积为a a = 积6 ，接触载荷为 e = 必e r “2 瓯3 心，其中瓯为第i 个微凸峰与平面接触时的相对法向接近量，e 为等效弹性模量，根据概率理论任一高度z 。的微凸峰发生接触的概率为 p r o b z , d t - c ( z ，比 ( 1 2 1 ) 如果名义接触面积为a ，则发生了接触的微凸峰总数的数学期望值为 n 一棚c 驴0 ，) a z ， ( 1 2 2 ) 总的接触面积的期望值为 4 = 聊肷c ( z s - d n z ，比 ( 1 2 3 ) 总的接触载荷的期望值为 f = 伽州r 佗f ( 乙一d ) 3 佗妒q ， ( 1 2 24 ) 通过如上提出表面参数，如粗糙峰高度、曲率和粗糙峰斜率的概率密度函 数，来假设一种理想的接触几何形状以求得解析解。但是，这种概率的模型是 半经验的，需要假设粗糙峰参数的概率分布，并且没有包含粗糙峰的相互作用。 2 0 0 5 年，m o l i n a r i 和r o b b i n s 【2 】2 等建立了真实粗糙表面的有限元接触模型， 模型考虑了粗糙弹塑性表面在载荷下和刚性平面的接触情况，并分别采用了4 节点和1 0 节点两种4 面体有限元单元，得出的结果表明实际接触面积a 和载荷 成线性增长关系，且其增长率随着屈服极限的增加而减小；并得出了由于塑性 的影响使接触压力和接触区域和弹性的情况相比有了质的区别。 1 2 3 摩擦引起温度变化的接触情况 自从b o l k 【3 】和j a e g e r 【4 】的研究后，对于摩擦引起温度变化的研究就成为 摩擦学研究主要主题之一。1 9 9 4 年，t i a n 和k e n n e d y 【5 】分析了在不同的p e c l e t 3 第1 章引言 数和给定不同形状热源情况下的表面温升。1 9 9 8 年q i u 和c h e n g 研究了混合润 滑接触中的温升问题，他们采用一个独立的等温接触模型来计算压力分布，以 决定热源。g a o 等提出了粗糙表面的瞬态f l a s ht e m p e r a t u r e 模型，用傅立叶变换 技术求热分割问题，同时得到压力分布。2 0 0 3 年，n y e 和k k o m v o p o u l o s 【6 】6 给 出了在受到热机载荷时弹塑性涂层介质的三维有限元分析，模型采用了弹性球 面和弹塑性光滑平面的滑压接触模型，并考虑了完全的热机耦合，模型通过和 弹性均匀半平面的解析解相比，在温度和应力分布上有着良好的一致性，结果 表明摩擦剪切牵引力和由于发热引起的载荷促使了应力集中和塑性变形，尤其 是在涂层薄且热传导率低的时候。但模型没有考虑表面粗糙度的影响。 1 2 4 同时考虑粗糙度和摩擦生热时的表面接触情况 两个物体发生接触运动时粗糙度与摩擦是同时存在的，但是由于模型的复 杂化和数值技术的难度，把两者结合起来的研究不是很多，而且以前的研究对 模型作了太多的简化。一些关于径向轴承混合润滑的研究分析了热传导和粗糙钿 峰的接触。在对接触问题的热分析中主要集中在粗糙度的f l a s ht e m p e r a t u r e 上。 在分析接触问题的热弹性分析时主要集中在赫兹接触和单个粗糙锋接触上或者 热不稳定上。很少有文章同时结合了热弹分析和粗糙接触时粗糙度在摩擦生热 过程中的变形分析。 近些年来，对粗糙表面接触时的摩擦发热问题的研究逐渐完善起来。1 9 9 9 年，n e d e r 和v a r a d i 等【7 】分别采用了数值方法和有限元法分析了真实粗糙表面 在接触时的温度分布，其中有限元采用了8 节点六面体单元，通过计算分别得 到了低速( 1 m s ) 和高速( 1 0 m s ) 情况下不同接触材料的接触温度分布，结果 表明接触温度将随滑动速度的增加而增加，最后还考虑了粗糙度方向对结果的 影响，得出结果表明粗糙度方向平行于滑动方向情况下的接触温度中是高于粗 糙度方向垂直于滑动方向的情况。但模型采用的是室温情况下的材料参数，没 有考虑其随温度的变化。 1 3 本文的研究目的和意义 汽车中存在着大量的运动摩擦副，如气缸一活塞环、轴承一轴瓦以及制动 4 第1 章引言 器等等，他们都不可避免的要产生磨损问题。要研究这些运动副的磨损情况， 破坏失效以及动力学特性等内容就必须对两表面的接触特性进行分析，建立接 触程度和接触力之间的关系，将这些研究结果应用在进一步的分析中。 从近几十年表面接触状态研究的发展过程来看，解析法由于其考虑函数求 解的难易程度对模型的假设过多，导致计算结果和真实情况出入较大，从而数 值解法由于其对模型的假设较少则被更多的采用。近年来，有限元法由于其对 复杂几何构形和各种物理问题的适应性、严格的数学理论基础和适合计算机实 现的高效性等特点，被越来越多的应用于表面接触问题的研究中，同时随着计 算机技术的发展和大型方程组解法研究的进展，使得求解几十万甚至几百万个 自由度的大型有限元问题成为可能。一批由专业软件公司研制的大型通用商业 软件( 如n a s t r a n ，a s 队s 觚a n s y s ，m a r c ，a b a q u s ，j i f e x 等) 已为工 程技术界广泛应用，并成为a 旧洲系统不可缺少的组成部分。 尽管不少国内外文献利用有限元法研究表面接触问题，但其中多数存在如 下问题： ( 1 ) 没有采用完全的热机耦合，即采用温度场和热应力等的解耦方法。接触 表面产生的摩擦热将对接触体的温度分布和热变形产生影响，同样热变形又将 反过来影响两表面的接触压力和面积从而对摩擦热产生影响。所以这两个相互 影响的因素必须同时考虑才更为合理，而并不应该分开进行迭代( 每一迭代步 中的摩擦热和温度分布是通过恒定的接触压力和接触面积得到的) 。 ( 2 ) 没有考虑材料本构关系随温度变化的影响。大多研究都把材料特性设为 常数，而试验表明，随着温度的提高，材料将有明显的软化现象，所以必须考 虑材料的本构关系随温度变化的关系。 ( 3 ) 对塑性材料没有考虑变型路径的影响。 5 第1 章引言 _ o 0 2o 0 0 2o 0 4o 0 6o o 1 0o 1 20 1 40 ，1 60 1 8o 2 0 s t r a i n 图1 2l v - 1 2 铝在不同温度下的应力一应变曲线 由于以上缺陷的存在，目前的接触模型并不能真实地模拟运动表面的接触 过程。 本文分为两个部分：第一部分将利用u g 的二次开发接口编写c 程序，建立 有粗糙表面的接触体的三维几何模型，转化成通用数据格式读入非线性有限元 分析软件m s c m a r c 进行网格划分、计算以及结果处理。利用m s c m a r c 的热机耦合分析模块以及先进的接触算法，实现了应力应变场和温度场的同时 迭代；通过m s c m a r c 为用户提供完善的t a b l e 功能( 即各输入参数可建立 函数关系) ，可以建立函数表来实现材料本构关系以及材料特性随温度的变化趋 势，得到合理的接触模型，提高计算精度；计算得到不同的滑压模式对应力应 变场和温度场的影响。对这些影响因素的正确考虑使对表面接触状态的描述更 为合理和真实。然后对接触模型在不同滑压模式下的计算结果进行了比较。第 二部分则是利用多体动力学分析软件a v l e x c i t e 对试验台和发动机轴承的油 膜压力和接触压力等润滑特性进行了分析，并利用其结果进一步进行的疲劳寿 命计算。 6 黜 瑚 黜 蕈耋 稻 季堇 董若 锕 。 问题，人们可以给出它们 微分方程) 和相应的定解 条件。对于大多数问题，不可能用解析方法求出精确解，只能采用数值方法求 解，2 0 世纪6 0 年代以来，随着电子计算机的出现，一种偏微分方程的数值分析 方法一有限单元法得到了迅速发展。 从有限元方法建立的途径来看，它不是直接从问题的微分方程和相应的定 解条件出发，而是从与其等效的积分形式出发。下面对有限元的基本理论做个 简单的介绍。 2 1 弹性力学问题的有限元方法 2 1 1 弹性力学的基本方程 弹性体在载荷作用下，体内任意一点的应力状态可以由6 个应力分量吒， 仃) ，仃：，z 矽，z 声，z 。来表示。其中吒，仃y ，仃：为正应力；z 可，f 弦，f 。为 剪应力。用矩阵表示为 ；b ， 0 y o z 气搿气y z 7 z 。】r 工 y z y 譬 茗办巩办锄和缸 第2 章有限单元法 图2 1 应力分量 跖2 三 2 l ， w 】r 弹性体内任一点的应变可以由6 个应变分量f ，占，：，矽，) ，弦，) ，荔表 示。其中，y ，f ：为正应变；) ，掣，) ，声，y 荔为剪应变。它的矩阵形式为 工 y f z ，捌 y 妒 y 。 ；l ， y zy 碍 弘 ，。】r 1 半衡方；i 茔 弹性体v 域内任一点坐标轴x ，y ，z 方向的平衡方程为， 堡+ 笠+ 监+ 7 。0 毗 秒 a z “ 鲁+ 鲁+ 鲁+ 7 ，- 0 汜， a x酗a z 3 堡+ 笠+ 竺+ 7 ，。0 融 谚 a z 其矩阵形式为 8 问题，人们可以给出它们 微分方程) 和相应的定解 条件。对于大多数问题，不可能用解析方法求出精确解，只能采用数值方法求 解，2 0 世纪6 0 年代以来，随着电子计算机的出现，一种偏微分方程的数值分析 方法有限单元法得到了迅速发展。 从有限元方法建立的途径来看，它不是直接从问题的微分方程和相应的定 解条件出发，而是从与其等效的积分形式出发。下面对有限元的基本理论做个 简单的介绍。 2 1 弹性力学问题的有限元方法 2 1 1 弹性力学的基本方程 弹性体在载荷作用下，体内任意一点的应力状态可以由6 个应力分量吼， 仃y ，仃：，z 弦，z 。来表示。其中吒，仃) ，仃：为正应力；z 叫，f 弦，f 。为 剪应力。用矩阵表示为 ；ka y o z t w 1 y z 7 吃】r 工 y z y 2 x以西办锄和缸 第2 章有限单元法 图2 1 应力分量 “2 三】2 k ， w 】r 坤 弹性体内任一点的应变可以由6 个应变分量s ，s ：，y 叫，) ，弦，y 。表 示。其中，s ，f ：为正应变；) ，秒，) ，弦，) ，。为剪应变。它的矩阵形式为 j y g z y 叫 y 粥 y 。 ；l 。 ， ： ) ，秒) ，曙 7 ，：。】7 1 平衡万程 弹性体v 域内任一点坐标轴x ，y ，z 方向的平衡方程为， 堡+ 生+ 堡+ 7 ，。o 瓠 时d z 笠+ 堕+ 堕+ 7 。- - 0 ( 2 1 1 ) 融 谚d z j 堡+ 笠+ 堡+ 7 ，；o 觑 砂 弛 其矩阵形式为 8 ( 2 1 2 ) ( 2 1 3 ) 7 是体积向量，7 ， _ ，7 y7 ： r 其张量形式为 ，j + 厂j = 0 ( i ；l 2 ，3 ) ( 2 1 4 ) 2 几何方程应变位移关系 在微小位移和微小变形的情况下，略去位移导数的高次幂，则应变向量和 位移向量间的几何关系为 o uo v o w x 2 一o x ) ，2 一o y f z o z ( 2 1 5 ) y 叫。面+ 面y 声2i + 面) ，a 2 i + i 几何方程的矩阵形式为 e = l u ( 2 1 6 ) 其中l 为微分算子，即 l = 4 t( 2 1 7 ) 几何方程的张量形式为 勺= 要g u + u j i ) ( “。1 , 2 ，3 ) ( 2 1 8 ) 3 物理方程应力应变关系 弹性力学中对于各向同性的线弹性材料， 阵表示为 盯一d e 其中 9 应力通过应变的表达式可以用矩 ( 2 1 9 ) 第2 章有限单元法 。= 踹 z ， 1 一u 1 对 2 ， 1 一 u 1 一 1 称 0 o 0 1 一动 。2。(。1。-。v) d 称为弹性矩阵，它完全取决于弹性体材料的弹性模量e 和泊松比矽 物理方程的张量形式可以写成 o - 玎一2 g e 玎+ a 6 盯胜一d o , r e 盯( f ，j ，k ，z 一1 , 2 ，3 ) 其中d 埘称为弹性张量， d m = 2 g 6 饿6 # + 入6 q 6 h 铲鹰弼g = 硐e 4 力的边界条件 ( 2 1 1 0 ) ( 2 1 1 2 ) a 一e v ) ( 2 1 1 3 ( 1 + v ) ( 1 - 2 v ) 。)j 弹性体在边界上单位面积的内力，弓，t ，在边界s 。上已知弹性体单位 面积上作用的面积力为z ，乃，丁：，根据平衡有 t = r x ，l = 丁，t = 丁z ( 2 1 1 4 ) 设边界外法线的方向余弦为n x , 以，万：，则边界上的内力为 _ 1 工= _ ，z x 仃工+ ，z v f 惯+ n z f 硝 t n x f 矽+ 万y 仃y + 以z f 刁， ( 2 1 1 5 ) l = 咒j f + ，l y f 弦+ n z 盯z 写成矩阵形式为 t 一，l 仃一t 1 0 ( 2 1 1 6 ) 秒一u o o o o o丝舡 生计 o o o o丝舡 1 7 ) 效积 值) ， 1 8 ) = 0 。 可以 得到 正觑r ：y i j , j d v - 0 丢g 可+ “，j ) c y i j d v + 正。d u r a o n j d s ( 2 1 1 9 ) 通过几何方程( 2 1 8 ) 并将上式代回( 2 1 1 8 ) 式，就得到了它经分部积分后的 “弱”形式 上【- 嘞+ 踟j 厂l 弦y + 摩t f d s = 0 ( 2 1 2 0 ) 上式体积分中的第一项式变形体在虚应变上所做之功，即内力的虚功：体积分 中的第二项及面积分分别表示体积力和面积力在虚位移上所做之功，即外力的 虚功。外力的虚功和内力的虚功总和为零，这就是虚功原理 虚位移原理的力学意义是：如果力系是平衡的，则它们在虚位移和虚应变 上所做之功的总和为零；反之，如果力系在虚位移和虚应变上所做之功的总和 为零，则它们一定是满足平衡的。所以虚位移原理表述了力系平衡的充分必要 条件。其矩阵形式为 正协t 一d u t 厂砂一肛t 嬲= 0 ( 2 1 2 1 ) 有限元法正是通过求解式( 2 1 2 1 ) ( 平衡方程的等效积分弱形式) 来得到节点 位移向量的。 2 1 2 弹性力学的有限元格式 为简单起见，以平面问题3 节点三角形单元为例说明弹性力学问题有限元 分析的表达格式和一般步骤。 第2 章有限单元法 典型的3 节点三角形单元节点编码为f ，_ ，m ，每个节点有两个位移分量 如图2 2 所示。 图2 2 三角形单元应力分量 每个单元有6 个节点位移即6 个节点自由度，其向量为 a e 2 k f 屹“，j “。v 。j r 一 在有限元方法中，为了建立单元节点力和单元节点位移的关系，必须先假 设一个单元位移模式，以此来得到单元内任一点的位移和节点位移的关系。虽 然整个弹性体内的位移函数很难假设，但单元较小，在每个单元内部可以用比 较简单的函数来近似表达单元内的真实位移。再将各单元的位移函数联合起来 就可以逼近整体的位移函数了。 3 节点三角形单元的位移模式选取一次多项式 “展+ , a 2 x + 3 y ( 2 1 2 2 ) v = | 3 4 + p s x + p b y 岛凡是待定系数，将节点f ，j ，m ，的坐标和位移分别代入式( 2 1 2 2 ) 得到6 个代数方程，则局一成可由节点坐标和位移表示，将求得的表达式代回 式( 2 1 2 2 ) ，可将位移函数表示成节点位移的函数，即 “= n i “f + n j u j + n r n u r a ( 2 1 2 3 ) v = n 以+ ，1 ，+ n r a v 肼 其中n ，。称为单元插值函数或形函数，它是坐标x ，y 的一次函数， 取决于单元的3 个节点坐标，( 2 1 2 3 ) 式的矩阵形式是 1 2 = 讹。( 2 1 2 4 ) 称为插值函数矩阵或形函数矩阵，a 。称为单元节点位移列阵。 由几何方程( 2 1 6 ) 可得到单元应变和节点力的关系为 e - - l uil n a e = b a e 曰称为应变矩阵，l 是平面问题的微分算子。 由物理方程( 2 1 9 ) 可得到单元应力和节点力的关系为 仃：d e ；d b a 。鼬。 ( 2 1 2 5 ) ( 2 1 2 6 ) 其中s 称为应力矩阵。 由弹性力学的虚功原理我们可以建立平衡力系中外力与应力应变的关系， 由于有限元法假设单元只在节点上传递作用力，所以体积力厂和边界力z 都需通 过静力等效原则转化为等效节点载荷，同时应力可以通过函数由节点位移表示， 我们可以得到节点力和节点位移的关系，由虚位移方程( 2 1 2 1 ) 并代入等效节 点载荷f 8 得 f 6 t u d v 一6 a8 1 f 8 ：0 ，y 式中，。分别是单元内任一点虚应变和节点虚位移， 应变全部表示为节点位移的函数，得 面8 1 f 8 ；面。1f b t d b d v a 。 ( 2 1 2 7 ) 将( 2 1 2 7 ) 式中的应力 ( 2 1 2 8 ) 由于缸8 是任意满足约束条件得节点位移得变分，所以式( 2 1 2 8 ) 等价于 f 。：f b t 删玩8 ；k 。a 。 ( 2 1 2 9 ) j y 其中k 8 称为单元刚度矩阵，它反映了节点力和节点位移得关系。 根据节点位移相同，节点受力平衡这两个原则将所有方程联立起来，将离 散的单元联合成一个整体。最后将得到一个由总体结构刚度矩阵( 由单元刚度 矩阵集成) 、单元节点位移以及节点外载荷组成的总体刚度方程 1 3 所n晰以如蛎 第2 章有限单元法 k a ；f( 2 1 3 0 ) k 和f 分别称为总体结构刚度矩阵和结构节点载荷列阵，a 为结构节点自由度列 阵，这样平面弹性力学问题就转化为求解一个线性代数方程组，通过迭代等数 值计算方法可以求出各节点的位移从而进一步得到各单元的其他信息。 以上表述得便是弹性力学的有限元求解的一般原理和步骤。由于总体结构 刚度矩阵的奇异性，需要引入位移边界条件对总刚进行处理，以消除其奇异性， 一般采用乘大数法，这样不改变矩阵结构且适用于非零位移约束。这里不作详 述。 2 2 接触问题的有限元法 接触问题的有限元法常常同时涉及三种非线性，即除材料非线性和几何非 线性以外，还有接触界面的非线性，这是接触过程所特有的。接触界面的非线 性来源于两个方面： 矗 ( 1 )接触界面的区域大小和相互位置以及接触状态不仅是事先未知的， 而且是随时间变化的，需要在求解过程中确定。 ( 2 )接触条件的非线性。接触条件的内容包括：接触物体不可互相侵 入：接触力的法向分量只能是压力；切向接触的摩擦条件。这 些条件的特点是单边性的不等式约束，具有强烈的非线性。 接触界面的事先未知性和接触条件的不等式约束决定了接触分析过程中需要经 常插入接触界面的搜寻步骤。 2 2 1 接触问题求解的一般过程 接触过程通常是依赖于时间，并伴随着材料非线性和几何非线性的演化过 程，特别是接触界面的区域和形状以及接触界面上运动学和动力学的状态也是 事先未知的。这些特点决定了接触问题通常使用增量方法法，并且需要采用试 探校核的迭代方法求解。每一增量步的试探校核过程可一般性地表述 如下。 ( 1 ) 根据前一步的结果核本步给定的载荷条件，通过接触条件的检查和搜索， 1 4 培 第2 章有限单元法 假设此步第1 此迭代求解时的接触面的区域和状态( 这里是指物体a 和b 在接 触界面上有无相对滑动。无相对滑动的接触状态称为“粘结”，有相对滑动的接 触状态称为“滑动 ) 。 ( 2 ) 根据上述关于接触面区域和状态所做的假设，对于接触面上的每一点， 将运动学或动力学上的不等式约束改为等式约束作为定解条件引入方程并进行 方程的求解。 ( 3 ) 利用接触面上和上述等式约束所对应的动力学或运动学不等式约束条件 作为校核条件对解的结果进行检查。如果物体表面( 包括原假设中尚未进入接 触的部分) 的每一点都不违反校核条件，则完成本步的求解并转入下一增量步 的计算；否则回到步骤1 再次进行搜索和迭代求解，直至每一点的解都满足校 核条件。然后在转入下一增量步的求解。 2 2 2 接触界面的定解条件和校核条件 为了应用增量方法求解，将接触界面条件表达成增量形式。现假设物体彳 和占在t 时刻的解已经求得，需求解f + at 时刻的解，则 1 不可贯入性条件 指物体彳和曰的位形在运动过程中不许相互贯穿。 l 图2 3 接触点对及点对间的距离 如图，设t 时刻物体彳的接触面s 一上任一点p 至物体b 接触面s 占上最接 近点q 的距离可表示为 g ；( x a _ t x 曰) ，l 口 ( 2 2 1 ) 1 5 第2 章有限单元法 式中x 4 ，z 口分别为点p ，q 的坐标值，万占为接触面s 曰的单位法线向量。对 于f 时刻的不可贯入性条件应有 g 苫0 ( 2 2 2 ) 对于时刻f + at 有 + 出z _ 昌x _ + u _+ 址z 口墨t x 口+ “口( 2 2 3 ) 其中u 4 和u 且是t 至f + at 时间间隔的位移增量，则不可贯入性条件可写成增量 形式 舢g 。0 4 一“丑卜址厅丑+ ( z at z 口) 舢，z 曰；“a 一“另+ 丐芝o ( 2 “2 ) 其中“a ：u ao t + a t 以b ，“等；u bo t + t 甩占，i ；( x at z 口卜出，z 丑。 2 粘结时无相对滑动条件 云r 一比r a 一“岁= 0 当| f + 覃i d 分 离 此条件时无接触力作用的自由边通过搜寻检查上列条件，若不满足，则转 条件为粘结，并给出接触点对的位置 2 4 3 引入接触约束的方法 有了上述接触问题的定解条件，可以通过约束变分原理，直接约束法等方 法将其引入系统的求解方程，简单介绍如下。 1 拉格朗日乘子法 一 如果原问题可以表示为求解泛函。的驻值，且未知函数h 还必须满足附加 的接触约束关系， c 0 ) z 0( 在q 中)( 2 2 8 ) 这时引入这些附加条件构造另外一个泛函 口= 。+ j = ；矿c p q = i i + 口c 工 ( 2 2 9 ) 其中入是q 域中一组独立坐标的函数向量，称为拉格朗日乘子，口+ 称为修正泛 1 7 第2 章有限单元法 函，口r ，是用拉格朗日乘子法引入接触约束条件的附加泛函。这样原泛函口。的 有附加条件驻值问题就转化为修正泛函口+ 的无附加条件驻值问题。 这种方法是把约束条件加在一个系统中最完美的数学描述。该方法增加了 系统变量数目( 入) ，并使系统矩阵主对角线元素为零。这就需要在数值方案的 贯彻中处理非正定系统，数学上将发生困难，需实施额外的操作才能保证计算 精度，从而使计算费用增加。另外，由于拉格朗日乘子与质量无关，导致这种 由拉格朗日乘子描述的接触算法不能用于显式动力撞击问题分析。 2 罚函数法 构造修正泛函 口”- 月乙+ 口f 1 = 。+