1000MW汽轮机滑参数停机优化研究之汽轮机转子温度场及应力场的有限元计算毕业论文.doc

1000MW汽轮机滑参数停机优化研究之汽轮机转子温度场及应力场的有限元计算毕业论文目录中文摘要IAbstractII目录III1.引言- 1 -1.1 课题的目的和意义- 1 -1.2 国内外研究发展状况- 2 -1.2.1 国外发展现状- 2 -1.2.2 国内发展现状- 3 -1.3 本文研究内容- 5 -2. 温度场、应力场模型的建立- 7 -2.1 概述- 7 -2.2 二维离散温度场模型- 7 -2.2.1 温度场模型的建立- 7 -2.2.2 各部分放热系数的计算- 10 -2.3 温度场的变量说明及运行程序- 11 -2.4 应力场数学模型的建立- 19 -2.5 汽轮机转子的合应力与强度指标- 21 -2.6 滑参数停机应力场分析- 22 -2.7 应力场变量说明及运行程序- 22 -2.8 本章小结- 32 -3 转子温度场与应力场的有限元计算- 33 -3.1 概述- 33 -3.2 转子有限元分析网格划分- 33 -3.3 几何模型的建立- 35 -3.4 有限元计算分析方案- 35 -3.5 转子温度场、应力场计算分析- 35 -3.6 本章小结- 40 -4 转子寿命分析- 41 -4.1 转子寿命损耗评价- 41 -4.1.1 转子寿命损耗- 41 -4.1.2 影响转子寿命的因素- 41 -4.1.2.1 材料因素- 41 -4.1.2.2 应力因素- 42 -4.1.2.3 结构尺寸因素- 42 -4.2 滑参数停机方案的比较- 42 -4.3 本章小结- 43 -5 结论与展望- 44 -5.1 本文主要工作- 44 -5.2 后续工作展望- 44 -致 谢- 45 -参考文献- 46 -III10011数停机优化研究1.引言1.1 课题的目的和意义能源是国家经济发展的重要支柱，电力在能源生产中起着重要的作用，电力生产是现代工业的的基础之一，其产量直接或间接的影响一个国家发展的现代化水平。很多地方由于受到资源环境制约，水利资源缺乏，形成以火电为主的电源结构。文献指出，在2003年，全国发电量18462亿千瓦，其中火电发电量14193亿千瓦，占全国发电总量的76.8%，所以工业和生活用电大部分依靠火力电站和核电站。绝大多数火力和核电站发电是以加热蒸汽推动安装于转子上的叶片旋转，给转子施以动力带动连接在转子上的发电机发电生产。随着现代电力资源日趋贫乏 ，电力生产中将逐步加大对新能源的开发利用，如太阳能发电。有资料一显示，国外在太阳能发电技术应用上，德国Schlaich教授1978年提出利用太阳能加热气流形成的动力驱动涡轮机带动发电机发电。另一方面汽轮机作为原动机械具有功率大、热经济性高、单位功率和制造成本低等一系列优点，所以在冶金、化工等部门广泛应用以作为驱动各种泵、风机等的动力机械。转子是汽轮机中是带动发电机旋转，间接把热能转化为机械能的部件，是汽轮机中的重要部件之一。大型汽轮机机组转子在运转过程中，工作条件恶劣月复杂，运行中不仅要承受自身、叶片的离心力、稳态汽流力以及交变应力的共同作用，还受到湿蒸汽的腐蚀作用以及高温给材料带来的蠕变 作用。近年来，用电量随着经济的飞速发展而迅速上升，能源产业作为经济发展支柱，也随之迅速发展。全国各地，相建立了大批大容量发电机组，为确保能源的合理利用，机组参与调峰运行也是必然。机组参与调峰运行，受启动频繁和大幅度负荷变动的影响，机组的主要部件经常承受大幅度的温差变化，由此引起汽缸、转子等部件承受交变热应力。启动或停机时部件承受热应力的大小，决定于部件内部的温度梯度。对于汽缸等厚壁部件，由于机组高压缸的设计普遍采用了双层缸结构，而且汽缸壁的金属厚度远较转子为薄，蒸汽对汽缸内壁的放热系数也远比转子小，因而启动时的径向温差及热应力远比转子小，汽缸热应力水平已显著下降。而转子的热应力却随着直径的增大而显著增加，况且汽转机转浙江大学硕士学位论文子长期在高温工作，受力情况复杂，转子上除了热应力外，还存在各种机械应力同时在高温高压工质中高速旋转，转子承受叶片、叶轮及转子自重产生的离心应力，以及蒸汽压力对转子产生的压应力。另外，还有转子传递扭矩引起的剪应力，自重引起的交变拉、压应力等。以一次启停的中心孔转子为例，在机组启动加负荷时，转子外表面的温度高于中心孔表面，外表面承受压应力，中心孔表面承受拉应力;当机组稳定运行阶段，转子外表面与中心孔的温度将逐渐趋于一致，热应力也将随之趋于零;而当机组停机降负荷时，转子温度场和应力场的分布将与启动时相反。这样，在机组一次启停过程中，转子外表面和中心孔表面都承受一次交变热应力，这种交变热应力将对转子产生疲劳损伤。使得金属表面特别是应力集中部位萌生裂纹并逐渐扩展以致断裂 。可见转子是工作条件最为恶劣的部件之一，在一定程度上汽轮机转子的寿命代表了整台汽轮机机组的寿命1。 因此对大容量机组转子热应力和寿命进行研究。合理控制机组的启停不仅关系到机组- 47 -的经济性，更涉及到机组的安全性。对转子进行热应力和寿命的有限元计算正是根据传热学，弹性力学和疲劳强度等基本原理，通过温度场和应力场的计算，采用局部应力计算法，结合金属材料具体的疲劳试验曲线，而求得转子的低周疲劳损伤，不仅对转子的寿命进行正确的预测，更对在线寿命监测系统可靠性做出正确的评价。本学位论文主要针对汽轮机高中压转子寿命问题中采用有限元方法进行分析根据其热应力等参数的分布特性开展研究。结论，有其深远的理论意义和极大的实际应用价值：(1) 通过对汽轮机高中压转子工况的数值模拟，可以进一步加深对汽轮机高中压转子运行中的特点认识。(2) 通过对汽轮机高中压转子运行中的特性认识，可以对进行汽轮机高中压转子热应力分析有一定的理论指导作用。(3) 通过对汽轮机高中压转子的数值模拟研究，可以清楚地了解其存在现有产生不足问题的主要原因，找到容易产生裂纹或者疲劳的部位，并可以节省大量的实验过程和实验时间及花费；对提高汽轮机高中压转子在实际应用中的综合性能，减少疲劳损伤以增加寿命有着非常重要的价值。1.2 国内外研究发展状况1.2.1 国外发展现状由于目前转子的温度和应力尚不能直接进行测量，只能通过间接方法，建立相应的数学模型，测量相关参数，求出转子金属温度和应力的变化及寿命损耗。现在转子应力的数学模型大多数是采用一维温度场理论解的简化式，其计算精度较低，只能反映应力的变化趋势，而不能得到应力的精确值。若在此基础上计算转子在启停和变负荷过程中的寿命损耗，将会产生较大的误差。国外机组寿命管理的应用在日本、美国和欧洲较为普遍。美国自60年代Gollin电站汽轮机失事以后，一些大的公司和研究机构如GE、WESTINGHOUSE、EPRI等对转子的安全性更为重视；进行了深入的研究。他们将有限元等先进数值方法由于汽轮机转子的分析计算，对转子材料的低周疲劳、高温蠕变、低温脆性和裂纹扩展规律等诸多方面的问题进行了大量的研究，并在汽轮发电机组上安装了应力及寿命损耗指示器以指导机组运行。日本在汽轮机寿命管理方面也做了很多工作，除了预测可能出现裂纹的寿命外，还对转子剩余寿命做出计算。日本的Kagawa University的Ebara等对汽轮机动叶片采用的12Cr钢和Ti-6AI-4V合金的疲劳特性进行研究，Fujiyama，Kazunari；Takaki，Keisuke；Nakatani，Yujiro等根据统计损伤和随机损伤仿真研究，对汽轮机设备进行寿命评估，采用先进技术设计汽轮机流通部分，以提高机组的性能和设备的可靠性。另外，日本在无损探伤的研究方面处于世界先进水平，日立、三菱重工、东芝、富士机电等著名大公司相继提出脆化一腐蚀法、硬度法、金属组织法、电极化法等无损探伤方法作为改进转子寿命评估的手段。德国的Wichtmann，Andreas研究了高温对汽轮机部件的蠕变损伤；Zaviska，Reichel研究了汽轮机冷态启动过程中的转子温度变化，在此基础上建立了冷态启动仿真模型； Scheefer，M；Knodler；Scarlin，B等对电厂抗高温、高压材料进行了探讨，一方面是发展新的材料，一方面是在已有的材料表面喷涂抗氧化性能强的图层；以及关于机组安全经济运行方面也进行了大量的研究.由美国萨金伦迪公司负责人总体设计，并由美国燃料公司和瑞士苏而寿公司提供锅炉，由瑞士ABB提供汽轮机发电机组的超临界600MW机组配有与自动启动序接口的热应力控制装置，它安装在高、中压转子进汽端的转子温度探针和一套计算回路组成。在启动升速和加负荷过程中，该装置都参与控制。运行部门反映该机组的最大优点是启动速度快，但也提出不能以应力作为控制的唯一依据，还应考虑到汽缸的膨胀是否均匀。英国GEC汽轮发电机公司在对大型汽轮机设计的最新进展进行推论时认为，在不远的将来提高目前的蒸汽温度和采用两次再循环的可能性都很小，因而循环热效率不会有很大提高。现代叶型效率已达到如此水平，将来的改进仅是有限的。可预见的主要发展是继续增大机组容量来提高经济性。研制更大的转子、更长的末级叶片和更大的排汽环状面积会有利于提高单机容量。然而，如果提高单机功率，采用大转子长叶片的同时，可靠性达不到较小机组的同样水平，对大型汽轮机的考虑就会降低，因为发电厂大型汽轮机设备的启停费用是极大的。目前国外先进大机组十分重视根据热应力状态和疲劳寿命自动控制机组的动态运行，提高机组动态运行的经济性。美国EPRI早在十年前就提出其进一步的研究方向，其中前两项是：(1) 对机组动态特性进行在线监视；(2) 机组起停过程中的经济性能分析，调峰机组的性能监测和机组运行寿命分析。1.2.2 国内发展现状目前，我国大部分电力由火力发电机组承担，其中将近三分之一的火力发电机组，运行时间已接近或超过其设计年限，属于老化设备，从目前我国用电负荷现状来看，仍需这机组继续发挥效益，即使在发电量相当大的发达国家如美、苏、英等国家，也都在努力探讨如何对这些设备进行寿命管理，进一步实现根据设备状态进行的维修管理。从运行方面看，这些火电机组原先承担的基本负荷已让位给大容量火电机组和核电机组，它们要参与调峰运行。本来由于多年运行，转子因热疲劳、蠕变等因素产生了损伤积累和材质退化现象，那么现在的工作条件更为苛刻，势必又要加速其时效老化过程。为确定这些老化设备能否安全的继续运行，必须对老化设备进行剩余寿命监测和预测，制订更新失效老化部件的计划，实施适当的维修管理，这对改善这些火电机组的老化状况，延长运行寿命，是非常重要的.目前有关机组调峰运行过程中的热应力变化和寿命管理方面还有若干问题没有彻底解决。例如在进行机组非稳态温度状态和热应力计算中蒸汽参数和换热系数的确定，寿命预测中我国转子用钢高温疲劳曲线的确定，都有待进一步的研究和完善。汽轮机在高温、高压和高转速的条件下工作，实际运行中参数的变动、负荷的波动与设计工况差别很大，若用理论值和设计值来分析汽轮机的热应力和寿命损耗，很难真实的反映机组的实际状况。只有用实测参数来进行分析计算，才能保证其结果的真实可信。但计算中许多所需要的参数，实地的测取有一定的困难，必须根据运行的实际情况来进行合理的处理。我国从上世纪80年代初开始进行转子寿命损耗预测和寿命分配研究。多年来，我国有关研究机构、高等院校以及制造部门、电厂针对机组调峰的需要，以国产机组为对象，研究了汽轮机主要零部件在非稳态下的温度及热应力分布、变化规律、金属材料的疲劳特性以及部件的寿命损耗。对国产大容量机组参与调峰运行的可行性、调峰运行的安全性和经济性、调峰幅度进行了深入的探讨，对低负荷、少汽无功和两班制等不同的调峰方式在经济性和安全性方面进行了理论分析和实验研究，很多单位都相继开展了汽轮机转子应力监测和寿命损耗计算的研究工作。清华大学的黄仙、倪维斗等提出汽轮机转子热应力的复频法建模，根据转子热应力的形成机理，通过建立“蒸汽温度变化率转子体全温差与“转子体全温差转子体内外表面热应力的复频域数学模型，然后由控制理论方便的得到“蒸汽温度变化率转子体内外表面热应力”的数学模型。倪维斗、黄以腾等在原有汽轮机转子热应力解析计算式的基础上，给出了递推求解热应力值的算法。利用递推式，不仅加快了计算速度，还使得在监控过程中用来保存历史信息所需的存储空间维持不变。支小牧、寇可新建立了有中心孔转子热应力解析递推计算法。该方法计入了换热系数的变化，对汽温非线性变化采用分段线性叠加的分析方法，这个对于边界介质温度变化和换热情况复杂的有中心孔的受热部件的应力计算和研究具有普遍意义，并应用雨流法实时监测和管理转子的寿命损耗。华中科技大学的王坤、黄树红等对125MW汽轮机转子的启停调峰实验进行了研究，计算了启停调峰过程中的应力谱和低周疲劳寿命损耗，探讨了如何对转子寿命损耗进行优化控制的问题。刘华堂、李树人等对国产200MW汽轮机参加调峰运行的寿命进行了分析计算。浙江大学的陈坚红、阮伟等采用大型结构有限元分析软件，对冷态启动温升率对转子寿命的影响进行了研究。郭乙木、鲁祖统在工业汽轮机转子材料低周疲劳实验及对实际工况下转子瞬态温度场、应力场分析的基础上，对转子的有效寿命进行了预测旧。华北电力大学的张光、张宝衡等提出一种计算大型汽轮机转子热应力的二维离散模型，这种二维离散模型也可以进一步被用来计算转子轴向热膨胀和监测机组胀差。黄仙、杨昆等对汽轮机热应力自适应模型进行了研究。崔健、安利强等对转子热应力计算中的有限元方法的精度问题进行了探讨。苑素玲、葛永庆等应用MATLAB软件中SIMULINK软件包对汽轮机转子热应力进行了实时仿真，得到实时监控下汽轮机转子热应力对温升率的仿真曲线。安利强、王璋奇提出了有限元热应力分析的前后处理可视化的实现方法。刘静静、杨昆通过直接测出的高中压缸内缸壁的温度来模拟高中压转子表面温度，作为第一类边界条件，导出了汽轮机转子温度和热应力的计算公式。安江英、杨昆、张保衡对汽轮机转子长期在高温高压状态下材质软化现象进行了研究，并建立了硬度与蠕变和应力松弛性能之间的解析关系式。西安交大对汽轮机转子抗疲劳的模糊可靠性理论进行了分析，主要对汽轮机叶片的疲劳寿命预测方法进行了研究。国家电力公司热工研究院也对汽轮机在变负荷运行的动叶片和疲劳寿命特性进行了相关的研究工作。2002年云南电力实验研究院与上海成套设备研究所合作研究了“汽轮机高压转子在蠕变和疲劳作用下寿命的设计和评定体系，为汽轮机高压转子的寿命设计、寿命评定和寿命延长奠定了坚实的基础2。近年来，我国电厂节能工作受到高度重视。增加高参数大容量机组的装机容量，汽轮机通流部分的改造，利用各种节能理论对火电厂热力系统进行节能诊断等等，都取得显著的成效。与此同时，在火电大机组参与电网两班制调峰和改善大机组动态运行的安全性与经济性方面，也做了大量工作。有关专家学者在该领域的理论研究中，以汽轮机转子疲劳寿命为重点，对国产汽轮机转子冷、热态启动及停机工况热应力及疲劳寿命的计算方法与寿命合理分配原则，国产汽轮机转子不同启停工况下温度场、应力场、寿命损耗率分析，国产机组调峰性能评价，不同调峰运行方式安全经济性分析等等，进行了多方面的开拓。在汽轮机转子热应力和疲劳寿命在线监测方面，我国自行研究的监测装置己在国产200MW汽轮机组上投入运行，它使运行人员直观地看到各种参数的连续变化趋势，迅速地得出正确地判断及时指导运行操作，并兼有运行参数巡检、记录打印等功能，受到运行人员的欢迎。我国力行业中，在努力实现机组最佳动态运行，充分发挥其经济效益上己开始取得共识，并从来的起步状态逐步走向更加深入的探索。我国按转让技术制造和配套的机组，具有上世纪80年代初的世界水平。其中300MW和600MW机组均可实现集中控制，用计算机监视。报警和计算，运行人员可在控制室内实现机组的启动，停机、正常运行以及一些事故处理，整套热控系统都采用了国际上比较先进的技术。与国产大机组相比，有经济性较高，安全性好，调峰性能较强以及自动化程度高等特点，使我国大型火电机组从设计制造到运行管理都有较大的提高，带动了全行业的快速发展。1.3 本文研究内容寿命管理是指以机组经济地实现其服役全寿命为目标，在对设备状态进行监测和评估的基础上优化设备运行与维修管理的新技术:汽轮机转子是发电机组的关键部位，汽轮机转子的寿命问题是机组寿命管理的主要问题3。汽轮机在启停或负荷变动时，转子金属内部将产生较大的温度梯度并由此而产生热应力，这种过度状况下的热应力是影响机组对负荷变动的适应能力和机组寿命损耗的重要因素。尽管大型汽轮机转子的寿命管理技术已在世界范围内取得了很大的进展，技术日见成熟，但仍存在许多问题需要从理论实验方面加以研究。这些问题主要集中于转子寿命消耗的计算方法和寿命预测的方法，比较焦点的问题包括:(1)有限单元法虽然被广泛采用于温度场和应力场的计算，但如何更精确更快速的进行计算和分析以建立更有效的计算模型仍有待研究，包括其他数值计算方法的运用研究和实际软件应用研究；(2)材料的低周疲劳曲线的实验研究和合理运用；(3)材料老化(包括软化、脆化)对转子低周疲劳损耗、蠕变损耗和裂缝演变的影响，从而对转子寿命计算和断裂预测的影响；(4)低周疲劳损伤和蠕变损伤的交互影响和他们之间耦合计算方法，以及对寿命预测的影响；(5)损伤力学的运用，可靠性研究以及断裂力学的发展。针对这些潜在的研究发展方向和问题，我的课题的主要研究内容是： (1)滑参数停机转子热应力数学模型的建立。 1)有中心孔转子外表面、中心部位应力计算公式的推导(无中心孔转子应力理论计算公式的推导)。 2)转子有限元建模。 3)转子温度场的计算，转子应力场的计算。(2)考虑转子复杂应力状态下低周疲劳寿命模型的建立并编写程序。(3)开发以转子热应力约束条件的超超临界汽轮机滑停优化的计算程序。(4)于汽轮机转子在实际运行中总是出现疲劳破坏，最后导致断裂，这将影响整个机组的运行。根据上述计算和分析结果，推断出汽轮机在实际运转过程中由于温度变化引起的应力集中的主要部件，提出一些有建设性的改进，延长汽轮机的运行寿命。(5)于汽轮机转子在实际运行中总是出现疲劳破坏，最后导致断裂，这将影响整个机组的运行。根据上述计算和分析结果，推断出汽轮机在实际运转过程中由于温度变化引起的应力集中的主要部件，提出一些有建设性的改进，延长汽轮机的运行寿命。2. 温度场、应力场模型的建立2.1 概述 汽轮机在启动、停机或负荷变动时，转子金属内部将产生较大的温度梯度并由此而产生热应力，这种过度状况下的热应力是影响机组对负荷变动的适应能力和机组寿命损耗的重要因素。由热弹性理论可知，部件中的热应力与部件金属梯度成正比。温度梯度是由于汽轮机启动、停机或负荷变动过程中，部件被连续加热或冷却，使部件处于热不稳定状态而造成的，汽轮机在稳态工况运行时，不存在或者只有轻微径向温度梯度和热应力。汽轮机转子的前轴封段和前几级，在启停级负荷变动过程中，气温的变化最为剧烈，导致该部位温度梯度和热应力最大，成为最危险的部位，由于设计时间有限，本文只研究了汽轮机高压转子的调节级叶轮根部的温度场及应力场的分布情况，由于汽轮机转子是高速旋转的部件，目前尚无直接测量其金属温度和热应力的有效手段，所以需要理论计算来解决。然而，对于汽轮机转子这样的几何形状不规则、边界条件变化较大的轴对称体，要获得精确的计算结果，需要应用数值解法。本文所应用的有限元法是数值分析方法的一种，它将连续弹性体加以离散化，把受载物体看成是由各个单元所组成的网络。单元的形状可以是各种各样，本论文所介绍的是三角形单元，各单元在他们的顶点处相连接，这些连接点称为结点。有限元法在处理具有复杂形状和条件的物体时有较好的效果，而且单元形状和稀疏程度可以是任意变化的，它可以用较少的节点而使区域达到更好的近似，因此是较好的分析转子的温度场和应力场的方法。对于转子来说，可以将转子考虑为轴对称二维计算模型，将复杂的几何形状连续体离散化，并将介质对转子表面的放热系数和材料物性作为随启动时间的变数来处理。有限元法是最常用于热应力计算的数值解法，它是在能量法和差分法的基础上发展而来的，兼具了两者的优点。有限元法的实质是通过单元内的变分计算并总体合成代替微分方程求解，只要划分单元足够小，就能用简单的线性插值函数得到足够精度的结果，而且网格划分灵活，对不规则边界的适应性好。2.2 二维离散温度场模型2.2.1 温度场模型的建立计算转子不稳定温度场时，将转子视为一个均匀、各向同性、无内热源的物体，属于轴对称非定常温度函数，温度T(z，r，r)在区域D中满足下列偏微分方程式： (2.1)式中 l 转子钢材的导热系数，W/(m)； r 材料的密度，kg/m3； Cp 材料的比热，kJ/(kg)。确定上面微分方程式(2.1)的解，需要满足的初始条件。，此外还应满足一定的边界条件。对于汽轮机转子来说，外表面的边界条件由蒸汽对转子表面的换热速度来确定，属于传热学中的第三类边界条件，即边界与介质的热交换条件为已知: (2.2)式中 Tf 与转子表面接触的气温； a 蒸汽与转子表面的放热系数。由变分原理，公式(2-1)的第三类边值问题可以等价转化为求一泛函数J的极值问题即：。采用有限单元法对区域D进行离散化，在分成的小单元上以上泛函数式是成立的。这样尽管T在D中的分布很复杂，但在每个小单元可以近似看作是线性分布。由于D是全部单元的总体，所以： (2.3)式中 Je 对单元e进行计算； E 在区域D内划分单元总数； n 待求温度值的节点个数。边界单元的泛函数形式：+ (2.4) 内部单元泛函数形式：若假设温度T(z，r，)在单元e中呈线性分布，在单元作变微分计算，并经过一系列的处理，最后得： (2.5)式中 i、j、k 三角形单元三个顶点编号，如图2.1所示。 ij j k图 2.1 i、j、k三角形单元三个顶点编号将式(2.5)带入式(2.4)中，可得到n阶线性代数方程组，以求出n个节点的温度，其形式为： (2.6)式中 刚度矩阵(nn阶正定对称方阵)； 非稳态温度矩阵(nn阶正定对称方阵)；、和分别为温度、热载荷和温度变化率列向量(n1阶方阵)。其中 项用有限元差分法展开。本文通过对各种差分格式的比较，选用伽辽金格式，则式(2.6)变为： (2.7)式(2.7)稳定且不振荡的条件是：0对于轴对称有限元分割，可近似地得出第一类边界条件时时间步长Dt 限制条件：式中Dx 单元的平均边长。伽辽金格式的特点是在超过上式的范围时发生的震荡比较轻微。详细的推导在文献上都有详细的介绍，由于篇幅所限，不再详述。二维对称温度场的计算精度关键在于正确处理边界条件。温度场计算过程中，根据边界位置及当时的蒸汽温度和压力，实时地计算该处蒸汽的各项物性参数及换热系数。温度场计算分为稳态温度场计算和瞬态温度场计算4。2.2.2 各部分放热系数的计算在启停过程中，转子外表面上蒸汽温度和放热系数随不同的轴向位置和时间而变化，亦即是时间和空间的函数。有关汽轮机转子表面放热系数的计算公式，国内外尚未统一定论，放热系数对热应力的计算影响不大，当启动未进入准稳态时，不同计算公式所得热应力有一定差值，进入准稳态之后，放热系数的大小对热应力值无明显影响。本文采用哈汽一南工公式，转子各部位放热系数公式如下：(1)叶轮两侧的放热系数： W/(m2) (2.8)式中 Nu 努赛尔特征数，可按下式计算：当Re2.4时，Nu=0.675当Re时， 式中 n 气流的运动粘度系数，m2/s； u 叶轮外圆Rb处得圆周速度，m/s； Rb 叶轮的外圆半径，m； Re 介质雷诺数； l 蒸汽导热系数，W/(m)。(2)汽封处的放热系数(曲齿)： W/(m) (2.9)式中 d 汽封间隙，mm； S 汽封齿距，mm； H 转子表面到汽封环的高度，m； wi 缝隙中气流的平均速度，m/s。(3)光轴处的放热系数： W/(m) (2.10)式中 Ra 光轴外圆半径，m。(4)轮缘处的放热系数： W/(m) (2.11)式中 l 叶片材料的导热率 ； r0 叶轮顶部的半径，m。也可以近似处理，高、中压转子叶轮外缘处的蒸汽放热系数349 W/(m)、160 W/(m)【5】。(5)调节级两侧的放热系数： 当Re2.4时， 当Re时，式中 r0 叶轮半径，m；W 叶轮旋转角速度，rad/s； m 叶轮切向速度，m/s；r 叶轮的平均半径，m； u 蒸汽的粘度，kg/(ms)；Cp 蒸汽的平均比热，kJ/(kg)； r1 蒸汽的密度，kg/m3；n 运动粘度，kg/(ms)。2.3 温度场的变量说明及运行程序温度场程序框图如图2.1所示，采用Basic语言对程序进行编制。标识符说明：(1)简单变量A8表示单元边界上的对流换热系数a(W/m2)；A9表示单元边界上对流换热系数a和介质温度Tf()的乘积；BI，BJ，BM分别表示公式(5.19)中的bi、bj、bm；B6，B7分别第二、三类边界单元的分组数；CI，CJ，CM分别表示公式(5.19)中的ci、cj、cm；C4，C5分别表示图纸尺寸放大倍数和时间变步长倍数；开 始形成温度总刚度矩阵求解温度矩阵求解节点方程组计算单元，节点温度，计算体积平均温度输出某一时刻节点温及体积平均温度判断是否计算温度场温度场结 束判断是否计算稳定厂IST=？0输入单元信息，节点坐标，赋初值输入各部位压力，温度及变化率判断是否结束温度场JST=JST+1JST=0图2.1 温度场程序图D1表示三角形单元的面积；分别表示内部、第一类和第二类边界单元的最后一个编号； 表示单元总数，也是第三类边界单元的最后一个编号；表示公式中的；判别量，和2分别作稳态和非稳态温度场计算；判别量，只作单元三结点编号信息的分解；判别量，和1分别表示轴对称问题和平面问题；分别表示单元顶点i，j，m的编号；表示结点总数；表示第一类边界单元的边界结点数；分别矩阵的行号和列号；分别表示组成单元材料的总数和该单元的材料类别；判别量，分别表示向后差分格式，C-N格式和伽辽金格式；表示公式中的乘积；分别表示第二类和第三类边界单元的数目；RI，RJ，RM分别表示单元三顶点的纵坐标；表示系数矩阵下三角阵中非零带内元素的总数；SI表示边界单元的边长；表示物体的均匀初温()；表示初始时间步长Dt0(s)；表示非稳定导热过程所进行的时间(s)；表示变时间步长，若作定时间步长计算，则；表示规定时间推进积分的总次数；表示时间推进积分的累计次数，为初始时刻；XI，XJ，XM分别表示单元三顶点的横坐标。(2)组数 用来保存P数组在非稳态导热计算中的初值；用来保存E数组在非稳态导热计算中的初值；表示第二类边界单元各分组中的最后一个单元号，并记用来保存K数组在非稳态导热计算中的初值；表示式(5.63)中N矩阵下三角元素一维压缩排列的数组；表示第三类边界单元的对流换热系数；表示第三类边界单元的周围介质温度；表示第三类单元信息的数组；表示第三类边界单元各分组的周围介质温度Tf()；表示第三类边界单元各分组的对流换热系数a(W/m2)；表示公式(5.63)中K矩阵下三角元素一维压缩排列的数组；表示第一类边界单元的边界结点温度()；表示第二类边界单元的热流密度；表示主对角元素在一维压缩排列中的地址；表示第一类边界单元的边界的结点号；表示线性代数方程组的右端列向量，最后存放解向量；表示单元的单位体积热产生率；表示结点的纵坐标y或r，输入单位用；作为临时存放非零带内元素的一个数组；存放非稳态导热前一时刻的温度场；表示第二类边界单元各分组的热流密度，流出取正值，流入取负值；表示第三类边界单元各分组中的最后一个单元，并记；表示结点的x坐标(横坐标)，输入单位用；表示单元的导热系数(W/m2)；表示单元的容积比热。(3)子程序DIV用来分解单元信息；ADDRESS用来确定主对角元素在一维压缩排列中的地址；JUDGE矩阵消元或相乘时为查明零元素而用的判据；SOALG用消元法求解线性代数方程组；TH11用来计算系数矩阵元素及内热源对右端列向量的贡献；BOUND用来计算边界条件对右端列向量的贡献；DIFF表示非稳态导热计算中的三种时间差分格式；NSTEA作式(5.84)中非稳态导热右端列向量的加法和惩罚运算。10 REM TEMPERATURE FIELD SOLVED BY FEM12 READ L0，L1，E0，E1，E2，E3，M0，C4，B6，B7，H0，T0，T1，C5，U0 20 Q2=E2-E1:Q3=E3-E2 22 DIM A(L0)，N(L0，)P(L0)，R(L0)，S(L0)，T(L0)，X(L0)24 DIM F(Q3)，G(Q3).H#(E3)，I(B7)，J(B7)，W(B7) 26 DIM Q(M0)，Y(M0)，Z(MO)40 FOR A0=1 TO L0:READ X(A0)， R(A0)42 PRINT X(A0);R(AO):NEXT A044 FOR A0=1 TO E3:READ H#(A0) 46 PRINT H#(A0):NEXT A048 FOR A0=1 TO M0:READ Y(A0)，Z(A0)，Q(A0)50 PRINT Y(A0)，Z(A0)，Q(A0)52 NEXT A0 54 IF Ll=0 THEN 62 56 PRINTST B.C.NODE NUMBER AND TEMPERATURE58 FOR A0=1 TO Ll:READ，O(A0)，L(A0) 60 PRINT O(A0).L(A0):NEXT A0 62 IF Q2=0 THEN 76 64 PRINT2ND B.C.GROUP LAST ELEMENT NUMBEP.& HEAT FLUX66 FOR A0=0 TO B6:READ C(AO)，V(AO) 68 PRINT C(A0);V(A0):NEXT A0 70 FOR A0=1 TO B6:V(A0)=V(A0)*(.000001)72 FOR E9=C(A0-l)+l TO C(A0) ，74 M(E9-El)=V(A0):NEXT E9:NEXT A076 IF Q3=0 THEN 94 78 PRINT3RD B.C.GROUP LAST NO.FILM COEF. & MEDIUM TEMP80 FOR A0=0 TO B782 READ W(A0)，J(A0)，I(A07)84 PRINT W(A0);J(A0); I(A0 ) :NEXT A0 86 FOR A0=1 TO B7:J(A0):J(A0)*.00000l 88 FOR E9=W(A0-l)十1 TO W(A0)90 F(E9-E2)=J(A0):G(E9-E2)=I(A0)92 NEXT E9:NEXT， A094 FOR A0=1 TO L0 96T(A0)=T0:NEXT A0 98 IF HO=0 THEN 104100 PRTNTPLANE TEMPERATURE PROBLEM 102 GOTO 106104 PRINTAXIS SYMMEIRY TEMPRATURE PROBLEM106IF ABS(C4-1)99999! THEN 126124 Z(A0)*(1E-09)126 NEXT A0128 T2=0:T7=T1:U9=1 130 FOR A0=1 TO S0:K(A0)=0:E(A0)=0:NEXT A0 132 FOR A0=l TO LO:P(A0)=0:NEXT A0134FOR E9=l TO E3:G1=2:GO SUB 1000 13b IF Dl0 THEN l44 138 PRINTFAIL IN ELEMENT NUMBER E9= ;E9140 GOTO 260144 IF H0=0 THEN 148146 RI=I:RJ=1:RM=1，148 F0=Y(M2)*(RI+RJ+RM/(12*D1)152 N2=D1*(M2)153 GOSUB 1250154 I7=I1:J7=J1:B8=BI:B9=BJ156 C8=C1:C9=CJ:R7=RI:R8=BJ158 R9=RM:GOSUB 1200160 I7=J1:J7=M1:B8=BJ:B9=BM162 C8=CJ:C9=CM:R7=RJ:R8=RM164 R9=RI:COSUB 1200166 I7=M1:J7=I1:B8=BM:B9=BI168 C8=CM:C9=CI:R7=RM:R8=RI170 R9=RJ:GOSUB 1200200 NEXT E9202 IF T19999! THEN208 204 G0=2:N0=1:GOSUB1300206 T2=T2+T7:GOSUB1350 208 IF L1=0 THEN 216210 FOR A0=1 TO L1:I0=O(A0)212 K(N(I0)=K(N(I0)*1E+08214 P(I0)=K(N(I0)*L(A0):NEXT A0216 GOSUB 1100:IF T199999! THEN 224218 PRINTSTEADY TEMPRATURE FIELD220 FOR I0=1 TO L0:PRINT P(I0);224 PRINTNONSTEADY TEMPRATURE FIELD226 PRINTTIME=;T2;STEP SIZE=;T7;C5=;C5228 FOR I0=1 TO L0230 PRINTT(;A0;)=;P(A0)232 T(A0)=P(A0):P(A0)234 NEXT A0:U9=U9+1:PRINT236 FOR I0=1 TO S0:K(I0)=D(I0)238 E-(I0)=B(I0):NEXT I0240 IF C51.2 THEN 244242C5=C5-.05 244 T7=T7*C5:IF U9U0+.1 THEN 204260 END 300 DATA 5，0，3，3，3，4，1，1，0，1，1，20，1，1，3302 DATA 0，0，10，0，5，5，0，10，10，10304 DATA 0.0010030041，0.0010020031，0.0030050041306 DATA 0.0030020051，16，1E6，0308 DATA 3，0，0，4，1600，1001000 REM SUBROUTIME OF DIV1002 Z1#=H#(E9)*1000:I1=INT(Z1#)1004 Z1#=(Z1#-I1)*1000:J1=INT(Z1#)1006 Z1#=(Z1#-J1)*1000:M1=INT(Z1#)1008 IF G1J1 THEN R0=J11042 IF R0M1 THEN R0=J11044 IF N(I1)=I1-R0 THEN N(I1)=I1-R0+11050 IF N(J1)=J1-R0 THEN N(J1)=J1-R0+11052 IF N(M1)N(L7)-N(L7-1)-1:NEXT 10901084 IF G01.5 THEN 10881086 S(L7)=E(N(L7)-L5) GOTO 10921088 S(L7)=K(N(L7-L5) GOTO 10921090 S(L7)=01092 RETURN1200 REM SUBROUTINE OF TH111202 A7=N(I7)：K(A7)=K(A7)+F0*(B8*B8+C8*C8) 1204 IF 1999999! THEN 12181216 E(A7)=E(A7)+N2*(3*R7+R8+R9)/301218 IF I7J7 THEN 12221220 A7=N(I7)-(I7-J7): GOTO 12241222 A7=N(J7)-(J7-I7) 1224