第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 ·