一种六辊冷轧机的轧制稳定性预测方法与流程

1.本发明涉及轧制过程自动化生产技术领域，具体涉及一种六辊冷轧机的轧制稳定性预测方法。


背景技术：

2.轧机振动是板带材生产中普遍存在和亟待解决的问题。高速轧制薄规格高强钢时，由于工艺参数、设备状态和控制系统的强耦合、非线性，轧机常常出现各种通过调整工艺参数难以消除的异常振动，例如轧机垂直方向的自激振动。轧机垂直方向的自激振动会导致带钢和轧辊表面出现周期性振纹，严重影响了产品质量；而且轧机垂直方向的自激振动也会加剧轧辊和轴承的磨损，甚至造成断辊断带，威胁工作人员的生命安全。轧机垂直方向的自激振动产生的原因是在轧机结构动态变化和轧制过程相互作用中，工艺参数改变造成轧制界面等效阻尼和刚度发生了变化，导致界面等效阻尼变小。若轧件-轧机系统总阻尼为负，则系统不断从传动装置吸收能量，使轧辊振幅不断增大，将导致轧制过程失稳。
3.针对高速轧制过程中不断出现的轧机垂直方向的自激振动问题，研究人员做了许多的相关研究。但这些研究主要存在有两个方面不足：(1)假定振动过程中轧制界面的摩擦状态不变。然而，实际上轧机垂直方向的自激振动会造成界面润滑状态发生周期性变化，这种变化会随着振动加剧而更加明显，若不考虑此变化将会导致对自激振动判定的准确度降低。(2)利用卡尔曼微分方程计算振动发生时的轧制力，由于未考虑轧辊垂向振动速度的影响，所得结果不够精确。


技术实现要素：

