热带钢轧机轧辊热变形的快速有限元模型

第36卷第5期 北京科技大学学报 Vol.36 No.5 2014年5月 Journal of University of Science and Technology Beijing May 2014 热带钢轧机轧辊热变形的快速有限元模型 孔繁甫”，何安瑞四，邵健”，葛鲜玉) 1)北京科技大学高效轧制国家工程研究中心，北京1000832)中国第二重型机械集团公司设计研究院，德阳618013 ☒通信作者，E-mail:harui@ustb.cedu.cn 摘要分析了解析法与简化方法的缺点，采用有限元法建立轧辊热变形计算模型.针对轧制过程中轧辊热变形计算数据特 点，将计算任务分为负责静态数据准备的预计算和负责动态数据准备与热变形求解的更新计算，换辊时进行预计算，计算任 意时刻的热变形时只需进行更新计算，计算量远小于标准有限元程序.根据特殊处理的计算流程，编写了基于轴对称有限元 法的轧辊热变形程序，其计算结果与ANSYS结果一致，精度均高于简化方法约30%.自编轧辊热膨胀有限元程序计算精度 高，耗时少，满足在线热膨胀预报要求. 关键词带钢轧机：热轧：板形控制：轧辊：热变形：有限单元法 分类号TG333.17 Finite element model for rapidly evaluating the thermal expansion of rolls in hot strip mills KONG Fanfu",HE An-rui,SHAO Jian,GE Xian-yu? 1)National Engineering Research Center for Advanced Rolling,University of Science Technology Beijing,Beijing 100083,China 2)Design Research Institute,China National Erzhong Group Co.,Deyang 618013,China Corresponding author,E-mail:harui@ustb.edu.cn ABSTRACT The disadvantages of the analytical solution and the simplified method were analyzed and a thermal expansion model was established by the finite element method.Based on the data characteristics of the rolls'geometry in a rolling unit,the calculation task was divided into pre-ealculation responsible for static data preparation and renewal calculation responsible for renewing dynamic data and thermal expansion.The pre-calculation was called only after roll change and the renewal calculation was called when the thermal expansion data was needed to be updated.This special tailored flow contributes to a much smaller computation complexity in compari- son with standard finite element analysis.Furthermore,the corresponding program was compiled according to the tailored flow and its results are in agreement with those by ANSYS,with the accuracy of about 30%higher than the simplified method.Due to high accura- cy and short calculation time,this home-made finite element analysis program is promising to be implemented in online applications for evaluating the thermal expansion of rolls. KEY WORDS strip mills:hot rolling:shape control:rolls;thermal expansion:finite element method 在热轧带钢生产过程中，工作辊的热变形与辊 础部分之一，其计算精度对产品板形质量有直接 缝凸度在同一数量级，对带钢凸度及平坦度的影响 关系. 非常显著.同时，轧辊的热变形是影响辊间接触压 在轧辊热辊形的研究中，文献09]都将轧辊 力分布的主要因素之一，对工作辊挠曲有重要影响. 温度场与热变形放在一起研究，且研究重点皆为温 因此，轧辊热辊形的计算模型是热轧板形模型的基 度场部分，轧辊热变形计算均采用简化方法.陈宝 收稿日期：201302-18 基金项目：新世纪优秀人才支持计划资助项目(NCET-100223):中央高校基本科研业务费专项资金资助项目(FRF-TP-11O03A) DOI:10.13374/j.issn1001-053x.2014.05.016:http:/journals.ustb.edu.cn

