板带轧制过程考虑滑动与粘着摩擦的轧制力模型

李维刚，鲁凌云, 刘超，刘相华



板带轧制过程考虑滑动与粘着摩擦的轧制力模型

李维刚1，鲁凌云1, 刘超1，刘相华2

(1. 武汉科技大学信息科学与工程学院，湖北武汉，430081；2. 东北大学研究院，辽宁沈阳，110819)

将轧制变形区分为滑动摩擦区和粘着摩擦区，得到适用于带钢热轧过程的改进卡尔曼(KARMAN)微分方程，给出轧制力数值计算公式。结合变形区入口与出口的边界条件，采用龙格−库塔(Runge−Kutta)法交替求解前滑区与后滑区的单位压力分布，2个接触弧角相交处就是中性点，通过对单位压力与摩擦应力进行数值积分得到轧件的轧制力，再根据希齐柯克(Hitchcock)公式迭代计算多次得到最终轧制力。最后，分析新型轧制力模型应用于热连轧机组面临的问题，给出热轧过程摩擦因数对单位压力分布与轧制力的影响规律，并通过布伦特(Brent)方法实现对热轧过程摩擦因数的软测量，指出摩擦因数在1个轧制周期内的演变规律，为热连轧过程摩擦因素在线模型的建立与新型数值轧制力模型的在线应用提供参考。

热轧带钢；轧制力；微分方程；摩擦因数

轧制力是板带钢轧制过程控制的关键参数，其计算精度直接影响着整个轧制过程，是板形、板厚设定控制的基础[1−2]，国内外学者对轧制力计算进行了大量研究[3−11]。研究单位压力在接触弧上的分布规律，对准确计算板带钢轧制时的轧制力具有重要意义。计算单位压力经常采用Karman(卡尔曼)微分方程或Orowan(奥洛万)微分方程，两者均建立了轧制变形区内微元体上的力平衡关系，很多轧制力解析公式是以这2个方程为基础推导出来的。如通过求解Karman方程得到Tselikov(采利柯夫)公式、Stone(斯通)公式，通过求解Orowan方程得到Sims(西姆斯)公式、Bland(勃兰特)公式。其推导过程分别采用不同的假设条件或化简条件，对轧制力计算精度必然造成一定影响。随着数值计算方法的兴起与计算机性能的提高，通过数值方法在线求解单位压力微分方程成为可能。与目前广泛在线使用的解析式轧制力模型不同，此类模型为数值轧制力模型。由于其假设条件较少，模型原理及结构清晰，开放性强，为精细化计算轧制力提供了新方法。从理论上说，只要摩擦因数与流变应力准确，数值轧制力模型可以达到所希望的预测精度。由代数模型向微分模型发展，由宏观模型向微观模型发展，代表了轧制模型未来的发展方向，因此，开发数值轧制力在线计算模型很有意义。本文研究一种滑动摩擦与粘着摩擦共存条件下的数值轧制力模型，给出适合于在线计算的数值解法，探讨新型轧制力模型用于热连轧机组在线计算面临的问题，指出摩擦因数对单位压力分布与轧制力的影响规律，最后通过Brent方法实现对热轧过程摩擦因数的软测量。

1 滑动/粘着摩擦共存的轧制力模型

通过对变形区内轧件微元体所受应力的分析，Karman推导了轧制过程的力平衡微分方程，即Karman方程[3]。在其推导过程中，认为整个接触区都是滑动摩擦。当摩擦力不大时，Karman方程近似程度较好，被广泛用于冷轧轧制力计算。当轧制过程表面摩擦力等于或接近材料剪切屈服极限时，方程采用的基本假设偏差大，故不宜直接用于热轧过程。CHEN等[3]同时考虑滑动摩擦与粘着摩擦，得到了适用于热轧过程的力平衡微分方程。

1.1 改进Karman微分方程

在轧制变形区内任取1个微元，分析作用在微元上的作用力(见图1)，对微元进行受力分析得