4.针对上述现有技术存在的不足，本发明提供了一种六辊冷轧机的轧制稳定性预测方法，可根据现有轧制规程和轧机结构参数对轧制界面润滑状态变化和动态轧制力进行计算，并能够更精准地预测出轧制过程中的轧机稳定性，进而避免因轧制速度过高或轧制规程制定不合理引发的轧机垂直方向的自激振动，达到提高带钢表面质量和轧制过程稳定运行的目的。
5.为实现上述目的，本发明提供的技术方案为：
6.一种六辊冷轧机的轧制稳定性预测方法，该方法包括如下步骤：
7.步骤1：采集相关参数，包括带钢参数、润滑油参数、轧制工艺参数以及轧机结构参数；
8.步骤2：根据带钢参数、轧辊辊径和轧辊垂向振动速度计算变形区动态接触弧长，沿轧制方向将变形区进行离散化处理，并计算所述离散化处理后获得的各微元体的平均变形抗力；
9.步骤3：利用一维reynolds方程计算动态入口油膜厚度；
10.步骤4：结合步骤3中获得的动态入口油膜厚度和christensen粗糙度分布假设，计算变形区摩擦应力分布；
11.步骤5：对卡尔曼微分方程进行改进，将步骤4中获得的摩擦应力分布代入改进的卡尔曼微分方程求解由轧辊垂向振动引发的轧制力波动量δp1；
12.步骤6：根据机架间张力关系计算振动导致的后张力变化量以及由后张力变化引发的轧制力波动量δp2；
13.步骤7：根据各轧机各部件间的受力关系和轧制力波动总量δp＝δp1 δp2，建立轧机系统的垂向振动动力学方程并求解，得到轧辊位移和速度响应，从而预测出轧制过程的稳定性。
14.进一步地，根据所述的六辊冷轧机的轧制稳定性预测方法，所述带钢参数包括：带钢牌号、带钢宽度和热轧来料厚度；所述润滑油参数包括润滑油黏度和黏压系数；所述轧制工艺参数包括：机架间前后张力、各道次轧制速度、各道次带钢入口速度、各道次带钢出入口厚度；所述轧机结构参数包括：轧辊质量、轧辊材质、轧辊辊径、轧辊长度、轧机刚度系数、轧机各部件阻尼系数、牌坊质量；所述的轧机刚度系数包括轧辊刚度系数和牌坊刚度系数。
15.根据所述的六辊冷轧机的轧制稳定性预测方法，所述步骤2进一步包括如下步骤：
16.步骤2.1：根据带钢出入口厚度、轧辊辊径和轧辊垂向振动速度求解变形区动态接触弧长；
17.步骤2.2：沿轧制方向将变形区进行离散化处理，获得若干微元体；
18.步骤2.3：根据带钢材质和微元体厚度利用变形抗力模型计算得到各微元体的平均变形抗力。
19.进一步地，根据所述的六辊冷轧机的轧制稳定性预测方法，所述动态接触弧长l的计算公式如下：
[0020][0021]
上式中，l为动态接触弧长；r为轧辊压扁半径；y
in
和y
out
为带钢入口和出口厚度；θ为咬入角变化量；vy为轧辊垂向振动速度，速度向上为正；v
out
为带钢出口速度。
[0022]
进一步地，根据所述的六辊冷轧机的轧制稳定性预测方法，所述步骤3中所述利用一维reynolds方程计算动态入口油膜厚度的方法为：
[0023]
考虑挤压效应的一维reynolds方程如下式所示：
[0024][0025]
上式中，h1为入口区油膜厚度；xf为距入口出距离；p为变形区轧制压力分布；η为不同压力下的润滑油黏度；v为带钢和轧辊的平均速度；t为时间；
[0026]
由于入口区油膜压力较小，润滑油黏压公式采用barus公式，如下所示：
[0027]
η＝η0e
γp
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0028]
引入无量纲参数φ：
[0029]
[0030]
上式中，γ为barus公式的黏压系数；p为变形区轧制压力分布；η0为大气压下的润滑油黏度；
[0031]
对一维reynolds方程积分并用vy＝0时的稳态结果替换积分常数，可以得到：
[0032][0033]
上式中，α为咬入角；为咬入角变化速率；σb为带钢后张力；当xf趋于无穷大时，润滑油油膜压力趋于0，因此可得到如下边界条件1：
[0034]
xf＝∞,h1＝∞，φ＝1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0035]
将边界条件1代入的表达式并整理可得：
[0036][0037][0038]
根据在入口区和变形区交界处可根据tresca屈服准则p＝σ
s-σb计算润滑油油膜压力，可得到如下边界条件2：
[0039][0040]
将边界条件2代入φ的表达式并整理可得：
[0041][0042][0043]
其中，为入口油膜厚度变化速率；h
0,d
为动态入口油膜厚度；δt为时间步长；h0为稳态时的入口油膜厚度，可由下式确定：
[0044][0045]
其中，v
in
为带钢入口速度；vr为轧制速度；l0为稳态时的变形区长度。
[0046]
进一步地，根据所述的六辊冷轧机的轧制稳定性预测方法，所述步骤4中所述的计算变形区摩擦应力分布的方法为：
[0047]
根据步骤3中获得的动态入口油膜厚度以及体积不变原理，变形区油膜厚度分布h(x)可以下式表示：
[0048][0049]
上式中，h(x)为变形区油膜厚度分布；vr为轧制速度；v
in
为带钢入口速度；vs为带钢沿轧制方向速度分布；h
0,d
为动态入口油膜厚度；
[0050]
由christensen粗糙度分布假设可知，实际接触面积比ac和平均油膜厚度h
t
可分别表示为：
[0051][0052][0053]
上式中，δ为粗糙度分布；z＝h/3rq为无量纲参数；f(δ)为概率密度函数，可表示为：
[0054][0055]
上式中，rq为带钢和轧辊的综合表面粗糙度；
[0056]
最后，变形区摩擦应力分布τ可表示为：
[0057][0058]
上式中，τ为变形区总摩擦应力分布；τa为粗糙接触产生的摩擦应力；τf为流体润滑产生的摩擦应力；k为材料的剪切强度。
[0059]
进一步地，根据所述的六辊冷轧机的轧制稳定性预测方法，所述步骤5进一步包括如下步骤：
[0060]
步骤5.1：对变形区微元体的受力进行分析后对变形区微元体沿轧制方向列静力平衡关系方程，由静力平衡关系方程获得改进的卡尔曼微分方程；
[0061]
步骤5.2：将步骤4中获得的摩擦应力分布代入改进的卡尔曼微分方程，积分后得到由轧辊垂向振动引发的轧制力波动量δp1。
[0062]
8、根据权利要求1所述的六辊冷轧机的轧制稳定性预测方法，其特征在于，所述步骤7包括如下步骤：
[0063]
步骤7.1：根据六辊冷轧机的二分之一简化模型以及轧辊、轧件和牌坊间的受力关系，再结合机械振动理论，建立轧机系统的垂向振动动力学方程；
[0064]
步骤7.2：采用newmark-beta法对轧机系统的垂向振动动力学方程进行求解，并以轧辊垂向速度作为下一时刻计算过程的输入量，以轧辊垂向位移作为判断轧机稳定性的依据，若轧辊位移曲线收敛，则轧机稳定，若轧辊位移曲线发散，则轧机不稳定。
[0065]
总体而言，通过本发明所构思的以上技术方案较现有技术具有以下有益效果：首先，本发明方法考虑了轧制界面润滑状态变化和轧辊振动造成的轧制力变化，轧制力计算结果更精确；其次，本发明方法在现有轧制规程或实时采集数据的基础上，结合压靠实验，即可预测大部分六辊冷轧机的轧制过程稳定性，具有广泛的适用性；再次，采用本发明方法可以避免因轧制速度过高或轧制规程制定不合理引发的轧机垂直方向的自激振动，达到提高带钢表面质量和提升企业效益；最后，本发明方法是基于轧制理论和机械动力学的仿真模拟，可有效避免实验导致的设备损耗及破坏，降低企业成本。
附图说明
[0066]
图1为本实施方式六辊冷轧机的轧制稳定性预测方法的流程示意图；
[0067]
图2为轧辊与带钢间的油膜厚度分布示意图；
[0068]
图3为本实施方式微元体受力分析示意图，其中图(a)为轧辊垂向振动速度vy00时的微元体受力分析图；图(b)为轧辊垂向振动速度vy《0时的微元体受力分析图；
[0069]
图4中(a)图为ucm六辊冷轧机结构示意图；(b)图为(a)图所示六辊冷轧机的二分之一简化模型示意图；
[0070]
图5为不同工艺参数下工作辊位移响应的预测曲线图，其中图(a)为不同压下率下工作辊位移响应的预测曲线图；图(b)为不同后张力下工作辊位移响应的预测曲线图；图(c)为不同粗糙度下工作辊位移响应的预测曲线图；图(d)为不同黏度下工作辊位移响应的预测曲线图；图(e)为不同轧制速度下工作辊位移响应的预测曲线图。
具体实施方式
[0071]
为了使本发明的目的、技术方案及优势更加清晰，下面结合附图和具体实施例对本发明做进一步详细说明。此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0072]
本发明方法的核心思路包括：1、考虑入口油膜挤压效应，将轧辊垂向振动速度引入油膜厚度计算公式，获得动态入口油膜厚度，并结合粗糙度分布假设，计算变形区摩擦应力分布随时间的变化情况；2、考虑轧辊垂向振动速度的卡尔曼微分方程推导，并带入变形区摩擦应力分布，计算动态轧制力及由轧辊垂向振动引发的轧制力波动量；3、根据轧辊、轧件和牌坊间的受力关系，建立轧机系统的垂向振动动力学方程，然后采用newmark-beta法求解，并以轧辊垂向位移作为判断轧机稳定性的依据，若轧辊位移曲线收敛，则轧机稳定，若轧辊位移曲线发散，则轧机不稳定。
[0073]
本实施方式中以某厂的1450mm ucm六辊冷连轧机组为例，对其轧制稳定性进行预测，ucm六辊冷轧机结构如图4中(a)图所示，其中轧机轧辊均为平辊。图1是本实施方式六辊冷轧机的轧制稳定性预测方法的流程示意图，如图1所示，所述六辊冷轧机的轧制稳定性预测方法包括如下步骤：
[0074]
步骤1：采集相关参数，包括带钢参数、润滑油参数、轧制工艺参数以及轧机结构参数；
[0075]
所述带钢参数包括：带钢牌号、带钢宽度和热轧来料厚度；所述润滑油参数包括润滑油黏度和黏压系数；所述轧制工艺参数包括：机架间前后张力、各道次轧制速度、各道次带钢入口速度、各道次带钢出入口厚度；所述轧机结构参数包括：轧辊质量、轧辊材质、轧辊辊径、轧辊长度、轧机刚度系数、轧机各部件阻尼系数、牌坊质量；
[0076]
在本实施例中，具体是从冷轧产线上的一级控制系统和二级控制系统获取所有带钢参数、所有轧制工艺参数和轧机结构参数中的轧辊质量、轧辊材质、轧辊辊径、轧辊长度、牌坊质量。
[0077]
所述的轧机刚度系数包括轧辊刚度系数和牌坊刚度系数，通过现场的压靠实验获得所述轧机刚度系数，然后利用hertz接触理论计算获得轧辊刚度系数，最后根据hooke定律计算获得牌坊刚度系数。
[0078]
在本实施例中，首先通过在现场对某1450mm ucm六辊冷连轧机进行压靠实验，获得轧机刚度系数为k＝4.4
×
109n/m。然后根据hertz接触理论，两轧辊压缩时的轧辊刚度系数ki可用式(1)表示：
[0079][0080]
上式中，ki为轧辊刚度系数，单位n/m；p为轧制力，单位kn；e为轧辊弹性模量，
[0081]
经验取值为2.1
×
1011pa；v为轧辊泊松比，经验取值为0.3；d1和d2为两轧辊辊径，单位mm。
[0082]
轧辊刚度系数确定后，在轧机二分之一简化模型中上半部分牌坊刚度系数ks可根据式(2)所示的hooke定律计算：
[0083][0084]
上式中，k为轧机刚度系数；kw为工作辊刚度系数，k
im
为中间辊刚度系数，kb为支撑辊刚度系数，ks为上半部分牌坊刚度系数，单位n/m。
[0085]
所述的轧机各部件阻尼系数根据rayleigh阻尼公式及轧辊刚度系数和牌坊刚度系数获得。rayleigh阻尼是常见的结构阻尼构造方法，其假设结构的阻尼矩阵c是质量矩阵m和刚度矩阵k的线性组合，即：
[0086]
c＝β1m β2k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0087][0088][0089]
上式中，c为阻尼矩阵，单位n
·
s/m；m为质量矩阵，单位kg；k为刚度矩阵，单位n/m；β1和β2为比例系数；ω1和ω2为固有频率，经验取值分别为100hz和500hz；ξ1和ξ2为阻尼比，经验取值为0.03。
[0090]
在本实施例中，带钢参数和润滑油参数如表1所示，轧制工艺参数如表2所示，轧机结构参数如表3所示。
[0091]
表1带钢参数和润滑油参数
[0092][0093]
表2轧制工艺参数
[0094][0095][0096]
表3轧机结构参数
[0097][0098]
步骤2：根据带钢参数、轧辊辊径和轧辊垂向振动速度计算变形区动态接触弧长，沿轧制方向将变形区进行离散化处理，并计算所述离散化处理后获得的各微元体的平均变形抗力；
[0099]
步骤2.1：根据带钢出入口厚度、轧辊辊径和轧辊垂向振动速度求解变形区动态接触弧长；
[0100]
由于轧辊垂向振动会造成变形区变化，因此在传统轧制理论的基础上进行修正，得到动态接触弧长l的计算公式如下：
[0101][0102]
上式中，l为动态接触弧长，单位mm；r为轧辊压扁半径，单位mm；y
in
和y
out
为带钢入
口和出口厚度，单位mm；θ为咬入角变化量，单位rad；vy为轧辊垂向振动速度，速度向上为正，单位m/s；v
out
为带钢出口速度，单位m/s。
[0103]
轧辊压扁半径根据式(7)所示的hitchcock公式计算获得：
[0104][0105]
上式中，r0为轧辊初始半径，单位mm；ew为工作辊弹性模量，取值为2.1
×
10
11
pa；pi为单位长度上的轧制力n/m。
[0106]
步骤2.2：为提高计算精度，沿轧制方向将变形区进行离散为1000份，由于每一份很小，因此称之为微元体；
[0107]
步骤2.3：根据带钢材质和微元体厚度利用变形抗力模型计算得到各微元体的平均变形抗力；
[0108]
在本实施例中，带钢的牌号为q195，平均变形抗力σs采用如下式(8)至式(9)计算：
[0109][0110][0111][0112]
其中，σs为平均变形抗力，单位mpa；各系数的经验取值分别为a＝498mpa、b＝136mpa、c＝0.2、d＝5；ε
σ
为累计变形量；y0为热轧来料厚度；为道次平均厚度。
[0113]
步骤3：利用一维reynolds方程计算动态入口油膜厚度；
[0114]
轧辊与带钢间的油膜厚度分布如图2所示，其中，y
in
和y
out
分别为带钢入口和出口厚度，单位mm；r0和r为轧辊初始半径和压扁半径，单位mm；h0和h1分别为入口处和入口区油膜厚度，单位mm；α为咬入角，单位rad；xf为距入口出距离，单位mm。利用一维reynolds方程计算动态入口油膜厚度具体按照如下方法进行；
[0115]
考虑挤压效应的一维reynolds方程如下式所示：
[0116][0117]
由于入口区油膜压力较小，润滑油黏压公式采用barus公式，如下所示：
[0118]
η＝η0e
γp
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0119]
引入无量纲参数φ：
[0120][0121]
上式中，t为时间，单位s；v为带钢和轧辊的平均速度，单位m/s；γ为barus公式的黏压系数，单位mpa-1
；p为变形区轧制压力分布，单位mpa；η0为大气压下的润滑油黏度，单位pa
·
s；η为不同压力下的润滑油黏度，单位pa
·
s。
[0122]
对一维reynolds方程积分并用vy＝0时的稳态结果替换积分常数，可以得到：
[0123][0124]
其中，为咬入角变化速率，单位rad/s；σb为带钢后张力，单位mpa。当xf趋于无穷大时，润滑油油膜压力趋于0，因此可得到如下边界条件1：
[0125]
xf＝∞,h1＝∞，φ＝1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0126]
将边界条件1代入的表达式并整理可得：
[0127][0128][0129]
又因在入口区和变形区交界处，润滑油油膜压力可根据tresca屈服准则p＝σ
s-σb计算，因此可得到如下边界条件2：
[0130][0131]
将边界条件2代入φ的表达式并整理可得：
[0132][0133][0134]
其中，为入口油膜厚度变化速率，单位mm/s；h
0,d
为动态入口油膜厚度，单位mm；δt为时间步长，单位s；h0为稳态时的入口油膜厚度，单位mm，可由下式确定：
[0135][0136]
其中，v
in
为带钢入口速度，单位m/s；vr为轧制速度，单位m/s；l0为稳态时的变形区长度，单位mm。
[0137]
步骤4：结合步骤3中获得的动态入口油膜厚度和粗糙度分布假设，计算变形区摩擦应力分布；
[0138]
根据步骤3中获得的动态入口油膜厚度以及体积不变原理，变形区油膜厚度分布h(x)可以式(22)表示：
[0139][0140]
上式中，h(x)为变形区油膜厚度分布，单位mm；vs为带钢沿轧制方向速度分布，单位m/s。
[0141]
由christensen粗糙度分布假设可知，实际接触面积比ac和平均油膜厚度h
t
可分别表示为：
[0142]
[0143][0144]
上式中，ac为实际接触面积比；h
t
为平均油膜厚度，单位mm；z＝h/3rq为无量纲参数；f(δ)为概率密度函数，可表示为：
[0145][0146]
上式中，δ为粗糙度分布，单位μm；rq为带钢和轧辊的综合表面粗糙度，单位μm。
[0147]
最后，变形区摩擦应力分布τ可表示为：
[0148][0149]
上式中，τ为变形区总摩擦应力分布，单位mpa；τa为粗糙接触产生的摩擦应力，单位mpa；τf为流体润滑产生的摩擦应力，单位mpa；k为材料的剪切强度，单位mpa。
[0150]
步骤5：将步骤4中获得的摩擦应力分布代入修正的卡尔曼微分方程求解由轧辊垂向振动引发的轧制力波动量；
[0151]
步骤5.1：对变形区微元体沿轧制方向列静力平衡关系方程；
[0152]
变形区微元体的受力分析如图3所示，其中图(a)为轧辊垂向振动速度vy00时的微元体受力分析图；图(b)为轧辊垂向振动速度vy《0时的微元体受力分析图。图中，dx为微元体宽度，单位mm；y，dy和δy为微元体厚度、厚度增量和轧辊垂向振动引发的厚度变化量，单位mm；σ
x
和dσ
x
为微元体所受应力和应力增量。以(a)图为例，将微元体每条边上的力做轧制方向的投影可得到：
[0153]fabx
＝(σ
x
dσ
x
)(y dy)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)
[0154]fefx
＝-σ
x
(y 2δy)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(28)
[0155][0156]
整理后可得：
[0157][0158]
根据几何关系：
[0159][0160]
经过整理后，改进的卡尔曼微分方程如下式：
[0161][0162]
其中，k
p
＝1.155σs为材料的平面变形抗力，单位mpa；“ ”为后滑区，
“‑”
为前滑区。
[0163]
步骤5.2：将步骤4中获得的摩擦应力分布代入改进的卡尔曼微分方程，积分后得到由轧辊垂向振动引发的轧制力波动量δp1；
[0164]
将步骤4中获得的摩擦应力分布代入步骤5.1中获得的改进卡尔曼微分方程，对其沿变形区积分后，可得到由轧辊垂向振动引发的动态轧制力pd，则由轧辊垂向振动引发的轧制力波动量δp1可表示为：
[0165]
δp1＝p
d-psꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(32)
[0166]
其中，ps为vy＝0时的稳态轧制力，单位是kn。
[0167]
步骤6：根据机架间张力关系计算振动导致的后张力变化量以及由后张力变化引发的轧制力波动量δp2，具体按照如下方法进行：
[0168]
步骤6.1：根据机架间张力关系计算振动导致的后张力变化量；
[0169]
轧辊振动会导致带钢入口速度发生改变，进而通过机架间张力关系导致后张力发生变化，后张力变化量δσb的表达式如下所示：
[0170][0171]
其中，es为带钢弹性模量，取值为2.1
×
10
11
pa；l为机架间距离，取值为5m；v
in,i
为当前机架入口速度，单位m/s；v
out,i-1
为前一机架出口速度，单位m/s。
[0172]
步骤6.2：计算由后张力变化引发的轧制力波动量δp2；
[0173]
后张力变化量引发的轧制力波动量δp2为：
[0174][0175]
其中，q
p
为应力状态系数，w为轧件宽度，取值为1000mm。
[0176]
步骤7：根据各轧机各部件间的受力关系和轧制力波动总量δp＝δp1 δp2，建立轧机系统的垂向振动动力学方程并求解，得到轧辊位移和速度响应，从而预测出轧制过程的稳定性。
[0177]
步骤7.1：根据图4中(b)图示出的六辊冷轧机的二分之一简化模型，轧辊、轧件和牌坊间的受力关系，再结合机械振动理论，建立轧机系统的垂向振动动力学方程如下：
[0178][0179]
上式中，ms、mb、m
im
和mw分别为轧机上半部分牌坊质量、支撑辊质量、中间辊质量和工作辊质量，单位kg；cb、c
im
和cw分别为支撑辊阻尼、中间辊阻尼和工作辊阻尼，单位n
·
s/m；ks、kb、k
im
和kw分别为轧机上半部分牌坊刚度、支撑辊刚度、中间辊刚度和工作辊刚度，单位n/m；xs、xb、x
im
和xw分别为轧机上半部分牌坊位移、支撑辊位移、中间辊位移和工作辊位移，单位m；和分别为轧机上半部分牌坊速度、支撑辊速度、中间辊速度和工作辊速度，单位m/s；和分别为轧机上半部分牌坊加速度、支撑辊加速度、中间辊加速度和工作辊加速度，单位m/s2；δp为轧制力波动总量，单位kn。
[0180]
步骤7.2：采用newmark-beta法对轧机系统的垂向振动动力学方程进行求解；轧辊
垂向速度作为下一时刻计算过程的输入量；轧辊垂向位移作为判断轧机稳定性的依据，若轧辊位移曲线收敛，则轧机稳定，若轧辊位移曲线发散，则轧机不稳定。
[0181]
不同工艺参数下的轧机稳定性预测效果如图5所示。在图5(a)中，压下率增大导致轧辊位移曲线由收敛变为发散，轧机不稳定且会发生自激振动；在图5(b)中，后张力增大则使轧辊位移曲线变得收敛，提高了轧机的稳定性，但是效果并不显著；在图5(c)和(d)中，粗糙度减小和黏度增大均会提高轧制界面的润滑性能，进而降低界面的摩擦耗能能力，导致轧机不稳定；在图5(e)中，随着轧制速度的不断升高，轧辊位移曲线从收敛变为恒振幅再到发散，表明轧制速度的升高会使轧机失稳。由此可见，本发明方法可以通过生产前的模拟计算来判断所制定的轧制规程是否合理，从而避免生产事故和设备损伤等问题。
[0182]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解；其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；因而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。