第 36 卷 第 5 期 2014 年 5 月 北京科技大学学报 Journal of University of Science and Technology Beijing Vol． 36 No． 5 May 2014 热带钢轧机轧辊热变形的快速有限元模型 孔繁甫1) ，何安瑞1) ，邵 健1) ，葛鲜玉2) 1) 北京科技大学高效轧制国家工程研究中心，北京 100083 2) 中国第二重型机械集团公司设计研究院，德阳 618013  通信作者，E-mail: harui@ ustb． edu． cn 摘 要 分析了解析法与简化方法的缺点，采用有限元法建立轧辊热变形计算模型． 针对轧制过程中轧辊热变形计算数据特 点，将计算任务分为负责静态数据准备的预计算和负责动态数据准备与热变形求解的更新计算，换辊时进行预计算，计算任 意时刻的热变形时只需进行更新计算，计算量远小于标准有限元程序． 根据特殊处理的计算流程，编写了基于轴对称有限元 法的轧辊热变形程序，其计算结果与 ANSYS 结果一致，精度均高于简化方法约 30% ． 自编轧辊热膨胀有限元程序计算精度 高，耗时少，满足在线热膨胀预报要求． 关键词 带钢轧机; 热轧; 板形控制; 轧辊; 热变形; 有限单元法 分类号 TG 333. 17 Finite element model for rapidly evaluating the thermal expansion of rolls in hot strip mills KONG Fan-fu1) ，HE An-rui1)  ，SHAO Jian1) ，GE Xian-yu2) 1) National Engineering Ｒesearch Center for Advanced Ｒolling，University of Science ＆ Technology Beijing，Beijing 100083，China 2) Design ＆ Ｒesearch Institute，China National Erzhong Group Co． ，Deyang 618013，China  Corresponding author，E-mail: harui@ ustb． edu． cn ABSTＲACT The disadvantages of the analytical solution and the simplified method were analyzed and a thermal expansion model was established by the finite element method． Based on the data characteristics of the rolls’geometry in a rolling unit，the calculation task was divided into pre-calculation responsible for static data preparation and renewal calculation responsible for renewing dynamic data and thermal expansion． The pre-calculation was called only after roll change and the renewal calculation was called when the thermal expansion data was needed to be updated． This special tailored flow contributes to a much smaller computation complexity in compari￾son with standard finite element analysis． Furthermore，the corresponding program was compiled according to the tailored flow and its results are in agreement with those by ANSYS，with the accuracy of about 30% higher than the simplified method． Due to high accura￾cy and short calculation time，this home-made finite element analysis program is promising to be implemented in online applications for evaluating the thermal expansion of rolls． KEY WOＲDS strip mills; hot rolling; shape control; rolls; thermal expansion; finite element method 收稿日期: 2013--02--18 基金项目: 新世纪优秀人才支持计划资助项目( NCET--10--0223) ; 中央高校基本科研业务费专项资金资助项目( FＲF--TP--11--003A) DOI: 10． 13374 /j． issn1001--053x． 2014． 05． 016; http: / /journals． ustb． edu． cn 在热轧带钢生产过程中，工作辊的热变形与辊 缝凸度在同一数量级，对带钢凸度及平坦度的影响 非常显著． 同时，轧辊的热变形是影响辊间接触压 力分布的主要因素之一，对工作辊挠曲有重要影响． 因此，轧辊热辊形的计算模型是热轧板形模型的基 础部分之一，其计算精度对产品板形质量有直接 关系． 在轧辊热辊形的研究中，文献［1--9］都将轧辊 温度场与热变形放在一起研究，且研究重点皆为温 度场部分，轧辊热变形计算均采用简化方法． 陈宝