对微元进行水平受力分析，由力平衡得

式中：p为单位压力；p′为单位压力的垂直分量；τ为摩擦应力；θ为任意微元的接触角；h为任意微元的高度；hin和hout分别为轧件入口厚度与出口厚度；σin和σout分别为入口张力与出口张力。式(1)和(2)构成了微元在垂直方向与水平方向上的力平衡方程。

图1 变形区内微元体的应力平衡关系

根据轧制过程中的板厚变化，由几何关系得：

式中：′为轧辊压扁半径。在轧制过程中，使用剪切变形功理论作为材料的屈服条件，有

式中：2为流变应力。轧件与轧辊间的接触摩擦非常复杂，对此所提出的假设或理论很多，如认为遵循干摩擦定律、假设接触表面上单位摩擦不变、采用液体摩擦定律、应用经验公式等。本文将变形区分为滑动摩擦区和粘着摩擦区，分别应用不同摩擦规律，推导得到改进Karman微分方程。

1.1.1 滑动摩擦区

根据滑动摩擦假设，认为在变形区内轧件与轧辊的接触表面产生滑移，遵循干摩擦定律，其接触表面的摩擦应力为

式中：为摩擦因数，正号表示后滑区(入口侧)，负号表示前滑区(出口侧)。根据式(1)~(6)，可推导得到滑动摩擦区内的微分方程为

式中：

式(8)和(9)中的分母取正号表示后滑区，取负号表示前滑区。

1.1.2 粘着摩擦区

根据粘着摩擦假设条件，认为在变形区内轧件与轧辊的接触表面产生粘着，可用剪切屈服极限代替接触表面的摩擦应力：

=±;∙≥(10)

根据式(1)~(5)及(10)，可推导得到粘着摩擦区内的微分方程：

式中：

式(12)取负号代表后滑区，取正号代表前滑区。式(7)和(11)是结合力平衡方程、几何条件、塑性方程与摩擦规律推导出来的结果，即改进Karman方程。由于在轧制变形区内同时考虑了滑动摩擦与粘着摩擦，故改进Karman方程可用于分析热轧过程，克服了原Karman方程在整个变形区仅考虑滑动摩擦的 缺陷。

1.2 轧制力计算公式

对轧件进行轧制力分析，假设接触弧变形后仍为圆弧，在接触弧上作用有垂直于辊面切线的单位压力和平行于辊面切线的摩擦应力。轧制力由单位压力和摩擦应力的垂直分量合成：

1.3 数值求解

在式(13)所示的轧制力计算公式中，需要事先获得单位压力与摩擦应力沿轧辊接触弧的分布，因此，要先求解微分方程(7)与(11)。

求解微分方程(7)和(11)，要先确定微分方程的边界条件。在接触弧的入口和出口平面，边界条件可通过微元的平衡方程予以确定。入口处的单位压力in如下。

∙in≤in时的摩擦为滑动摩擦，

∙in＞in时的摩擦为粘着摩擦，

出口处的单位压力out为

入口处、出口处的接触角分别为：

式(7)和(11)为带边界条件的微分方程，要得到其解析解非常困难，可通过数值方法进行求解。龙格−库塔方法(Runge−Kutta，简称为R−K法)是一种在工程上应用广泛的高精度单步算法，其理论基础来源于泰勒展开公式和使用斜率近似表达微分，它在积分区间多预计算几个点的斜率，然后进行加权平均，用于计算下一点低斜率的依据，从而构造出精度更高的数值积分计算方法。本文采用5阶R−K方法求解改进Karman微分方程，得到单位压力沿轧辊接触弧的分布。

