（工程力学专业论文）大型水轮机转轮焊接残余应力热弹塑性有限元分析.pdf

西南交通大学硕士研究生学位论文第1 页 摘要 转轮作为大型结构，叶片与上冠和下环交界处均为焊接连接，焊接过程中 产生的残余应力是导致转轮疲劳裂纹萌生的重要因素之一。本课题采用三维热 弹塑性有限元方法计算大型转轮的焊接残余应力，取叶片与上冠和下环的连接 体作为分析模型，从而完整考虑叶片与上冠和下环之间在焊接过程中的相互影 响。在热分析过程中，采用集中质量矩阵方法以减小热分析结果的波动；选取 合适的载荷步以避免热分析最初加载时出现的温度振荡现象；在高温时选取合 适的材料性能参数以适应收敛性要求，结合顺序耦合场热应力有限元分析方法， 求得转轮的焊接残余应力分布和残余变形规律。将计算得到的残余应力分布与 实测的残余应力数据进行比较，两者符合较好，而且高温回火有效的降低了焊 接残余应力的峰值。焊接残余应力计算为后续的转轮疲劳寿命分析奠定基础， 同时为转轮的设计阶段考虑焊接残余应力的影响提供依据。 关键词：残余应力，焊接，转轮，有限元 西南交通大学硕士研究生学位论文 第l l 页 a b s t r a c t r u 肌e ri sal a e g es t r u c t u 阳，t 1 1 ec o n n e c t i o nb e t w e e nl a m i n aa n dc r o w ni sw e l d e d ，m e c o n n e c t i o nb e t w e e nl a n l i n aa t l db a n di sw e l d e d ，m ew e l d i n gr e s i d u a ls l r e s sa p p e a r e d i np r o g r e s so f 、v e l d i n gi so n eo ft h ei m p o n a l l tf k t o r sw h i c ha f r e c tt h er u n n e rf 乱i g u e 1 i f c t h i sp a p e rc a l c u i a t et h ew e l d i n gr e s i d u a ls t r e s so fl a r g er u n n e rb yt 1 1 r e e d i m e n s i o n t h e n a le l 臧i cp l a s t i cf i n i t ee l e m e m ，s e l e c t i n gt h ec o 皿e c t i o nw i t c hi n v o l v cl 锄i n a a n dc r o w na n db a r i da st h ea n a l y s i sm o d e l ，s of u l l yc o n s i d e r i n gt h ei n n u e n e e b e “v e e nl a m i n aa n dc r o 、v r la n db a n d i nt h ep r o g r e s so f t h e 皿a 1a n a l y s i s ，u s el u m p e dm a s s e l e m e n tm a i r i xt or e d u c et h eo s c i l l a t i o no f t h e r m a lr e s u l t ；s e t u pt l l ea p p m p r i a t el o a ds t e pt oa v o i d l h ep h e n o m e n ao ft e m p e r a t u r e s u 唱e a t o r i g j n a l i y l o a ds t e p ；c h o o s ea p p r o p r i a t em a t e a 1 p a 阳m c l e rt 0 e n s u r et h ec o n v e r g e n c e ，w i t ht h eo r d e rc o u p n gm e r m a ls t r e s sf i n i t ee l e m e m ， c a l c u l a t et h ew e l d i n gr e s i d u a ls t r e s sa n dr e s i d u a l d e f o r m a t i o n c o m p a r e dt h ew e i d i n g r e s i d u a ls t r e s sw i t ht c s t i n gd a 乜，也eb o t l la c c o r de a c ho m e r a n dh i g ht e m p e r a t u r e t e m p e rr e d u c et h ep e a kv a l u eo fw c l d i n gr e s i d u a ls 他s se 艉c t i v e i y1 1 1 ec a l c u l a t i n g o f w e l d i n gr e s i d u a ls t r e s se s t a b l i s ht h eb a s eo ff a t i g u el i f ea n a l y s i s ，o 髓r st h eg i s tf o r c o n s i d e 啦t h ea 腩c to f w c l d i n gr e s i d l l a ls t r e s si nt h ed e s i g no f m n n e r k e y w o r d s ：r e s i d u a ls t r e s s ；w e l d i n g ；r 吼n e r 二f i n i t ce i e m e n t 西南交通大学硕士研究生学位论文第l 页 1 1 问题的提出 第l 章绪论 近年来，随着我国水力资源的开发和利用，大型水电机组不断涌现，大直 径的水轮机转轮被广泛应用。但在国内已建成的大型电站中，比如岩滩电站， 李家峡电站，二滩电站，五强溪电站，小浪底电站，这些直径超过6 米的大直 径转轮，在投运后都不同程度的出现了疲劳裂纹，转轮是水轮机的关键部分， 裂纹的出现极大她影响了发电机组设备运行的安全性和可靠性。 现有的事故分析报告和文献资料大部分指出川，这些裂纹是由于转轮所受 的工作应力、交变应力和转轮焊接过程产生的残余应力共同作用的结果，残余 应力是引发转轮裂纹的重要原因之一，由于残余应力是焊接工艺过程中产生的， 无法彻底消除，而转轮残余应力的实测由于技术繁杂，且带有一定的破坏性， 不可能在转轮上大范围施行，因此十分需要从数值分析的角度弄清转轮中焊接 残余应力的分布规律和应力值大小，从而为转轮的抗疲劳裂纹设计奠定坚实的 基础。 1 2 焊接数值模拟技术的发展 焊接是一个牵涉到电弧物理、传热、冶金和力学的复杂过程。焊接现象包 括焊接时的电磁、传热过程、金属的熔化和凝固、冷却时的相变、焊接应力与 变形等。由于焊接过程的复杂性及影响因素很岁2 1 ，要想准确获得焊接应力和 变形的分布规律并加以调节和控制是非常困难的事，而有限元技术的发展和计 算机性能的提高为解决焊接应力和变形提供了可行途径。 2 0 世纪7 0 年代，日本大阪大学的上田幸雄等人首先以有限元方法为基础， 西南交通大学硕士研究生学位论文第2 页 提出了考虑材料力学性能与温度有关的焊接热弹塑性分析理论，导出了分析焊 接应力应变过程的表达式，从而使复杂的焊接应力过程的分析成为可能【3 】；r i k a r l s s o n 和b lj o s e f s o n 运用a d i n a 软件的三维有限元方法分析了圆管对焊时 残余应力、残余应变和变形，分析过程中舍弃了轴对称假设，材料为s i s 2 1 7 2 钢，焊接为坡口对接方式。结果发现较大的载荷步和较粗糙的单元划分仍然可 以得到满意的结果，在焊缝上存在很大的残余压应力和较小的拉应力，和考虑 轴对称情况下的结果比较后，总结了两种情况下偏差的意义，为以后只分析轴 对称模型奠定基础，减少计算时间鸭p t e k r i w a l 和j ，m a z u r n d e r 以a b a q u s 软 件分析了三维结构焊接的热传导问题，分析过程中，网格的划分随着焊接过程 的变化而逐渐细化，以适应焊缝区域的热传导，模型中考虑了热参数随温度的 变化、潜热影响、热传导和辐射边界条件，将热影响区和熔池区的温度变化结 果与试验结果进行了比较。结果符合很好巴s b ，b r o r i l 和h s o n g 认为：以往 的焊接残余应力分析一般集中在焊缝上，研究焊接过程中焊缝上的材料微观结 构变化、温度变化、焊缝扭曲等，数值模拟只取焊缝附近的一小片区域，是不 全面的，实际上，周围的广大弹性区域对焊缝附近也是有影响的，如热膨胀、 整体结构的固定方式等，工程上的焊接也需要考虑整体结构在焊接时的扭曲状 态，因此s b b r o w n 和h s o n g 以薄壁圆筒和外接圆环的焊接为例，分别以二 维和三维情况下的有限元模拟，比较了两种情况下是否考虑广大弹性区域对结 果的影响，结果表明二维和三维模型中，考虑广大弹性区与不考虑弹性区相比， 结果有很大的差别，相比而言，二维情况下，影响更大一些【6 】； s o h f i n g 和 h j l u c t 用三维瞬态数值方法模拟焊接熔池在焊接过程中的变化以估计熔池的 物理性质，焊接材料用流体金属喷射在熔池表面的方法模拟焊接材料施加在焊 西南交通大学硕士研究生学位论文第3 页 缝上，同时伴随着蒸发和潜热吸收，以3 0 4 不锈钢为例验证了数值分析方法p j 。 国内在2 0 世纪8 0 年代开始了关于焊接的热弹塑性理论及在数值分析方面 的研究工作，上海交通大学出版了一本数值分析在焊接中的应用【8 】；9 0 年 代以来，汪建华等对焊接热应力应交进行了分析，鉴于一般通用有限元软件难 以解决焊接问题，汪建华等与日本大阪大学上田幸雄等人合作开发了用于高温 热加工过程的三维热弹塑性有限元软件，采用热弹塑性有限元分析了压缩机的 焊接变形和管板接头焊接变形1 9 1 ，对压力容器多层焊问题进行了研究1 1 0 】；西安 交通大学的张建勋采用热弹塑性有限元法，分析计算了c o 合金静叶片电子束焊 接时的焊接工艺对焊接残余应力的影响，分析了焊接方向和焊件的装夹方式对 接头中的残余应力的影响】；西安石油学院的李霄和李栋才采用热弹塑性有限 元方法对超载拉伸消除焊接残余应力过程进行了数值模拟，对三种拉伸载荷下 的宽板和窄板进行了比较，发现拉伸载荷越大，残余应力消除得越好，但同时 也要考虑材料的塑性能力、构件允许的变形量和构件允许的残余应力的影响1 2 】； 蔡志鹏等忽略弹、塑性变形能和相变引起的能量变化，推导出焊接温度场相似 应满足的条件，并利用平板堆焊的结果进行了验证【1 3 l ；赵海燕等研究了焊接移 动热源、动态可逆的自适应网格技术、并行计算和材料性能在高温时的处理等 影响热弹塑性计算的因素，并以实际结构为例进行了验证，表明合理的处理方 法能有效的改善分析结构1 1 4 】；陈芙蓉等对b t 2 0 钛合金薄板电子束焊接后焊态、 焊后炉内整体热处理和电子束局部热处理三种焊接接头不同区域的力学性能和 断裂韧性进行了探索性试验研究，得出结果，电子束局部焊接可以在一定程度 上改善焊接接头的断裂韧性【l5 】；杨建等分析了陶瓷和金属接头的残余应力，得 出不同形状接头对残余应力的影响，介绍了以残余应变能来判断接头抗剪强度 西南交通大学硕士研究生学位论文第4 页 的思路。引。 目前已有很多使用热弹塑性有限元计算焊接残余应力的成功实例，但这些 例子都属于比较简单的问题，如平板对焊、可以近似为二维情形的薄壁环焊结 构、和其它一些可以将焊缝视为二维结构的模型。关于复杂三维模型下的焊缝 残余应力模拟，文献【1 7 】认为：随着数值模拟的发展和计算机容量和速度的提高， 三维下的模拟逐渐成为可能。以热弹塑性有限元方法计算三位复杂结构焊接残 余应力的最大困难在于：如果希望精确得出焊缝上的残余应力分布，需要将焊 缝划分为很密的网格和很小的载荷步，导致计算需要的时间是个天文数字。一 般的处理方法是选取一段焊缝及其附近部分区域作为研究对象，但这种方法没 有考虑选取模型周围的结构对结果的影响。 1 3 本文的研究内容 大型水轮机转轮的焊接残余应力数值分析在目前尚属于空白给结构和工 艺设计人员带来了极大不便。现在只能靠已发生的裂纹事故判断结构的危险部 位。本文的目标是分析转轮的焊接残余应力的分布，在完整考虑叶片、上冠和 下环之间相互影响的基础上，对涉及到的问题作一些简化，然后运用热弹塑性 有限元方法得出转轮的残余应力分布，和已有的实测残余应力数据进行验证， 确定转轮其它部位的残余应力分布，为后续的疲劳分析奠定基础。 西南交通大学硕士研究生学位论文第5 页 第2 章焊接残余应力涉及的问题 焊接过程的分析涉及到许多复杂问题，包括热影响区组织转变、析出相的 分布及其形态、硬化、软化、脆化和韧化的形成机制，焊接时氢的扩散、聚集 和溢出途径，弹塑性状态下的应力和应变，相变的影响，三维约束的影响，金 属结晶形态对组织和性能的影响，微量元素的变质作用等【1 8 1 。本课题主要使用 热弹塑性有限单元法计算热过程产生的残余应力和变形，不单独考虑相变等微 结构的变化对残余应力的影响。 目前在数值上对于复杂和大型焊接结构残余应力和残余变形的预测还存在 许多困难，主要表现在数值模拟中结构的自由度规模庞大、材料高温性能参数 缺乏、材料高温时的性能参数非线性变化影响了求解过程的精度和效率，大变 形和扭曲的影响在二维模型中属于可忽略因素，在三维模型中则有重要的影响 。焊接过程的复杂性导致在精确的数值模拟过程中计算量很大。以下将讨论 三维焊接过程数值模拟涉及到的问题。 2 1 材料高温性能参数的处理 材料的高温性能参数对焊接过程数值模拟的精度和计算效率均有较大影 响，目前材料高温性能参数方面的数据较少，弹塑性模型不能精确描述材料高 温时的本构关系，只能作一些经验估计；在熔化、半熔化状态下，材料的弹性 模量和塑性阶段的切变模量很小，导致有限元刚度矩阵的奇异性增加，解的收 敛性和稳定性变差，计算时间相应增加【”。 2 2 焊缝金属熔敷过程的描述 在进行焊接过程数值模拟时会经常出现多层焊和焊缝金属的熔敷现象，解 在进行焊接过程数值模拟时会经常出现多层焊和焊缝金属的熔敷现象，解 西南交通大学硕士研究生学位论文 第6 页 决这一问题的方法是采用单元生死技术【2 l 】。单元生死的本质是：在“死”单元 对应的刚度矩阵位置乘以一个很小的系数，默认值是l e 6 ，也可以修改此系数， 乘以此系数以后，“死”单元对应的刚度矩阵很小，相比于其它位置的刚度矩阵 可以忽略不计，相当于“死”单元不参与分析，以此模拟单元不存在的情况， 以后可以通过对该刚度矩阵乘以系数1 e 6 的方法重新使该单元活过来。在分析 过程中，这种技术可以通过使单元变“死”而有效的将这些单元从模型中消去， 在适当的时候激活这些单元则可以使这些单元重新参与分析。在焊接过程中， 随着热源的移动，填充金属形成熔池，这个过程的实现方法是，在划分网格时 将焊缝划成独立单元，并且在焊接前定义为“死”单元。模拟焊接过程时，焊 接处的单元按顺序被激活，与焊缝的填充过程同步，从而更真实的模拟焊接过 程。 2 3 分段移动热源 焊缝热源在时问和空间上分布是不均匀的，目前，比较好的热源模型为高 斯热源模型【19 1 ，假设在热源中间温度最高，向四周逐渐降低，成高斯分布，本 课题假设热源温度均匀分布。为了模拟焊接热源，需要采用很细的网格，时间 上也分成很细的时间步。如此一来，运算的规模与计算时间将会急剧增加，使 问题的求解难以实现。一般来说，对于一条焊缝，考虑焊接热源的移动速度， 可以将焊缝上施加的热源近似变换为等效的、提前一定时间作用的带状热源。 张建强等人比较了分段热源计算的焊接残余变形和残余应力与移动热源的变形 与应力，结果相差不大口2 1 。 西南交通大学硕士研究生学位论文第7 页 第3 章焊接过程残余应力数值模拟 在热弹塑性有限元分析时有如下假定：材料各向同性；材料的屈服服从 m i s e s 屈服准则；塑性区内的行为服从流变法则并显示出应变硬化；弹性应变、 塑性应变、与温度应变是可分的；与温度有关的机械性能、应力应变在微小的 时间增量内线性变化。 3 1 热场分析方法 热传导问题的控制方程如下【2 i 】： 胛2 口+ 孑= o ( 3 一1 ) i 彳 对稳态温度场 虿2 荦一e 票 对瞬态常温度场 。之 【譬一。百 嗣辩态吊温发场 k 为热传导率，口为温度，虿为物体内热源强度，c 为系数，等于比热和密度的 3 1 1 边界条件和初始条件 为了得到控制方程的解需要边界条件和初始条件。初始条件是物体在初 始时刻的温度分布情况，即： 6 1 = 岛g ，_ y ，z ) ( 3 3 ) 边界条件分为三类： ( 1 ) 如果物体部分边界0 上的任意点p 处，给定各时刻的温度为歹，则： 一0 ，) = 万0 ，f )p k ( 3 4 ) 西南交通大学硕士研究生学位论文第8 页 ( 2 ) 如果在物体部分边界l 上的任意点p 处其法向热流强厦为虿，则： g o ，f ) ：t 掣：虿o ，) p l ( 3 5 ) 咖 ( 3 ) 如果在物体的部分边界l 上的任意点p 处，与外界进行热量对流，或在边 界r 上接受辐射源的宿舍时，则： q ：七堂譬型： 晚o ，r ) 一只( p ，r ) 】 p l ( 3 6 ) 硼 和 g ：七掣：r e o ，f ) 一统o ，f ) 】p r ， ( 3 7 ) 册 h 为对流系数，k 为辐射系数，见、啡为外界温度。以为物体温度。 求解物体的温度场就是寻找满足控制方程和定解条件的温度函数。现给出 二类和三类边界条件的统一形式： 七掣：爵 ( 3 - 8 ) 3 1 _ 2 热传导问题的变分原理 取泛函为： 兀= 三f i 七 ( 署 2 + 愕) 2 + ( 署 斗矿一j i 口缸矿一护爵打 c ，圳 对于瞬态温度分析，孑中包含詈，其泛函变分问题目前尚未很好的解决。 因此，在下面证明时，是将时间变量t 暂时固定，把f 警看作和内部热源相类 似的量。待变分后，再考虑t 的变化。 西南交通大学硕士研究生学位论文第9 页 羽= i 七 警占( 警 + 考攻考 + 警6 ( 髫 d 矿一j j f 甜缸矿一f 瓣盯 由于 f f 宴删矿 j j ，苏2 将( 3 - 1 l 3 1 3 ) 式带入( 3 ，l o ) 中等号右边的第一个积分式得： 根据高斯公式有： 将( 3 - 1 5 3 1 7 ) 代入( 3 1 4 ) 式 册 ( 3 1 0 ) ( 3 1 1 ) ( 3 1 2 ) ( 3 1 3 ) ( 3 ，1 4 ) ( 3 1 5 ) ( 3 1 6 ) ( 3 1 7 ) 矿 y 剃 涮 目丁 口 矿一砂扩一如 m m 一 一 一 矿 矿 y p p p 船 册 阳 加一缸卯一砂卯一瑟 a一苏a一砂a一出 m m m = = = 矿 矿 矿 p p p 卯一缸卯一砂阳一良 哎 阳一缸卯一砂卯一如 m 小m 门川 翌出 塑锄 + 、 塑砂 r 科l 塑砂 + 、， 塑撕 哎 塑缸 r-j叫l m 塑出 ，l a 一瑟 ， 、，订甸蚴 挚牛馐纠 文盈瑟 斗 ，一： 卜冶一旷叫鱼砂 翌舐殳一 几u鱼缸 a一缸，一 m 1 f j 锄一苏翌锄 和 一一 锄瓦翌苏 和 = 甜 望知a 一缸 肚 陪 塑砌 归 锄一砂塑砂 归 叫 卯一砂峙 跏石塑锄 和 = 锄云塑出 妒 =甜 塑瑟畦 西南交通大学硕士研究生学位论文 第l o 页 一= 弹七署2+ ( 弘( 害 2j f 4 窘+ 窘+ 窘卜d 矿 = 弹t 警跏一m 4 害+ 雾+ 窘p d 矿 将( 3 1 8 ) 式代入( 3 1 0 ) 式有 铘= 一j 阿窘+ 寄+ 害 + 孑 删矿一弹( 一t 警+ 一r ( 3 1 8 ) ( 3 一1 9 ) 由6 n = o 和变分笳的任意性，可得到( 3 1 ) 和( 3 - 8 ) 式。可见，温度函 数如果满足式( 3 3 ) 和式( 3 4 ) ，而且使使( 3 1 0 ) 所示的泛函一阶变分为零 则此温度函数便是要求的真实解。 设： 边界面为 匿考钾l 苏 砂出j r = l + l + l + ( 3 2 0 ) ( 3 2 1 ) 温度函数在l 上满足第一类边界条件，如将式( 3 5 ) 、( 3 6 ) 和( 3 2 ) 代 入( 3 9 ) 中，并令瓤= o ，得 m 舱 。) 7 。p 矿= j f _ 一刚嘲盯+ f ，r n 一刚嘞r + m p c 詈 删y ( 3 2 2 ) 式( 3 2 2 ) 中等号右边为热生产率，左边为热传到率，它们是相等的。这样一 来，求解物体的温度场就是寻找满足热流平衡方程( 3 2 2 ) 和约束条件( 3 3 ) 、 泖 、， 西南交通大学硕士研究生学位论文 第l l 页 ( 3 4 ) 的温度函数o 3 1 3 热传导问题的有限元公式 用m 个单元将物体进行离散。对于第m 个单元，节点数为n ，温度插值时 的形函数矩阵为防】，且相应的形函数为石，五，无， 眵】- 匦，五，五 ( 3 2 3 ) 则单元内任一点的坐标值扛，y ，z y 为： x = 弘妊) j ，= 眵抄) ( 3 2 4 ) z = 瞄mj 其中扛_ y z ) 为单元节点的坐标列阵，温度插值公式为： 目= 眵) = 舻 7 眵】 ( 3 2 5 ) 设 阱 旦茜昙 2 b z s ， 将式( 3 2 5 ) 代入( 3 2 0 ) 有： 。 = p 痧) ( 3 2 7 ) 设 【岛】= p 痧】 ( 3 2 8 ) 将式( 3 2 8 ) 代入式( 3 2 7 ) 得： 。) = 阮】p ( 3 2 9 ) 由式( 3 2 5 ) 得： 西南交通大学硕士研究生学位论文 第1 2 页 署= 嚼p ) = 根据式( 3 ，2 2 ) ，对于第m 个单元有： m 。，尼占 。y 。扣矿= 几。，虿潮+ f f ，。，a 他一见p 簖 + n k 一只胁+ j 卜c 箬 将式( 3 。2 5 ) ，( 3 2 7 ) ，( 3 3 0 ) 代入式( 3 3 1 ) 得 ( 3 3 0 ) ( 3 3 1 ) 占p y f i 。【岛r “岛p 矿p ) = j p y 飓。眵r 虿打+ 万矽 r 几。事j r | b 眵弦( 姚一p u ，。 + 万p ) 7 肌，眵i r f 眵 ( 姚一+ j 坩f j f 计虿d y 一占p 7 f 。眵r c 眵p 矿p 设单元传导矩阵忙。1 为： 单元对流矩阵【七。】为 单元热容矩阵k 】为： 节点热流列阵为： k 】= f 。，阮n 阪矽 陋。】= i 咖，计a 眵k k 。】= j f i 计c 鼬矿 k ) = 儿。，计虿打 ，2 ) = f f 眵r 虿d y 将式( 3 3 3 ) ( 3 3 7 ) 代入( 3 3 2 ) 式得 ( 3 3 3 ) ( 3 3 4 ) ( 3 3 5 ) ( 3 3 6 ) ( 3 3 7 ) y剃 、1， 西南交通大学硕士研究生学位论文 第1 3 页 【k 】p ) ( ：“) + k ) 一陋。】p y + k 。如) ：一p 妒，) + 【七，如) ，一p ，) ( 3 3 8 ) 式( 3 3 8 ) 为第m 个单元的热流平衡方程。 如果将式( 3 3 8 ) 对所有单元集合，则可得整体热流平衡方程 【颤】舻 ：翻。) + 伍： 一【x 。】矽 + k 。】( ，一和1 ) + k 胎) ，一p l ) ( 3 3 9 ) 式( 3 3 9 ) 中，p ) 、p 分别为由所有单元节点的温度和温度对时间导数组成 的列阵：忸。) 为因边界热流输入而形成的整体节点热流列阵； 月：) 为内部热源 而形成的整体节点热流列阵，其它矩阵与m 单元的意义相似。在不同的时刻求 解式( 3 3 9 ) 可得不同时间的节点温度值。 3 1 4 线性热传导问题的有限元解法 假设没有辐射，且热传导率k ，对流系数h ，系数c 均与温度无关，则问题 简化为线性热传导问题。 因为是瞬态问题，所以在不同时刻，温度值是不同的，求解时需将时间域 进行离散，求各离散点处的温度值。根据对的处理形式不同，将求解方法分 为前差分法、后差法法、平均公式法。以前差分法为例。 设时刻t 的温度p 。已知，温度对时间的导数在t 时刻的数值为： 址掣 b 。 将式( 3 4 0 ) 代入( 3 3 9 ) 得到t 时刻的热流平衡方程： 古k 。】p l + 。；忸1 + j r ：l + 古k l + k 】 p ) 。一p l ，) 一k p l ( 3 - 4 1 ) 设 怫幽 o1 r ( 3 4 2 ) 西南交通大学硕士研究生学位论文第1 4 页 忙 = + 姒+ 掣舻l + k 。x 乩一一【醚m ( 3 _ 4 3 ) 陋j 、陋】称为有效热传导矩阵和有效热流列阵。式( 3 4 1 ) 可简化为： k 如h = 伍l ( 3 州) 式( 3 4 4 ) 为线性代数方程组，求出k 】和t 时刻的忙) ，便可求出f + f 时刻的 温度场，如此不断重复下去，可得到任意时刻的温度场分布。 后差分法的p 表示为： 其它的与前差分法思想一致。 3 2 应力场分析 ( 3 4 5 ) ( 3 4 6 ) 3 2 1 材料力学性质依赖于温度变化时的增量型本构关系 ( 1 ) 弹性状态时，由广义虎克定律可知： 叠l = 【d 】_ 1 p ) ( 3 4 7 ) 【d l 为弹性矩阵，但是材料常数时温度日的函数。对式( 3 4 7 ) 求导可得： 雌i = 簪枷耐+ 时j ( 3 - 。7 ) 设温度引起的变形为和b ： p ) 。= k 。占一。，。k 。7j ( 3 4 8 ) l 一 留一 一 一一 二， 竺王出 学 掣 她 扣一 西南交通大学硕士研究生学位论文 第l5 页 忙 。= 缸p ( 3 4 9 ) 其中： 缸) = 口【1 llo o o l ( 3 5 0 ) 在此只考虑线性热传导问题，线膨胀系数口与温度无关，对式( 3 4 9 ) 求导得： d 扛k = 向弦臼 ( 3 5 1 ) 全部应变扛) 的微分为： d p ) = d 扛l + d 忙k ( 3 5 2 ) 将式( 3 4 7 ) 和式( 3 5 1 ) 代入式( 3 5 2 ) 可得： 州= ( 譬缸归州制 s ， 变换式( 3 5 3 ) 可得： 州埘p ( 等扛州 s 4 ， 设： 删。- 簪矗炉9 协s s ， 则有： d p ) = 【d 1 0 叠 一d 譬k ) ( 3 5 6 ) ( 2 ) 塑性状态时，全部应变包括弹性分量、塑性分量和由温度引起的应变分量三 部分。 d 扛 = d 叠) 。+ d 叠 ，+ 谁) 。 ( 3 5 7 ) 西南交通大学硕士研究生学位论文 第1 6 页 其中塑性应变分量推导如下。根据普朗特一柔埃斯塑性应力应变关系得： 喵= 等s ， ( 3 5 8 ) 将等效应力彳展开有： 孑= 击肛i 闩i 玎石i 而再虿砑 s 。， 将式( 3 6 0 ) 对各应力分量求导，可得 a 万 3 s ， 0 盯2 孑 将式( 3 6 1 ) 的矩阵形式代入式( 3 5 8 ) 可得： d p ) 一2 磊晖 将式( 3 4 7 ) 、( 3 5 1 ) 和( 3 6 2 ) 代入式( 3 5 7 ) 有； d 阱簪p 枷+ m p 嗝肌矗 式( 3 6 3 ) 经过变换并利用式( 3 5 5 ) 可得： 州城( 荆一为拒l 俐0 因为硬化曲线与温度有关，可设： 厅= h ( 肛，臼) ( 3 6 1 ) ( 3 6 2 ) ( 3 6 3 ) ( 3 6 4 ) ( 3 6 5 ) 肌( 磊卜阳镛一十等d 口 s s ， h b 3 一 一一 里 将加删两边乘以( 为) 7 酬川，确，代入得 一圃坐! ：域竺! 塑。， + l 磊j 【d l 为 d 云一可用来判断单元是否处于加载状态。 其中 将式( 3 6 7 ) 代入式( 3 6 4 ) ，整理得 d 扛 。鼎蒯，燃i d 毙) + d 爷k 6 8 ) = 【d l p 忙 一d 舻 。) + d 忙b ” 俐一2 辗 ( 3 6 9 ) ( 3 7 0 ) 式( 3 - 6 8 ) 为塑性状态时增量型本构关系。其中【d l 、d 爷 。均可以同样的方式 展开，现描述如下 在三维问题中 a 于 3 s ， 羽2 l i 3 1 7 i j ( 3 7 1 ) ，旦堕纠靛 叫i 西南交通大学硕士研究生学位论文第1 8 页 在此，有s ，+ s 。+ s ：= 0 ，所以： 【。l 磊= 詈b ， ( 3 7 2 ) 隅m = ( p l 磊) 7 = 等k 】 ， ( 为m 为= 器雠+ $ 亿托如z r 三) 瑚b ，a ， 将式( 3 7 2 ) 、( 3 7 3 ) 和( 3 7 4 ) 代入式( 3 6 9 ) 整理可得 吼= z s ：s 。sys l sz s ，ws 。f 口s ，。 s j s y s ：sy f w s y f 口s ，。 s ： s ，vs ，。sz i 。 f 刍f 删r ff 删f 玎 对称 f 三乇 3 1 2 2 应力场分析的变分原理 构造泛函： ( 3 7 5 ) 肚吉埋砟y d 和沙：j d 扛y 阮d p ) ( 3 - 7 6 ) + f f d p 7 d 伶l j 矿一几d c 厂扫沁 一 西南交通大学硕士研究生学位论文 第1 9 页 丌= 主l 。肛) r 时阮怫d 矿 一肌妇) r 时荆。d 矿 ， + i k 妇耐”垆旧7 酬爿。d 肛儿神p 田“垆纠7 矾p 汪 由变分原理可知，节点真实位移增量应使泛函的一阶变分为零，即： 薹! 陋r f d l 陋p y ) ：_ r 一薹c 。，陋r 【d l 秽k d y ( ，9 ) + 薹i 。州陋】7 痧) 。d 矿一薹儿。， 扫沁= 。 善m 。，陋r f d b 陋p 矿= k 】 ( 3 8 0 ) 薹m 。，陋【。k p k d 矿一萋j f f 。，陋i r 侈b d 矿十羹地。渺i r 扫洒= 球 式( 3 - 8 2 ) 在形式上与一般有限元方程一样，不同在于载荷项忸 中多了温度 西南交通大学硕士研究生学位论文 第2 0 页 在弹塑性有限元中，非线性方程的求解需要将方程线性化，把非线性问题 转化为线性问题大多数采用逐次加载法。如果每一个加载步所加载荷足够小， 它引起的应力微分和全应变微分，可用相应的增量值扫 和扛) 代替，则扛 和料之间成线性关系。 假设现为第i 步求解，【d k 是此次加载之前应力值p ) 。的函数，k 】也应 为b l 的函数因此将式( 3 8 2 ) 改写为： k p l 一。】p ，= 忸) ( 3 _ 8 3 ) 式( 3 - 8 3 ) 的计算过程可如下进行：设物体开始进入屈服时的节点位移，应变 和应力分别为p ) o 、忙 0 和p ) 。，在此基础上再加载时，应用p ) 0 计算【d 】。和 k r ，并集合为k 】，再代入( 3 8 3 ) 式求得p l ，由p ) l 解出扛 和p ) ， 此次加载完后的应力值为： p ) 。= p ) 。+ p 1 ( 3 - 8 4 ) 扫 。用来计算下一次加载时的【d 】。和k 】”的计算。如此重复计算，直到全部 载荷求解完毕。 此方法再计算每次加载时的节点位移增量时都必须重新计算刚度矩阵， 称之为变刚度法或切线刚度法。因为每次重新计算刚度矩阵需要大量的时间， 一种简化的方法为：可以每5 次计算一次刚度矩阵2 3 1 。 3 3 。2 加权系数法 在结构中，一部分单元是在此次加载完了仍不屈服或开始加载前已屈服， 还有一部分单元是在此次加载过程中屈服，后者称为过渡单元。对于弹性区， 西南交通大学硕士研究生学位论文 第2 1 页 可直接用弹性矩阵【d l 计算刚度矩阵：对于过渡单元，无论用【d 】。，还是【口k 计 算刚度矩阵，都会带来误差。因此，用如下定义的加权平均弹塑性矩阵陋k 来 计算刚度矩阵： 陋l ：p l d l + ( 1 一p ) 【d 1 ， ( 3 8 5 ) 过渡单元的刚度矩阵为： 医 = m 。【b 】r 陋k 陋p 矿 ( 3 _ 8 6 ) 式( 3 8 5 ) 中的p 为加权系数，是单元达到屈服所需要的等效应变增量与此次 加载在单元内所引起的等效应变增量的比值。对于过渡单元，o p l 。当p = l 时，表示此次加载后不屈服的单元，当j d = o ；表示此次加载前便已屈服的单元。 因此，式( 3 ，8 5 ) 对于过渡单元与非过渡单元均适用，也可以用来计算弹性单 元与己屈服单元的刚度矩阵。 主要计算步骤如下： ( 1 ) 将全部载荷趣 力在离散后的物体中，并按线弹性计算，得节点位移p ) ， 应变p ) ，应力p 。 ( 2 ) 计算出每个单元的等效应力云，取最大值厅_ m 。如果氏。盯，则材 料尚未屈服，则计算结果为最终结果。否则，设：垒l ，由于在弹性时，各 仃1 物理量均与载荷成线性变化，故加载过程中，物理量按同一比例变化，物理开 始屈服对的节点位移p k 、应变黏 。、应力扫) 0 以及所加的载荷 足j 0 应为： 西南交通大学硕士研究生学位论文 第2 2 页 p o ：掣 l 乩：掣 l b o = 譬 尺 0 ：警 ( 3 ) 物体进入屈服后应加的载荷为仁) 一掣。如果按n 次均匀加完，每次加 的载荷为： 州= 粤( ，一圭 订、l ， ( 4 ) 施加忸) 。 ( 5 ) 由忙 。计算出等效应变，并计算出单元的p 值。 ( 6 ) 由式( 3 - 8 6 ) 计算单元刚度矩阵医 。 ( 7 ) 将单元刚度矩阵整合为整体刚度矩阵。 ( 8 ) 解线性方程组式( 3 8 3 ) ，得到节点位移增量。计算各单元的应变和等 效应变，对各单元计算为进入屈服需要的等效应变增量，修改p 值。重复( 5 ) ( 8 ) 三次左右，可得到比较理想的p 值。 ( 9 ) 计算应力增量。将得到的位移增量，应变增量和应力增量叠加到加载前 的相应值上，得到此次加载后的位移，应变和应力。 0 国重复( 4 ) ( 9 ) 步，直至载荷施加完毕，结束。 3 4 有限元软件的应用 3 4 1 热弹塑性分析方法 西南交通大学硕士研究生学位论文第2 3 页 热弹塑性分析采用a n s y s 有限元通用分析软件完成，其基本分析过程如 下。 当一个结构加热或冷却时，会发生膨胀或收缩。如果超静定结构各部分之 间膨胀收缩程度不同，或结构的膨胀、收缩受到限制，就会产生热应力。热应 力的分析是一个耦合问题。焊接过程的有限元模拟是一个高度非线性的求解过 程。a n s y s 提供了三种进行热应力分析的方法【1 8 】： ( 1 ) 在应力分析中直接定义节点的温度。如果节点温度已知，则可以通 过b f 念令直接定义节点温度，各个载荷步节点温度作为载荷施加在模型上，以 此来求得结果。相当于第二种顺序耦合方法省略了热分析的过程。 ( 2 ) 顺序耦合方法( 问接耦合方法) ：首先进行温度场的分析，得到各个 载荷步时温度场的分布。将第一步温度分析的结果作为载荷施加在应力分析模 型上，得到第一步应力分析的结果，将其它后续温度分析载荷步的结果变化施 加在应力分析模型上得到应力场变化后的结果，如此持续下去，可得到最终的 应力场的分布。它通过把温度分布场的分析结果作为应力场分析的载荷来实现 两种场的耦合，但应力场的结果不影响热场分析，属于单向耦合。 ( 3 ) 直接耦合方法：利用包含温度和位移自由度的耦合单元求解，即可 同时得到温度场和应力场的分析结果。这种情况下，耦合通过计算包含所有必 须自由度的单元矩阵或单元载荷向量实现。 如果节点温度是己知。适合用第一种方法但节点温度一般情况下是不知 道的。对于大多数问题，使用第二种方法。顺序法可以使用所有热分析和应力 分析的功能。如果热与应力的耦合是双向的，即热分析影响应力分析，应力分 析也影响热分析，可使用直接耦合方法。直接耦合法在解决耦合场相互场作用 西南交通大学硕士研究生学位论文 第2 4 页 具有高度非线性时更具优势，并且可以一次性计算得到更好的结果。求解这类 问题均有专门的单元供直接利用，但在a n s y s 的热应力问题中，三维耦合单 元不具有塑性功能，只能进行弹性范围内的计算，二维情况下，使用耦合单元 比较方便。应力场对温度场的影响不严重的情况下，顺序耦合解法更为有效和 方便，可以独立进行两种场的分析。严格来说，对焊接问题来讲，应力场和温 度场应该是相互耦合的，不过应力场对温度场的耦合效果很小，而且焊缝附近 的温度场变化很大，材料的各种物理性能也变化很大，这种影响与上述耦合效 应相比要大得多，所以就焊接的热弹塑性问题说，可以采用非双向耦合场的分 析，即温度场和应力场单独进行计算的间接法，首先进行焊接过程温度场的分 析，然后根据温度场分析所得结果进行热弹塑性分析。 3 4 2 分析步骤 ( 1 ) 进行热分析。可以使用热分析的所有功能，包括传导、对流、辐射 等进行稳态或瞬态热分析。注意划分单元要充分考虑结构分析的要求，比如在 可能有应力集中的地方，网格要密一些。如果进行瞬态分析，需要注意热梯度 最大时的时间点或载荷步，此功能由a n s y s 自动完成。 ( 2 ) 新建一个文件，读入刚才的热分析模型，进入前处理，将热单元转 换为对应的结构分析单元。热单元的轴对称不能直接转到结构单元中，需要手 工设置。 ( 3 ) 设置结构分析，包括温度的材料属性及前处理细节，如节点耦合和 边界条件。读入热分析的节点温度，将节点温度作为载荷施加在模型上。 ( 4 ) 进行求解和后处理。 需要注意的是a n s y s 本身提供了一个用于三维热结构耦合计算的耦合单 西南交通大学硕士研究生学位论文 第2 5 页 元，通过此单元，可用直接耦合法计算，很方便，但是此单元没有塑性功能， 即只有在处于材料弹性范围内时，计算才算有效。对于温度比较高的情况下， 可能导致塑性变形的情况下，不能使用这种方法。 3 5 提高计算精度和稳定性的若干途径 焊接时局部温度很高，材料性能随温度剧烈变化。由于熟源的移动，整个 焊件的温度、应力应变随时间和空间急剧变化，同一时刻可能存在加热和冷却、 加载和卸载现象。焊接问题包含着高温塑性、材料性能随温度剧烈变化以及大 变形等高度非线性因素，使得要得到好的计算精度和保证解的收敛性相当困难， 本次分析为改善三维焊接有限元分析精度采用了以下方法： ( 1 ) 材料性能在高温时的处理。在温度很高时，缺乏材料性能的实验数 据，此时假定材料的弹性模量和屈服应力均为很小数值。然而，如果假定不合 理，可能产生卸载时应力反而增大的不合理现象。为防止该现象。可认为由于 弹性模量随温度改变而导致的应力变化应小于热膨胀引起的应力变化2 4 1 ，即 竺盯 e 以丁 ( 3 9 0 ) e ( 2 ) 适当的计算步长和网格划分。在热分析加载的最初阶段，会出现结 果振荡，最大值和最小值均超过规定的范围，违背热传导规律2 5 1 。在加载的最 初阶段，传热需要时间，若步长太小，在热量尚未传到较冷区域时，该区域应 保持不低于下限的温度；若步长太大，在热源或介质相邻区域的温度会急剧上 升，该区域应保持不高于上限的温度。当待解的参数不受约束时，解的范围可 在实数全空间取值，为了使泛函最小，结果尽快收敛，程序让应保持下限温度 的区域降温，或者让应保持上限温度的区域升温，产生了违背热规律的解。因 西南交通大学硕士研究生学位论文 第2 6 页 此，对载荷步长需要进行限制。焊接开始时的温度步长r 不能太大，最好控制 在o 1 s 之内，也不能太小，此时不要打开自动时间步长功能，待温度下降到一 定范围时，可增大时间步长并打开载荷步自动搜索功能。网格划分应考虑焊接 接头处较密而远离焊缝处较疏。以满足温度和应力计算的精度。 ( 3 ) 在瞬态温度场的分析中，发现在开始起步的若干时间内，会出现温 度场结果违背热传导规律的情况，热从商温处向低温处传导时后，温度的变化 应该是单调的，但结果有时是波动的。 单元热容矩阵一般为协调质量的热容矩阵脚】，是指在单元内，竽与r 的插 讲 值函数是一样的，一般均采用此方法。采用协调质量热容矩阵会使节点热流平 衡方程失去明确的物理意义。单元热容矩阵通常是满阵，累加后的热容矩阵必 然有非对角线元素，假设热容矩阵的第m 行第n 列元素为七。o ，则总体热平 衡方程( 3 3 9 ) 中第m 行除了含譬外，还有婴，其第m 个方程应是m 节点 u h lo n 的热平衡方程，流入m 节点的热量与m 节点的温度升高保持平衡，而含有答， 册 说明m 节点的热平衡与n 节点的温度升高有关，这是不合理的。如果把单元质 量离散到节点上，即取集中质量的单元热容矩阵，则单元热容矩阵为对角阵， 累加后的热容矩阵也为对角阵，便不存在上述问题，a n s y s 中的命令为l u m p m 。 西南交通大学硕士研究生学位论文 第2 7 页 第四章转轮焊接残余应力分析 4 1 建立转轮计算模型 计算采用的模型是三峡水轮机右岸转轮，如图4 1 所示 图4 - l 转轮模型 4 1 1 单元划分和对称条件的使用 转轮由上冠