第5期 孔繁甫等：热带钢轧机轧辊热变形的快速有限元模型 ·675· 官等0、孔祥伟等四和高建红☒采用有限元法计 空间轴对称的热弹性力学问题，其几何方程和平衡 算轧辊温度场与热变形，均使用商用有限元软件 方程均与一般轴对称弹性问题一致，而物理方程有 (ANSYS,NASTRAN),并着眼于热变形的离线分析 所区别.热弹性变形的物理方程中，应变一部分由 研究而非在线预测.当前，在边界条件准确的情况 应力引起，一部分由温升引起.求解空间轴对称热 下，轧辊温度场计算模型基本成熟，二维差分模型的 弹性问题时，采用位移法表示平衡方程可得其基本 计算速度与精度能达到在线应用水平，德国 微分方程： Siemens、日本TMEIC等公司提供的板形控制系统 中，轧辊温度场模型均为二维差分模型. u- 1e=2(1+以.0 7+1-2m=1-2m dr' 然而，在轧辊温度场己给定的情况下，解析法求 (2) 解热变形困难。工程上所采用的简化计算模型没有 1e-20+》a9 0+1-2'=1-2w dz 考虑热应力对热变形的轴向平滑作用，存在较大计 式中，“为径向位移，为轴向位移，v为材料泊松 算误差.对计算结果进行格林函数滤波)或跑偏 概率滤波的方法没有从引起计算误差的根 比，a为热弹性系数，V为Laplace算子，e为体 应变 源一热应力的牵制作用的角度考虑，均为唯象方 上述微分方程为线性非齐次微分方程，其解由 法，缺乏理论支撑 以下两部分组成：非齐次微分方程的特解和相应齐 本文研究在给定温度场作用下轧辊的热变形问 次微分方程(Laplace方程)的通解.其中，特解可采 题.首先分析解析法和传统热辊形简化计算方法的 用热弹性位势函数Φ来求得，通解可采用乐甫位移 缺点及原因，然后针对轧制过程中轧辊热变形计算 函数L求得，基本微分方程转化为： 的数据特点，建立了基于特殊流程的轴对称有限元 模型，最后采用自编有限元程序计算了给定温升下 2 1+卫a0: 的轧辊热变形.高精度与高效率的热变形计算程 1-v (3) 序，可实现轧辊热膨胀在线计算，推动板形设定模型 22L=0. 精度的提高. 因此，对于轴对称热弹性问题，找到适当的调和 1 轧辊热弹性变形模型 函数Φ及重调和函数L使得由此给出的位移分量 和应力分量能够满足边界条件，就得到该问题的正 在轧制过程中，轧辊温度场是一个非稳态系统， 确解答.然而，上述解析解求解困难。乐甫位移函 轴向、径向和周向温度均随着轧制过程而变化.然 数的通解表达式由一系列Bessel函数组成，计算繁 而根据Stevens等n对工作辊的温度进行测量的发 琐.同时，由于工作辊的温度分布复杂、局部变化剧 现：在轧辊表层2~3mm深度内，温度呈周期性剧烈 烈，需要用高次二元多项式描述，所以特解需要用更 变化，而在轧辊内99%以上的部分，温度分布基本 高次二元多项式来描述，求解过程极其复杂 上是轴对称的.因此，在研究轧辊温度场与热变形 时，可以忽略轧辊温度场沿周向的变化，将轧辊温度 2工程简化计算方法 场简化为空间轴对称模型（见图1所示），温升的柱 由于解析解求解困难，众多学者习对轧辊热 坐标描述为 变形的研究及板形设定工程应用均将轧辊热变形模 ⊙(r,)=T(r,z)-T (1) 型进一步简化，采用温升分布与轴向无关的无限长 式中，和z为径向和轴向坐标，T。和T分别为初始 柱体的解析结果.得到轧辊温升⊙（π，z)之后，采用 温度和任意时刻温度. Timoshenko公式计算轧辊表面径向位移： Ar u(2l,=R=2(1+）.g. (r,2)rdr. (4) R T(r.:) u(r.= 式中，R为轧辊半径 在圆柱体的温升只与径向相关时，采用以上方 图1轧辊轴对称简化模型 法计算得到的热变形量是准确的.然而，轧辊温升 Fig.I Axial symmetrical simplified model for rolls 在轴向有较大变化，这种计算方法没有考虑热应力 轧辊热变形模型与温度场模型一致，也简化为 对热变形的平滑作用，因此计算结果会高估热变形

第 5 期 孔繁甫等: 热带钢轧机轧辊热变形的快速有限元模型 官等［10］、孔祥伟等［11］和高建红［12］采用有限元法计 算轧辊温度场与热变形，均使用商用有限元软件 ( ANSYS，NASTＲAN) ，并着眼于热变形的离线分析 研究而非在线预测． 当前，在边界条件准确的情况 下，轧辊温度场计算模型基本成熟，二维差分模型的 计算 速 度 与 精 度 能 达 到 在 线 应 用 水 平，德 国 Siemens、日本 TMEIC 等公司提供的板形控制系统 中，轧辊温度场模型均为二维差分模型． 然而，在轧辊温度场已给定的情况下，解析法求 解热变形困难． 工程上所采用的简化计算模型没有 考虑热应力对热变形的轴向平滑作用，存在较大计 算误差． 对计算结果进行格林函数滤波［13］或跑偏 概率 滤 波［14］ 的方法没有从引起计算误差的根 源———热应力的牵制作用的角度考虑，均为唯象方 法，缺乏理论支撑． 本文研究在给定温度场作用下轧辊的热变形问 题． 首先分析解析法和传统热辊形简化计算方法的 缺点及原因，然后针对轧制过程中轧辊热变形计算 的数据特点，建立了基于特殊流程的轴对称有限元 模型，最后采用自编有限元程序计算了给定温升下 的轧辊热变形． 高精度与高效率的热变形计算程 序，可实现轧辊热膨胀在线计算，推动板形设定模型 精度的提高． 1 轧辊热弹性变形模型 在轧制过程中，轧辊温度场是一个非稳态系统， 轴向、径向和周向温度均随着轧制过程而变化． 然 而根据 Stevens 等［15］对工作辊的温度进行测量的发 现: 在轧辊表层 2 ～ 3 mm 深度内，温度呈周期性剧烈 变化，而在轧辊内 99% 以上的部分，温度分布基本 上是轴对称的． 因此，在研究轧辊温度场与热变形 时，可以忽略轧辊温度场沿周向的变化，将轧辊温度 场简化为空间轴对称模型( 见图 1 所示) ，温升的柱 坐标描述为 Θ( r，z) = T( r，z) － T0 ． ( 1) 式中，r 和 z 为径向和轴向坐标，T0 和 T 分别为初始 温度和任意时刻温度． 图 1 轧辊轴对称简化模型 Fig． 1 Axial symmetrical simplified model for rolls 轧辊热变形模型与温度场模型一致，也简化为 空间轴对称的热弹性力学问题，其几何方程和平衡 方程均与一般轴对称弹性问题一致，而物理方程有 所区别． 热弹性变形的物理方程中，应变一部分由 应力引起，一部分由温升引起． 求解空间轴对称热 弹性问题时，采用位移法表示平衡方程可得其基本 微分方程: 2 Δ u － u r 2 + 1 1 － 2ν ·e r = 2( 1 + ν) 1 － 2ν ·α·Θ r ; 2 Δ w + 1 1 － 2ν ·e z = 2( 1 + ν) 1 － 2ν ·α·Θ z { ． ( 2) 式中，u 为径向位移，w 为轴向位移，ν 为材料泊松 比，α 为 热 弹 性 系 数， Δ 为 Laplace 算 子，e 为 体 应变． 上述微分方程为线性非齐次微分方程，其解由 以下两部分组成: 非齐次微分方程的特解和相应齐 次微分方程( Laplace 方程) 的通解． 其中，特解可采 用热弹性位势函数 Φ 来求得，通解可采用乐甫位移 函数 L 求得，基本微分方程转化为: 2 Δ Φ = 1 + ν 1 － ν ·αΘ; 2 Δ 2 Δ L = 0 { ． ( 3) 因此，对于轴对称热弹性问题，找到适当的调和 函数 Φ 及重调和函数 L 使得由此给出的位移分量 和应力分量能够满足边界条件，就得到该问题的正 确解答． 然而，上述解析解求解困难． 乐甫位移函 数的通解表达式由一系列 Bessel 函数组成，计算繁 琐． 同时，由于工作辊的温度分布复杂、局部变化剧 烈，需要用高次二元多项式描述，所以特解需要用更 高次二元多项式来描述，求解过程极其复杂． 2 工程简化计算方法 由于解析解求解困难，众多学者［1--9］对轧辊热 变形的研究及板形设定工程应用均将轧辊热变形模 型进一步简化，采用温升分布与轴向无关的无限长 柱体的解析结果． 得到轧辊温升 Θ( r，z) 之后，采用 Timoshenko 公式［16］计算轧辊表面径向位移: u( z) | r = Ｒ = 2( 1 + ν)·α Ｒ ·∫ Ｒ 0 Θ( r，z) rdr． ( 4) 式中，Ｒ 为轧辊半径． 在圆柱体的温升只与径向相关时，采用以上方 法计算得到的热变形量是准确的． 然而，轧辊温升 在轴向有较大变化，这种计算方法没有考虑热应力 对热变形的平滑作用，因此计算结果会高估热变形 · 576 ·

·676 北京科技大学学报 第36卷 量，轴向温升变化越大，误差越大.采用基于二维差 用节点位移表示的单元位移场为： 分法的温度场计算模型时，轴向相邻节点位置的温 [u=∑N,4: 度不可避免出现图2中左图所示差异，按照上述方 (6) 法计算的位移结果将产生右图虚线所示的阶梯，而 lw=∑No 理论实际位移应该是右图中实线表示的连续平滑曲 根据几何关系，得单元应变为 线，且理论实际位移量小于上述方法的计算结果 e=Le,E:Yn so]T= 为避免轴向温升不一致带来的热膨胀计算误 Bu ow r L ar (7) 差，Sauer、Jarrett和Allwood网在整个轧辊内采用 az ar Ritz法求解热变形.Ritz法的精度依赖于位移函数 将式(6)代入式(7)，可得用节点位移表示的单 的选择，轧辊温度分布复杂，当温度在径向表面层及 元应变： E=Bq. (8) 轴向局部变化剧烈时，很难找到在全域内与实际位 移吻合的位移函数 式中：B为单元应变矩阵；q为单元节点位移向 量，即 9=4101山202…山w].(9) 将温升引起的应变看成初始应变，将物理方程 写成矩阵形式，可得 =[o,T o]=D(8-E0).(10) 温度变化 -·一计算位移一实际位移 式中：D为材料矩阵；eo为温升⊙(r,z)所产生的初 图2轴向温差与位移结果示意图 始应变，即 Fig.2 Axial temperature difference and resulting expansion Eo=[aeae 0 a0]T. (11) 3 轴对称有限元方法 因此，单元应变能为： 有限元法的基本求解思想是把复杂的计算域划 U-ionv-f (De-De)dv- 分为规则的分析单元，在单元内选择位移函数来逼 人FDRqdV--gB'De,d业 (12) 近单元中的真解，整个计算域内的近似解由所有单 元解构成。相比而言，有限元法在单元内选择位移 由于单元不受外力作用，因此单元外力功为0. 函数进行逼近，而Rtz法在全域内选择位移函数来 根据虚位移原理，可得 进行逼近.因此，对于复杂问题，有限元法处理灵 'B'DBqdy='B'Deodv. (13) 活、假设条件少、精度高以及容易处理边界条件，逐 通过高斯数值积分可得单元刚度矩阵k。与单 渐成为工程数值分析的有力工具. 元温度载荷向量∫，分别为： 在有限元法中，将温升引起的应变看成初始应 变，即可将热弹性问题等效为一般弹性问题，统一求 k.2B'DBrdeuJdedn: (14) 解.小变形热弹性分析属于线性问题，有限元法计 算流程简单.针对轧辊热变形计算的数据特点，采 .BDeardejdedn (15) 用特殊流程的自编有限元程序可进一步减小计算 式中，J为单元雅可比矩阵 量，实现高效率、高精度计算，计算速度能达到在线 利用单元节点连接关系，可分别由单元刚度矩 应用要求 阵k.和单元温度载荷向量∫集成总体刚度矩阵K 3.1有限元计算模型 和总体载荷向量F,轧辊热变形的总体弹性势能可 轧辊热弹性有限元模型中采用四节点轴对称线 表示为 性等参单元，选取形态函数如下： n=20'K0-Qr (16) N,=0.25(1-)(1-n): 式中，Q为总体位移列阵，K为总体刚度矩阵，F为 N2=0.25(1-)(1+7); (5) 总体载荷向量 N3=0.25(1+)(1+n): 对于轧辊热变形问题，对称轴上节点有径向位 N4=0.25(1+)(1-n). 移约束，中心位置节点有轴向位移约束.采用罚函 式中，、刀为单元局部坐标 数法处理位移约束条件，总体势能为：

北 京 科 技 大 学 学 报 第 36 卷 量，轴向温升变化越大，误差越大． 采用基于二维差 分法的温度场计算模型时，轴向相邻节点位置的温 度不可避免出现图 2 中左图所示差异，按照上述方 法计算的位移结果将产生右图虚线所示的阶梯，而 理论实际位移应该是右图中实线表示的连续平滑曲 线，且理论实际位移量小于上述方法的计算结果． 为避免轴向温升不一致带来的热膨胀计算误 差，Sauer［17］、Jarrett 和 Allwood［18］在整个轧辊内采用 Ｒitz 法求解热变形． Ｒitz 法的精度依赖于位移函数 的选择，轧辊温度分布复杂，当温度在径向表面层及 轴向局部变化剧烈时，很难找到在全域内与实际位 移吻合的位移函数． 图 2 轴向温差与位移结果示意图 Fig． 2 Axial temperature difference and resulting expansion 3 轴对称有限元方法 有限元法的基本求解思想是把复杂的计算域划 分为规则的分析单元，在单元内选择位移函数来逼 近单元中的真解，整个计算域内的近似解由所有单 元解构成． 相比而言，有限元法在单元内选择位移 函数进行逼近，而 Ｒitz 法在全域内选择位移函数来 进行逼近． 因此，对于复杂问题，有限元法处理灵 活、假设条件少、精度高以及容易处理边界条件，逐 渐成为工程数值分析的有力工具． 在有限元法中，将温升引起的应变看成初始应 变，即可将热弹性问题等效为一般弹性问题，统一求 解． 小变形热弹性分析属于线性问题，有限元法计 算流程简单． 针对轧辊热变形计算的数据特点，采 用特殊流程的自编有限元程序可进一步减小计算 量，实现高效率、高精度计算，计算速度能达到在线 应用要求． 3. 1 有限元计算模型 轧辊热弹性有限元模型中采用四节点轴对称线 性等参单元，选取形态函数如下: N1 = 0. 25( 1 － ξ) ( 1 － η) ; N2 = 0. 25( 1 － ξ) ( 1 + η) ; N3 = 0. 25( 1 + ξ) ( 1 + η) ; N4 = 0. 25( 1 + ξ) ( 1 － η)        ． ( 5) 式中，ξ、η 为单元局部坐标． 用节点位移表示的单元位移场为: u = ∑Niui ; w = ∑Niwi { ． ( 6) 根据几何关系，得单元应变为 ε =［εr εz γrz εθ ］T = u r w z u z + w r u [ ] r T ． ( 7) 将式( 6) 代入式( 7) ，可得用节点位移表示的单 元应变: ε = Bq． ( 8) 式中: B 为单元应变矩阵; q 为单元节点位移向 量，即 q =［u1 w1 u2 w2 … un wn］T ． ( 9) 将温升引起的应变看成初始应变，将物理方程 写成矩阵形式，可得 σ =［σr σz τrz σθ ］T = D( ε － ε0 ) ． ( 10) 式中: D 为材料矩阵; ε0为温升 Θ( r，z) 所产生的初 始应变，即 ε0 =［αΘ αΘ 0 αΘ］T ． ( 11) 因此，单元应变能为: U = 1 2 ∫V εT σdV = 1 2 ∫V ( εT Dε － εT Dε0 ) dV = 1 2 ∫V qT BT DBqdV － 1 2 ∫V qT BT Dε0 dV． ( 12) 由于单元不受外力作用，因此单元外力功为 0． 根据虚位移原理，可得 ∫V qT BT DBqdV = ∫V qT BT Dε0 dV． ( 13) 通过高斯数值积分可得单元刚度矩阵 ke 与单 元温度载荷向量 fe，分别为: ke = 2π ∫ 1 －1 ∫ 1 －1 BT DBrdetJdεdη; ( 14) fe = 2π ∫ 1 －1 ∫ 1 －1 BT Dε0 rdetJdεdη． ( 15) 式中，J 为单元雅可比矩阵． 利用单元节点连接关系，可分别由单元刚度矩 阵 ke和单元温度载荷向量 fe集成总体刚度矩阵 K 和总体载荷向量 F，轧辊热变形的总体弹性势能可 表示为 Π = 1 2 QT KQ － QT F． ( 16) 式中，Q 为总体位移列阵，K 为总体刚度矩阵，F 为 总体载荷向量． 对于轧辊热变形问题，对称轴上节点有径向位 移约束，中心位置节点有轴向位移约束． 采用罚函 数法处理位移约束条件，总体势能为: · 676 ·

第5期 孔繁甫等：热带钢轧机轧辊热变形的快速有限元模型 ·677· n=00+C∑i-QR. (17) 次更新温升时仅执行更新计算：计算单元载荷向量、 集成并分解总体载荷向量、回代得到热变形结果 式中，C为惩罚因子.在满足边界条件的情况下，根 由于更新计算中不涉及总刚矩阵相关的计算，计算 据最小势能原理将势能表达式对Q求最小值，可得 量将大大减小 如下平衡方程组： KO=F 换捏 更新温升：，) (18) 式中，K为约束处理后的总体刚度矩阵，F为约束处 何格划分 总体载荷向量 理后的总体载荷向量 )计算单元线荷向量 通过计算机求解上述线性方程组，可求得各节 总体刚度矩阵 2)集成总体载荷向量 1)计算积分点DB阵 3)处理位移约束 点位移，即可得到轧辊的热变形量 2计算单元刚度矩阵 )分解总体载荷向量 3.2计算过程特殊处理 集成总体刚度矩阵 4处理移约束 对于轧辊热变形计算，可认为在换辊周期内轧 闻代求解方程组 辊几何尺寸不变，且网格划分固定，所以DB矩阵是 总刚矩阵三角分解 热变形结果 固定值，结构的总体刚度矩阵也是固定值.因此，在 预算 更新计算 每次换辊时，即可完成网格划分、单元刚度矩阵计 算、总刚矩阵集成等工作. 图3轧辊热变形有限元计算程序流程 同时，根据总刚矩阵的对称、带状特性，采用半 Fig.3 Finite element calculation flow for computing the thermal ex- pansion of rolls 带宽存储以节省内存，提高访问效率；采用带状矩阵 三角分解法直接求解式(18)所示的线性方程组.由 4对比与验证 于在换辊周期内总刚矩阵是固定值，变化的仅是温 度变化而引起的载荷向量变化，所以对于式(18)所 为对比简化方法与有限元方法的计算结果，本 示的多重右边项问题，每次总刚矩阵的三角分解相 文根据上述理论模型及轧辊热变形计算流程，采用 同，只有与载荷向量相关的分解与回代过程不同. C++语言编写基于轴对称有限元方法的轧辊热变 因此，每次换辊时，在完成总刚矩阵集成后，即可进 形计算程序RoEx邓，并与现有的基于显一隐交替差 行总刚矩阵的三角分解 分的轧辊温度场计算程序RollTemp整合为一体. 根据有限元法计算流程及上述特殊处理，将计 以表1所示参数为例，计算了轧制过程中工作辊的 算任务划分为负责静态数据准备的预计算和负责动 温度场与热变形，得到热凸度如图4所示.图中 态数据准备与热变形求解的更新计算，制定如图3 Timoshenko表示采用式(4)所示的简化方法计算得 所示的轧辊热变形有限元法计算流程图.从图中可 到的热凸度，RollExp表示自编轴对称有限元计算结 以看出，在轧制单位内，预计算工作只执行一次，每 果，两种计算方法的结果差异较大 表1几何与材料特性参数 Table 1 Geometry and material parameters 参数工作辊直径，D./mm工作辊长度，L./mm弹性模量，E/GPa 泊松比4 热膨胀系数，a八05℃1 单元数量 数值 700 2000 206 0.3 1.1 100(轴向)×25（径向） 实际生产中，轧辊内部温度无法测量，所以难以 因此图中数据重合：而由简化方法计算得到的热膨 用实际热膨胀数据进行精度验证.为对比简化方法 胀量较大，比RollExp与ANSYS的计算结果大约 与自编有限元程序RollExp的精度，本文将通用有 30% 限元软件ANSYS的计算结果作为精度标准.以轧 为对比ANSYS与RollExp的求解时间，在采用 辊温度场计算程序RollTemp计算得到的20min时 主频为2.83GHz的计算平台运行了上述由2400个 和60min时的温度场数据为例（如图5所示），分别 线性四节点单元组成的轧辊热变形模型，得到 采用简化方法、RollExp和ANSYS计算轧辊在相同 ANSYS与RollExp的求解时间如表2所示.从数据 工况下的热变形，结果如图6所示.从数据对比中 中可以看出，RoEx即的计算时间远小于通用有限元 可以看出：由于计算理论相同，RollExp计算结果与 ANSYS的计算时间.RollExp在换辊时的预计算（数 ANSYS计算结果一致，在最大差别处仅差约0.5%， 据准备)时间约为40ms,温度变化时的热膨胀更新

