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