一种能够精确计算冷轧薄硬带材滚动力的方法，并对工作辊的弹性压平进行了优化

《Applied Mathematical Modelling》：Method to accurately calculate rolling force of cold-rolled thin-hard strip with improvement for elastic flattening of work roll

　　•提出了一种用于表征工作辊弹性压扁现象的改进方法。•能够准确计算使用大直径辊子轧制薄硬带材时的轧制力。•该方法弥补了传统Hitchcock圆形轮廓假设的不足。

引言
轧制力计算是冷轧过程中变形量分配的基础[1,2]，如果计算值不准确，轧机中的计算机系统就会给出不合理的变形量分配策略，从而导致实际轧制过程中的轧制负荷（轧制力、轧制扭矩、轧制功率）出现异常[3]，甚至超过轧机的承载能力。当轧机处于超负荷状态时，生产参数的微小变化都可能导致轧制负荷大幅波动，这会严重影响轧制过程的稳定性，进而使得带材的冷轧变得十分困难[4,5]。为避免这一问题，必须准确计算冷轧过程中的轧制力[6,7]。

轧制力受多种因素影响[8,9]，工作辊与带材之间的接触状况是一个重要方面[10,11]。当工作辊受到较大的轧制力作用时，由于弹性变形，它会在压力作用下变“平”，这种现象被称为弹性压扁[12,13]。当高强度钢带（如双相钢和高等级硅钢）在六架高速轧机中轧制时，由于工作辊直径较大（约380毫米至430毫米），且带材的屈服强度很高（超过1000兆帕）[14,15]，工作辊的弹性压扁现象会十分明显，尤其是在多机架冷轧过程中，随着轧制的进行，带材逐渐变薄，加工硬化程度也越来越高，因此工作辊的弹性压扁现象会更加显著。所以，在工作辊直径较大且轧制力较大的情况下（由带材的薄硬特性导致），在轧制力计算中必须准确考虑工作辊的弹性压扁现象。如果弹性压扁的设定不合理，轧制力的计算结果就会不准确，与实际工业生产情况可能存在较大偏差[16,17]。

传统的计算方法是Hitchcock方程[18]，该方程假设压扁后的工作辊轮廓仍为半径增大的圆形，如图1所示，RD表示轧制方向，工作辊与带材最初的接触曲线为AB?，压扁后的接触曲线为AB*?，相应的压扁半径可通过公式(1)计算得出。轧制力的计算步骤如下：首先根据工作辊的原始半径，通过轧制力模型计算出轧制力P，然后通过公式(1)计算出压扁半径，再通过公式(2)计算出压扁后的变形区长度，最后将上述数值代入轧制力模型中，计算出压扁后的轧制力PF。许多研究人员采用这种方法来计算冷轧过程中的轧制力，Tang等人[19]基于轧制力方程和Hitchcock方程，提出了不对称冷轧过程中最小厚度的隐式表达式；Jacobs等人[11,20]则量化了冷轧过程中各向异性材料的屈服行为及其对轧制力的影响，他们的轧制力计算也采用了Hitchcock方程。然而，在不同的带材厚度和变形抗力条件下，这种方法的轧制力计算精度会有所差异[[21], [22], [23]]，因此Hitchcock圆形假设的适用性可能有限。

RF=R·[1+16(1?ν2)πE·PBs·Δh]
LF=(RF·Δh)