第 5 期 孔繁甫等: 热带钢轧机轧辊热变形的快速有限元模型 Π = 1 2 QT KQ + 1 2 C ∑ q 2 i － QT F． ( 17) 式中，C 为惩罚因子． 在满足边界条件的情况下，根 据最小势能原理将势能表达式对 Q 求最小值，可得 如下平衡方程组: ^ KQ = ^ F． ( 18) 式中，^ K 为约束处理后的总体刚度矩阵，^ F 为约束处 理后的总体载荷向量． 通过计算机求解上述线性方程组，可求得各节 点位移，即可得到轧辊的热变形量． 3. 2 计算过程特殊处理 对于轧辊热变形计算，可认为在换辊周期内轧 辊几何尺寸不变，且网格划分固定，所以 DB 矩阵是 固定值，结构的总体刚度矩阵也是固定值． 因此，在 每次换辊时，即可完成网格划分、单元刚度矩阵计 算、总刚矩阵集成等工作． 同时，根据总刚矩阵的对称、带状特性，采用半 带宽存储以节省内存，提高访问效率; 采用带状矩阵 三角分解法直接求解式( 18) 所示的线性方程组． 由 于在换辊周期内总刚矩阵是固定值，变化的仅是温 度变化而引起的载荷向量变化，所以对于式( 18) 所 示的多重右边项问题，每次总刚矩阵的三角分解相 同，只有与载荷向量相关的分解与回代过程不同． 因此，每次换辊时，在完成总刚矩阵集成后，即可进 行总刚矩阵的三角分解． 根据有限元法计算流程及上述特殊处理，将计 算任务划分为负责静态数据准备的预计算和负责动 态数据准备与热变形求解的更新计算，制定如图 3 所示的轧辊热变形有限元法计算流程图． 从图中可 以看出，在轧制单位内，预计算工作只执行一次，每 次更新温升时仅执行更新计算: 计算单元载荷向量、 集成并分解总体载荷向量、回代得到热变形结果． 由于更新计算中不涉及总刚矩阵相关的计算，计算 量将大大减小． 图 3 轧辊热变形有限元计算程序流程 Fig． 3 Finite element calculation flow for computing the thermal ex￾pansion of rolls 4 对比与验证 为对比简化方法与有限元方法的计算结果，本 文根据上述理论模型及轧辊热变形计算流程，采用 C + + 语言编写基于轴对称有限元方法的轧辊热变 形计算程序 ＲollExp，并与现有的基于显--隐交替差 分的轧辊温度场计算程序 ＲollTemp 整合为一体． 以表 1 所示参数为例，计算了轧制过程中工作辊的 温度场与热变形，得到热凸度如图 4 所示． 图 中 Timoshenko 表示采用式( 4) 所示的简化方法计算得 到的热凸度，ＲollExp 表示自编轴对称有限元计算结 果，两种计算方法的结果差异较大． 表 1 几何与材料特性参数 Table 1 Geometry and material parameters 参数 工作辊直径，Dw /mm 工作辊长度，Lwr /mm 弹性模量，E /GPa 泊松比，μ 热膨胀系数，α/10 － 5℃ － 1 单元数量 数值 700 2000 206 0. 3 1. 1 100( 轴向) × 25( 径向) 实际生产中，轧辊内部温度无法测量，所以难以 用实际热膨胀数据进行精度验证． 为对比简化方法 与自编有限元程序 ＲollExp 的精度，本文将通用有 限元软件 ANSYS 的计算结果作为精度标准． 以轧 辊温度场计算程序 ＲollTemp 计算得到的 20 min 时 和 60 min 时的温度场数据为例( 如图 5 所示) ，分别 采用简化方法、ＲollExp 和 ANSYS 计算轧辊在相同 工况下的热变形，结果如图 6 所示． 从数据对比中 可以看出: 由于计算理论相同，ＲollExp 计算结果与 ANSYS 计算结果一致，在最大差别处仅差约 0. 5% ， 因此图中数据重合; 而由简化方法计算得到的热膨 胀量较大，比 ＲollExp 与 ANSYS 的计算结果大约 30% ． 为对比 ANSYS 与 ＲollExp 的求解时间，在采用 主频为 2. 83 GHz 的计算平台运行了上述由 2400 个 线性四节点单元组成的轧辊热变形模型，得 到 ANSYS与 ＲollExp 的求解时间如表 2 所示． 从数据 中可以看出，ＲollExp 的计算时间远小于通用有限元 ANSYS 的计算时间． ＲollExp 在换辊时的预计算( 数 据准备) 时间约为 40 ms，温度变化时的热膨胀更新 · 776 ·