在R−K法求解过程中，将轧件与轧辊的入口接触弧角in分为等分，由此变形区内得到个微分单元，如图2所示。结合变形区入口的边界条件，逐步求出后滑区每个微分单元的单位压力，接触弧角由in逐步减小；结合变形区出口的边界条件，逐步求出前滑区每个微分单元的单位压力，接触弧角由0°逐步增大；两者相交处(单位压力相等)就是中性点，此时的接触弧角就是中性角N。前滑区与后滑区的求解过程可以交替进行，如设置当前滑区求得的单位压力比后滑区大时则启动后滑区计算，后滑区求得的单位压力比前滑区大时则启动前滑区计算。

图2 变形区内微分单元划分

采用R−K法求得轧辊接触弧上的单位压力分布，根据式(6)和(10)可求摩擦应力分布，根据式(13)可计算轧件单位宽度上的轧制力。需要注意的是：微分方程(7)和(11)的右式中均含有轧辊压扁半径′，而初始的压扁半径并不知道，初次计算可取原始轧辊半径。根据本文方法求得轧件的轧制力后，再根据希齐柯克(Hitchcock)公式更新计算轧辊压扁半径。以上过程经过多次迭代，当2次计算得到的轧制力偏差小于预先给定的正值小量时，认为满足收敛条件，结束 迭代计算。

2 新模型在热连轧中的应用

这里探讨新型数值轧制力模型应用于热连轧机组在线计算面临的问题。图3所示为新型轧制力模型的结构图。由图3可知：要采用新型轧制力模型在线计算热轧带钢的轧制力，需要事先建立流变应力和摩擦因数2个子模型；根据子模型计算出流变应力、摩擦因数后，采用本文给出的方法可求得变形区内单位压力及摩擦应力的分布，对单位压力与摩擦应力进行数值积分求得轧制力，再根据Hitchcock公式更新轧辊压扁半径，经过几次迭代计算可得到最终轧制力。

图3 数值型轧制力模型的结构

目前，国内外关于流变应力的研究很多，其在线计算模型也比较成熟，但对热轧过程摩擦因数的研究并不多见。以往研究多针对特定条件来设计实验或实验装置来测定摩擦因数[12−14]，能否直接用于热连轧在线计算有待进一步考证。研究表明：轧制时变形区内摩擦因数对整个轧制过程的影响关系较复杂，几乎所有工艺参数(轧制温度、轧材性质、压下制度、润滑剂等)和设备参数(辊面粗糙度、轧制速度等)的变化都会影响变形区内的摩擦状态，从而引起摩擦因数的变化；而摩擦因数的变化将导致轧制过程的力能参数发生变化。因此，有必要研究摩擦因数对单位压力与轧制力的影响。下面通过数值模拟来予以说明。

2.1 摩擦因数对单位压力分布的影响

选取典型热轧带钢F4机架为仿真对象，轧辊直径为824.90 mm，带钢宽度为1 543.50 mm，入口厚度为10.85 mm，出口厚度为8.19 mm，入口张力为 8.28 MPa，出口张力为9.77 MPa；轧件材质为普碳钢，轧制温度为933 ℃，流变应力为207.00 MPa，摩擦因数变化范围为0.15~0.50。图4所示为单位压力、摩擦应力随摩擦因数的变化曲线。从图4可见：

1) 单位压力、中性角随着摩擦因数的增大而增大，当摩擦因数增大到0.30以上时，单位压力的增加幅度减缓。这是因为变形区开始出现粘着摩擦，整个变形区逐步由滑动摩擦过渡到粘着摩擦。

2) 摩擦应力随着摩擦因数的增大而增大，当摩擦因数增加到0.40以上时，摩擦应力几乎不再增大。这是因为当摩擦因数较小时变形区内发生的是滑动摩擦，摩擦应力随摩擦因数呈线性增大；当摩擦因数较大时变形区内发生粘着摩擦，摩擦应力取材料剪应屈服极限(本例=103.5 MPa)。

摩擦因数：1—0.15；2—0.20；3—0.25；4—0.30；5—0.35；6—0.40；7—0.45；8—0.50。

图4 摩擦因数对单位压力和摩擦应力分布的影响