从工作辊弹性压扁对轧制力的影响机制来看，如图1所示，一方面，工作辊的弹性压扁会增大变形区的长度，进而增加轧制力；另一方面，压扁后接触曲线的形状发生变化，这会改变中性点的位置以及变形区内塑性区与弹性区的比例，进而改变轧制压力分布，最终影响轧制力。在使用传统的Hitchcock方程计算薄硬带材的轧制力时，虽然变形区的长度已经得到准确计算，但在超高强度钢的冷轧过程中，发现压扁后的接触曲线并非圆形[[24], [25], [26]]，实际上，如图2所示，它可能变为曲线2（介于原始曲线和圆形压扁曲线1之间），或者变为曲线3（比曲线1更平缓）。在这种情况下，轧制压力分布会发生变化，从而影响总的轧制力，因此使用Hitchcock圆形假设来计算薄硬带材的轧制力可能会出现较大偏差。一些学者持这一观点，并根据对变形区划分或辊子变形的不同处理方式，提出了新的弹性压扁计算模型，如Fleck模型[[27], [28], [29]]和Nakajima-Matsumoto模型[30]。因此，在薄硬带材的冷轧过程中，应精确处理接触曲线的形状。

本文针对冷轧薄硬带材，提出了一种用于计算工作辊弹性压扁的改进方法。具体而言，该方法改变了轧制变形区内工作辊的轮廓曲线，结合实际工业轧制条件，分别用该方法与传统方法计算不同薄硬带材的轧制力，然后对比分析两种方法的差异。同时，通过与实际测量的轧制力进行比较，验证了改进方法在计算轧制力方面的准确性。

本文的结构如下：第二部分以不同的薄硬带材为研究对象，评估在不同条件下传统辊子压扁方法计算轧制力的准确性；第三部分基于分析结果，提出一种用于计算工作辊弹性压扁的改进方法，以提高轧制力计算的准确性，并分析传统方法、改进方法与实际测量值之间的差异；第四部分总结主要结论。

工业背景为：某工厂的五架高速六辊冷轧机用于生产一种具有薄硬特性的钢带，该轧机的工作辊直径较大（385毫米至425毫米），成品带材的厚度较小（0.30毫米至0.50毫米），且在后续机架中带材的变形抗力很大，超过了1000兆帕。由于不同机架中带材的厚度和变形抗力存在差异，传统方法在计算轧制力时的适用性也有所不同，因此首先需要评估传统方法在计算不同薄硬带材轧制力时的准确性，这为改进弹性压扁计算提供依据。

以一个典型规格为例，该带材的宽度为1286毫米，进料厚度为2.598毫米，这也是S1机架入口处的厚度，成品厚度为0.347毫米，各机架出口处的带材厚度见表1，各机架分别简写为S1、S2、S3、S4和S5，各机架的带材变形抗力、工作辊半径以及实测轧制力也列于表1中，其中变形抗力为轧机计算机中变形抗力模型的计算值，半径为轧制前的工作辊实测值。

如引言中所述，在考虑工作辊弹性压扁的轧制力计算中，应首先根据工作辊的原始半径计算出初始轧制力P，因此本节首先建立基本计算模型，然后根据工作辊压扁后的Hitchcock圆形假设，计算出压扁后的轧制力PF。

为深入理解带材与工作辊之间接触曲线变化对轧制力影响的机制，轧制力计算模型采用了带区域划分的轧制压力分布模型，即根据带材的变形特征以及带材与工作辊在不同位置的速度差，将变形区划分为多个不同的区域，分别求解每个区域的轧制压力，再将各区域的压力积分后求和，从而得到轧制力。实际上，在冷轧过程中，变形区内同时存在弹性区和塑性区，带材先发生弹性变形，随后才发生塑性变形，而且在变形区出口处还存在弹性恢复现象。此外，以中性点为分界点，变形区还可以划分为前滑移区和后滑移区。因此，如图3所示，从变形区入口到出口，依次划分为4个部分：入口弹性区（区域a）、塑性-后滑移区（区域b）、塑性-前滑移区（区域c）以及出口弹性区（区域d），分别对应入口处的弹性变形、伴有后滑移的塑性变形、伴有前滑移的塑性变形以及出口处的弹性变形，其中区域a属于后滑移区，区域d属于前滑移区。轧制压力模型采用了Karman微分方程[31]和Lian公式[32,33]，其中Karman微分方程用于计算塑性区的轧制压力，Lian公式用于计算弹性区的轧制压力。然后将每个区域的轧制压力在其对应区域内积分后再求和，再乘以带材宽度，即可得到总轧制力，如公式(3)所示，积分的下限为每个区域起点的x坐标，上限为该区域的终点x坐标。