·678 北京科技大学学报 第36卷 500 -Timoshenko 5结论 RollExp 400 (1)传统的轧辊热膨胀计算方法采用基于温升 300 分布与轴向无关的无限长柱体的解析结果，其应用 200 的假设条件与实际情况相差较远.实际生产中，轧 辊温升在轴向有较大变化，这种简化计算方法没有 250 考虑热应力对热变形的平滑作用，计算结果会高估 30 354045 轧辊的热变形量. 60 90 120 150 时间in (2)有限元法将温升引起的应变看成初始应 图4工作辊的热凸度 变，即可将热弹性问题等效为一般弹性问题，统一求 Fig.4 Thermal crown of the work roll 解.计算结果本身涵盖了热应力对热变形的平滑作 用，可计算任意给定轴对称温度场下的轧辊热变形 20 min 问题，灵活、简便且精度较高 (3)针对轧辊热变形计算的数据特性，自编轴 对称有限元程序对计算过程做特殊处理，将计算任 60 min 务分为负责静态数据准备的预计算以及负责动态数 据准备和热变形求解的更新计算，计算量大大减少， 2025 3035404550556065 计算速度快且精度高，可实现轧辊热变形的在线 温度℃ 预报. 图5工作辊温度分布 Fig.5 Temperature distribution within the work roll 参考文献 --Timoshenko(20 min)-RollExp/ANSYS(20 min) Guo R M.Two-dimensional transient thermal behavior of work rolls --Timoshenko(60 min) -RollExp/ANSYS(60 min) in rolling process.J Manuf Sci Eng,1998,120(1):28 120 2]Guo J B,Lian J C.Tu Y H.The calculation of work roll's tem- perature field and thermal crown in hot strip mills.J Yanshan 80 Unim,1998,22(3):255 0 (郭剑波，连家创，涂月红.热带钢连轧机工作辊温度场和热 凸度计算.燕山大学学报，1998,22(3)：255) B]Zhang X M,Jiang Z Y,Tieu A K,et al.Numerical modelling of -40 1000 -500 the thermal deformation of CVC roll in hot strip rolling.J Mater 0 500 1000 轴向坐标/mm Process Technol,2002,130:219 4] Shi J,Kong X W,Wang G D.Effects of work roll shifting on its 图6工作辊表面半径热膨胀 temperature field and thermal crown.Steel Rolling,2002,19 Fig.6 Thermal displacement of the work roll surface (4):7 计算时间约为3ms,求解速度满足在线应用要求. (史静，孔样伟，王国栋.工作辊横移对工作辊温度场及热凸 度的影响.轧钢，2002,19(4)：7) 表2计算时间 Zhang X L,Zhang J.Li X Y,et al.Analysis of the thermal pro- Table 2 Calculation time file of work rolls in the hot strip rolling process.JUnir Sci Technol 计算方法 时间/ms Bemg,2004,11(2):173 ANSYS 6 781.250 Du F S,Yu H,Guo Z Y.Study on temperature field of hot strip work roller.J Plast Eng,2005,12(2)78 预计算，36.450 RollExp (杜凤山，于辉，郭振宇.板带热轧工作辊温度场的研究.塑 更新计算，2.783 性工程学报，2005,12(2)：78) Lin Z C,Chen CC.Three-dimensional heat-transfer and thermal- 综上所述，自编轴对称有限元轧辊热膨胀程序 expansion analysis of the work roll during rolling.I Mater Process RollExp计算精度高（与ANSYS计算结果一致）且 Technol,2005,49(1/2):125 求解时间短，可应用于热带钢轧机板形设定模型中， 8] Guo Z F,Xu JZ,Li C S,et al.Temperature field and thermal 对轧辊热变形进行在线预报. crown of work rolls on 1700 hot strip mill.J Northeast Unir Nat Si,2008,29(4):517