Fig. 4 Effect of friction coefficient on distribution of normal pressure and friction shear stress

图5所示为轧制力与应力状态系数p随着摩擦因数的变化。从图5可见：轧制力与应力状态系数p开始都随着摩擦因数的增大而增大，当摩擦因数增加到一定程度时，两者的增加效果都不再明显。这是因为轧制变形区内的摩擦发生变化，逐步由滑动摩擦过渡到粘着摩擦。本仿真实例中，当摩擦因数从0.15增加到0.50时，轧制力从12.7 MN增大到19.6 MN，变化非常显著。因此，若能适当采取措施提高轧制区的润滑效果，减小摩擦因数，则可以有效降低精轧带钢的轧制力。

1—轧制力；2—p。

图5 摩擦因数对轧制力和应力状态系数的影响

Fig. 5 Effect of friction coefficient on roll force and stress state coefficient

图5中的应力状态系数p由下式求出：

其中：为根据式(13)求得的单位宽度的轧制力；d为轧辊的压扁弧长。

2.2 热轧过程摩擦因数的软测量

上述分析表明：热轧过程摩擦因数对单位压力分布、轧制力的影响显著，有必要研究热连轧过程摩擦因数的变化规律。由于热轧过程摩擦因数的直接测量非常困难，寻找间接的软测量方法变得很有意义。本文给出了在已知摩擦因数与流变应力的情况下求解轧制力的方法，即建立了如下函数映射关系：

因此，若已知轧制力与流变应力，则可以采用基于非导数优化的方法反求出对应的摩擦因数。布伦特(Brent)方法[15]是一种高效可靠的求解有约束单变量极值问题的方法，它巧妙结合了黄金分割搜索和抛物线插值两者的优点：基于抛物线插值的思想，在优化迭代步长的选取方面，将抛物线插值与黄金分割法相结合，从而自适应限定移动步长并确定新的区间包含最优解。选择Brent方法的优点在于：无需提供导数信息，可同时考虑输入约束，具有很好的自适应性。

选取某热连轧生产线连续4个轧制计划的实际生产数据(钢种为中低碳钢，宽度为800~1 450 mm，厚度为2.0~6.0 mm；轧辊材质，前4机架为高速钢辊，后3机架为高Gr铸铁轧辊)，将数值轧制力模型计算方法与Brent方法相结合，求得热轧过程的摩擦因数的变化曲线如图6所示。

由图6可知：1) 在轧制过程中，各机架的摩擦因数均在变化，相对而言，F1~F5机架变化幅度较小，而F6~F7机架变化幅度较大；2) 在每个轧制计划初期，由于轧制长度较少、轧辊辊面较光滑，故摩擦因数稍偏小；随着轧制长度延长，摩擦因数有所增加但增加幅度不大，到达一定长度后，除F7机架外，其余机架几乎保持恒定。

F7机架摩擦因数的变化可能与其压下率较小、摩擦因数对轧制条件较敏感有关，具体原因需要深入 分析。

与实验条件下对热轧过程摩擦因数的测量不同，通过数值轧制力模型与Brent方法实现对热轧生产过程摩擦因数的软测量很有意义，它为热轧摩擦因数在线模型的建立提供了数据基础。收集对摩擦因数影响较大的因素如轧辊速度、轧制温度、压下率、轧制长度等参数，通过回归方法即可得到热轧摩擦因数的在线计算模型，这为数值轧制力模型在热连轧机组中的在线应用提供了参考。

机架：1—F1；2—F2；3—F3；4—F4；5—F5；6—F6；7—F7。

图6 热连轧精轧过程中摩擦因数的演变

Fig. 6 Evotution of friction coefficient during finishing rolling process in hot strip mills

3 结论

1) 对板带轧制过程，通过区分滑动摩擦与粘着摩擦，结合应力平衡关系、摩擦规律、几何条件与塑性方程等条件，得到了滑动摩擦与粘着摩擦共存条件下的改进Karman微分方程，该方程可用于热轧过程的分析计算。根据方程求得的接触弧上的单位压力与摩擦应力分布，给出了轧制力数值计算公式，提出了数值型轧制力模型的体系结构。