P=BS·(∫0Δx0paxdx+∫Δx0xcpbxdx+∫xcLpcxdx+∫LL+Δx1pdxdx)

此外，在图3中，点A为变形区入口，点C为入口弹性区与塑性-后滑移区的边界点，点D（中性点）为塑性-后滑移区与塑性-前滑移区的边界点。由于带材的弹性恢复作用，点B（即图1中与工作辊最低点接触的点）为变形区的理论出口位置，而点E为变形区的实际出口位置，可以看出实际出口厚度he略大于理论出口厚度h1。点B同时也是塑性-前滑移区与出口弹性区的边界点。

如图3所示，以半厚度带材为研究对象，以带材厚度对称线为x轴，带材厚度方向为y轴，点A位于y轴上，以此建立用于计算的坐标系。工作辊的轮廓曲线方程为公式(4)，在变形区内，由于带材的厚度分布曲线与工作辊的轮廓曲线完全重合，其方程为公式(7)。

Wx=h12+R?R2?(L?x)2
L=Δh·R
Δh=h0?h1
hx=2·Wx

计算过程包括三个步骤：确定各区域的边界及中性点位置、计算各区域的轧制压力、计算总轧制力，详细的计算过程已在之前的研究[34]中阐述。在确定各区域边界时，首先需要计算入口弹性区和出口弹性区的长度，如公式(8)、(9)所示。由于Lian公式的特点，其解需要与轧制压力分布的计算相结合，而中性点xc的位置则通过金属流动速度分布的解来确定。区域A和区域D的轧制压力是通过Lian公式计算的，分别如式(10)、(11)所示；区域B和区域C的轧制压力则是通过Karman微分方程计算的，对应式(12)、(13)。中间物理量的计算公式见附录。Δx0=L2+ξ02L21?ν1ψ0?LΔx1=ξ1L1?ν1ψpax=pC·B0x?(1?B0L?B0Δx0)·[exp(B0x)?1]1?B0L?(1?B0L?B0Δx0)·exp(B0Δx0)pdx=pB·Bx?(1?BΔx1)·[exp(Bx)?1]1?(1?BΔx1)·exp(BΔx1)dpbxdx+Khx·dhxdx?2·μpbxhx=0dpcxdx+Khx·dhxdx+2·μpcxhx=0在计算变形区的轧制压力分布时，将点A处的轧制压力设为0MPa，入口弹性区的轧制压力通过式(10)计算得出点C处的轧制压力pC，然后在塑性后滑移区，轧制压力通过式(12)计算，边界条件为px?=?pC。式(12)是一个带初始值的一阶微分方程，其一般形式如式(14)所示，该轧制压力可通过欧拉数值算法求解，该算法将微分方程转化为差分方程以便计算，从而也可得到中性点处的轧制压力pxc，该值作为式(13)的初始值，塑性前滑移区的轧制压力同样通过欧拉算法求解，得到点B处的轧制压力pB，出口弹性区的轧制压力则通过式(11)计算。使用欧拉法时，区域B和区域C中的节点数量取决于区域长度，且各区域的节点间距相等。最后，根据式(3)计算轧制力。{dpdx=f(x,p)p(x0)=p0该基本模型用于计算原始轧制力，将其代入式(1)可计算出压下半径，再代入式(4)确定压下后工作辊的轮廓曲线。依据基本模型的计算流程，还可计算出压下后的轧制力，此时轧制力和压下半径需要通过迭代方式求解。S3、S4、S5三种厚度（从大到小）和变形抗力（从小到大）的带材的轧制参数被设定为计算条件，如表1、表2所示。其中摩擦系数为轧机计算机中摩擦系数模型的计算值，前后张力则由张力表的测量值转换而来，轧制速度为激光测速仪的测量值。进而计算不同机架对应的轧制力以进行比较。不同薄硬带材（S3、S4、S5）的轧制压力分布计算结果如图4所示，与原始工作辊的对比可知，在S3和S4中，考虑到工作辊的压下情况，轧制压力分布的变化较小，变形区长度的增加也很少；而在S5中，考虑到工作辊的压下后，轧制压力分布发生显著变化，塑性后滑移区和塑性前滑移区的轧制压力以及变形区长度都有明显改变，变形区总长度从3.47毫米增加到4.94毫米，增幅为原始长度的42.36%。同时，还计算了S3、S4、S5中工作辊的压下半径，并与原始半径进行比较，如图5所示，S3和S4中两者的半径差异较小，相对原始半径的增幅分别为5.47%和8.30%；而S5的压下半径为402.173毫米，约为原始半径的两倍。因此，在薄硬带材轧制过程中，工作辊经过弹性压下后半径会大幅增加，这也是变形区长度和轧制压力分布显著增大的原因。不同薄硬带材的轧制力也进行了计算，如图6所示，S3、S4、S5的计算值分别为13584千牛、11570千牛和9613千牛，随后与测量值对比以验证准确性，可见对于相对较厚较软的S3、S4带材，计算值较为准确，计算值与测量值的差异分别为8千牛和160千牛，根据式(15)计算得的相对偏差分别为0.06%和1.40%；而对于最薄最硬的S5带材，计算值不够准确，二者差异达1317千牛，相对偏差为12.05%。通过上述对比分析可知，对于最薄最硬的S5带材，工作辊的弹性压下程度较大，考虑弹性压下后轧制力的变化十分显著；而相对较厚较软的S3、S4带材这种情况并不明显。因此，在工作辊弹性压下程度较大且轧制力变化显著的条件下，应精确处理变形区内工作辊的轮廓曲线，这样才能进一步提升计算精度，否则计算得到的轧制力可能与测量值存在一定偏差。EP=|PF?Pmea|Pmea×100%段落节选工作辊弹性压下的改进计算方法如2.2节所述，应精确考虑变形区内压下后工作辊的轮廓曲线，一些研究者发现将曲线设为二次函数可更准确地计算冷轧超高强度带材的轧制力[[24], [25], [26]]。在此，我们尝试将曲线设为二次函数，并分析其对轧制力的影响，同时提出了确定该曲线的改进方法，曲线方程如式(16)所示。结论本研究针对采用大直径工作辊的6辊轧机轧制具有薄硬特性的带材，由于传统的Hitchcock圆形假设（传统方法）可能在轧制力计算上不够准确，因此提出了一种改进方法，用以精确描述变形区内工作辊的轮廓曲线，可得出以下四个主要结论：(1)首先基于工业多机架生产情况对传统方法的准确性进行了评估利益冲突作者声明不存在利益冲突。作者贡献韩国民：数据整理；方法学；软件；验证；撰写初稿。李洪波：概念构思；写作审阅与编辑；监督；项目管理。刘玉金：数据整理；写作审阅与编辑；监督。张杰：概念构思；监督。颜东佳：写作审阅与编辑；监督。赵建伟：监督；项目管理。杨泉：写作审阅与编辑；监督。游学昌：数据整理。胡志远：数据整理。刘雷：CRediT作者贡献说明韩国民：写作——初稿，验证，软件，方法学，数据整理。李洪波：写作——审阅与编辑，监督，项目管理，概念构思。刘玉金：写作——审阅与编辑，监督，数据整理。张杰：监督，概念构思。颜东佳：写作——审阅与编辑，监督。赵建伟：监督，项目管理。杨泉：写作——审阅与编辑，监督。游学昌：数据整理。胡志远：数据韩国民|李洪波|刘玉金|张杰|颜东佳|赵建伟|杨泉|游学昌|胡志远|刘雷