北 京 科 技 大 学 学 报 第 36 卷 图 4 工作辊的热凸度 Fig． 4 Thermal crown of the work roll 图 5 工作辊温度分布 Fig． 5 Temperature distribution within the work roll 图 6 工作辊表面半径热膨胀 Fig． 6 Thermal displacement of the work roll surface 计算时间约为 3 ms，求解速度满足在线应用要求． 表 2 计算时间 Table 2 Calculation time 计算方法 时间/ms ANSYS 781. 250 ＲollExp 预计算，36. 450 更新计算，2. 783 综上所述，自编轴对称有限元轧辊热膨胀程序 ＲollExp 计算精度高( 与 ANSYS 计算结果一致) 且 求解时间短，可应用于热带钢轧机板形设定模型中， 对轧辊热变形进行在线预报． 5 结论 ( 1) 传统的轧辊热膨胀计算方法采用基于温升 分布与轴向无关的无限长柱体的解析结果，其应用 的假设条件与实际情况相差较远． 实际生产中，轧 辊温升在轴向有较大变化，这种简化计算方法没有 考虑热应力对热变形的平滑作用，计算结果会高估 轧辊的热变形量． ( 2) 有限元法将温升引起的应变看成初始应 变，即可将热弹性问题等效为一般弹性问题，统一求 解． 计算结果本身涵盖了热应力对热变形的平滑作 用，可计算任意给定轴对称温度场下的轧辊热变形 问题，灵活、简便且精度较高． ( 3) 针对轧辊热变形计算的数据特性，自编轴 对称有限元程序对计算过程做特殊处理，将计算任 务分为负责静态数据准备的预计算以及负责动态数 据准备和热变形求解的更新计算，计算量大大减少， 计算速度快且精度高，可实现轧辊热变形的在线 预报． 参 考 文 献 ［1］ Guo Ｒ M． Two-dimensional transient thermal behavior of work rolls in rolling process． J Manuf Sci Eng，1998，120( 1) : 28 ［2］ Guo J B，Lian J C，Tu Y H． The calculation of work roll’s tem￾perature field and thermal crown in hot strip mills． J Yanshan Univ，1998，22( 3) : 255 ( 郭剑波，连家创，涂月红． 热带钢连轧机工作辊温度场和热 凸度计算． 燕山大学学报，1998，22( 3) : 255) ［3］ Zhang X M，Jiang Z Y，Tieu A K，et al． Numerical modelling of the thermal deformation of CVC roll in hot strip rolling． J Mater Process Technol，2002，130: 219 ［4］ Shi J，Kong X W，Wang G D． Effects of work roll shifting on its temperature field and thermal crown． Steel Ｒolling，2002，19 ( 4) : 7 ( 史静，孔祥伟，王国栋． 工作辊横移对工作辊温度场及热凸 度的影响． 轧钢，2002，19( 4) : 7) ［5］ Zhang X L，Zhang J，Li X Y，et al． Analysis of the thermal pro￾file of work rolls in the hot strip rolling process． J Univ Sci Technol Beijing，2004，11( 2) : 173 ［6］ Du F S，Yu H，Guo Z Y． Study on temperature field of hot strip work roller． J Plast Eng，2005，12( 2) : 78 ( 杜凤山，于辉，郭振宇． 板带热轧工作辊温度场的研究． 塑 性工程学报，2005，12( 2) : 78) ［7］ Lin Z C，Chen C C． Three-dimensional heat-transfer and thermal￾expansion analysis of the work roll during rolling． J Mater Process Technol，2005，49( 1 /2) : 125 ［8］ Guo Z F，Xu J Z ，Li C S，et al． Temperature field and thermal crown of work rolls on 1700 hot strip mill． J Northeast Univ Nat Sci，2008，29( 4) : 517 · 876 ·