2) 给出改进Karman微分方程的数值解法，结合变形区入口与出口的边界条件，采用五阶龙格−库塔法交替求解前滑区与后滑区的单位压力分布，2个接触弧角相交处(单位压力相等)就是中性点；得到接触弧上的单位压力与摩擦应力分布后，通过对其进行数值积分计算轧件的轧制力，再根据Hitchcock公式更新轧辊压扁半径，经过多次迭代，当2次计算得到的轧制力偏差足够小时，认为满足收敛条件，得到最终的轧制力。

3) 热连轧精轧过程轧制力随摩擦因数的增大而增大；当摩擦因数增加到一定程度时，轧制力的增加效果不明显。借助数值轧制力模型与Brent方法，实现了对热连轧过程摩擦因数的软测量，指出摩擦因数在1个轧制周期内的演变规律，为热轧过程摩擦因数在线模型的建立与数值型轧制力模型在线应用提供了参考。

[1] 刘相华, 胡贤磊. 轧制参数计算模型及其应用[M]. 北京: 化学工业出版社, 2007: 36−46. LIU Xianghua, HU Xianlei. Rolling parameter calculation model and its application[M]. Beijing: Chemical Industry Press, 2007: 36−46.

[2] 曹建国, 张杰, 张少军. 轧钢设备及自动控制[M]. 北京: 化学工业出版社, 2010: 64−74. CAO Jianguo, ZHANG Jie, ZHANG Shaojun. Rolling equipment and automatic control[M]. Beijing: Chemical Industry Press, 2010: 64−74.

[3] CHEN Shuixuan, LI Weigang, LIU Xianghua. Calculation of rolling pressure distribution and force based on improved Karman equation for hot strip mill[J]. International Journal of Mechanical Sciences, 2014, 89: 256−263.

[4] 李维刚, 谭树彬, 李家波, 等. 热连轧机轧制力和轧制力矩模型研究[J]. 东北大学学报(自然科学版), 2011, 32(5): 622−625.LI Weigang, TAN Shubin, LI Jiabo, et al. Modelling rolling force and rolling torque in a hot strip mill[J]. Journal of Northeastern University(Natural Science), 2011, 32(5): 622−625.

[5] 李维刚, 陈水宣, 刘相华. 热轧带钢精轧过程考虑相变的轧制力模型研究[J]. 东北大学学报, 2013, 34(10): 1425−1429. LI Weigang, CHEN Shuixuan, LIU Xianghua. Rolling force model considering phase transformation for hot strips during finishing process[J]. Journal of Northeastern University (Natural Science), 2013, 34(10): 1425−1429.

[6] GUO Mingyuan, JIAN Wang, HONG Xiao, et al.Research on online model of vertical rolling force in hot strip roughing trains[J]. Advanced Materials Research, 2010, 145(5): 198−203.

[7] BAGHERIPOOR M, BISADI H. Application of artificial neural networks for the prediction of roll force and roll torque in hot strip rolling process[J]. Applied Mathematical Modelling, 2013, 37(7): 4593−4607.

[8] 周富强, 曹建国, 张杰, 等. 基于神经网络的冷连轧机轧制力预报模型[J]. 中南大学学报(自然科学版), 2006, 37(6): 1155−1160. ZHOU Fuqiang, CAO Jianguo, ZHANG Jie, et al. Prediction model of rolling force for tandem cold rolling mill based on neural networks and mathematical models[J]. Journal of Central South University (Science and Technology), 2006, 37(6): 1155−1160.

[9] 王健, 刘长青, 肖宏. 热连轧精轧机组预设定和自适应研究[J]. 中南大学学报(自然科学版), 2012, 43(7): 2607−2612. WANG Jian, LIU Changqing, XIAO Hong. Research on finishing setup and self-adaptive of hot strip mill[J]. Journal of Central South University (Science and Technology), 2012, 43(7): 2607−2612.

[10] 陈金山, 李长生, 曹勇. 不锈钢冷连轧轧制力显函数模型的开发与应用[J]. 中国冶金, 2015, 25(3): 6−12. CHEN Jinshan, LI Changsheng, CAO Yong. Development and application of explicit function model on rolling force for stainless steel in tandem cold rolling[J]. China Metallurgy, 2015, 25(3): 6−12.

[11] 杨景明, 顾佳琪, 闫晓莹, 等. 基于改进遗传算法优化BP网络的轧制力预测研究[J]. 矿冶工程, 2015,35(1): 111−115. YANG Jingming, GU Jiaqi, YAN Xiaoying, et al. Rolling force prediction based on BP network optimized by an improved genetic algorithm[J]. Mining and Metallurgical Engineering, 2015, 35(1): 111−115.

[12] 顾正秋, 于宗源, 蔡英翘. 用制动仪测定热轧摩擦系数[J]. 钢铁, 1982, 18(10): 31−34. GU Zhengqiu, YU Zhongyuan, CAI Yingqiao. Determination of the coefficient of friction in hot rolling by a new tribometer[J]. Iron and Steel, 1982, 18(10): 31−34.

[13] NILSSON M, OLSSON M. Microstructural, mechanical and tribological characterization of roll materials for the finishing stands of the hot strip mill for steel rolling[J]. Wear, 2013, 307(1/2): 209−217.

[14] CHENG Xiawei, JIANG Zhengyi, WEI Dongbin, et al. Oxide scale characterization of ferritic stainless steel and its deformation and friction in hot rolling[J]. Tribology International, 2015, 84(4): 61−70.

[15] BRENT R P. Algorithms for minimization without derivatives[M]. NJ: Prentice-Hall, 1973: 200−248.

(编辑 陈灿华)

Roll force model in consideration of slipping and sticking friction during strip rolling process

LI Weigang1, LU Lingyun1, LIU Chao1, LIU Xianghua2

(1. College of Information Science and Engineering, Wuhan University of Science and Technology,Wuhan 430081, China;2. Research Institute of Science and Technology, Northeastern University, Shenyang 110819, China)

An improved KARMAN equation for hot rolling process was deduced in consideration of slipping and sticking friction in the roll-bite, and the numerical formula of roll force calculation was given. Combined with boundary condition of inlet and outlet in the roll bite, the Runge−Kutta method was adopted to alternately solve the distribution of normal pressure in the forward and backward slip zones. The intersection of contact angle was just the neutral point. The roll force could be calculated by integrating the normal pressure and friction shear stress in the roll bite, and then the final roll force was obtained by repeated iteration calculation in the use of Hitchcock formula. Finally, the problem in the application of the numerical roll force model to hot strip mill was analyzed, and the influence of friction coefficient on the normal pressure distribution and roll force during hot rolling was shown. The soft measurement of friction coefficient during hot rolling process was carried out by Brent method, and evolution rule of friction coefficient in a rolling recycle was pointed out, which provides references for the establishment of online mathematical model of friction coefficient and the application of the proposed numerical roll force model in online process control of hot strip mill.

hot-rolled strip; roll force; differential equation; friction coefficient

10.11817/j.issn.1672-7207.2016.12.009

TG335.5

A

1672−7207(2016)12−4021−07

2015−12−25；

2016−02−06

湖北省教育厅科学技术研究计划重点项目(D20161103)；武汉市青年科技晨光计划项目(2016070204010099) (Project(D20161103) supported by the Ministry of Education Science and Technology Research Program of Hubei Province; Project (2016070204010099) supported by the Youth Science and Technology Program of Wuhan Minicipality)

李维刚，博士，教授，从事冶金工业过程控制与智能系统研究；E-mail：liweigang.luck@foxmail.com