CN114091308B - 基于三维模型的六辊冷轧机临界振动速度预测方法 - Google Patents

Abstract

本发明提供基于三维模型的六辊冷轧机临界振动速度预测方法，属于轧制过程自动化技术领域，基于三维的六辊冷轧机模型预测临界振动速度，考虑到轧辊应视为短粗梁且需考虑剪切变形的影响，选择Timoshenko梁，同时对节点位移矢量采用Hermite插值；通过对轧件、轧辊和牌坊间的受力分析，可以建立轧机‑轧件系统的垂向振动动力学方程，采用Newmark‑Beta法进行求解，可得到特定速度下的轧辊位移响应曲线，若位移响应曲线的幅值恒定，则该速度为轧机的临界振动速度；本发明不仅可以研究二维轧制工艺参数对轧制过程稳定性的影响，还可以分析弯蹿辊等宽向参数对临界轧制速度的影响，在轧制规程制定阶段就预测出轧机的临界振动速度，可以为工艺参数优化提供理论支撑。

Description

基于三维模型的六辊冷轧机临界振动速度预测方法

技术领域

本发明属于轧制过程自动化技术领域，具体涉及一种基于三维模型的六辊冷轧机临界振动速度预测方法。

背景技术

薄规格高速轧制时，由于工艺参数、设备状态和控制系统的多变量、强耦合和非线性等特点，轧机常常出现各种通过调整工艺参数难以消除的异常振动，是制约轧制过程高速度和高精度的瓶颈问题。保证轧制过程高效稳定运行的关键是轧制工艺参数的合理制定。在高速轧制过程中，当轧制速度超过某一临界值时，轧机就会发生剧烈的自激振动。由于轧制工艺参数种类繁多且相互耦合，还具有很强的非线性特点，不同轧制工艺参数下的临界振动速度往往相差很大。轧机振动产生的根本原因是工艺参数改变造成了轧制界面等效阻尼减小，导致轧件-轧机系统不断从传动装置吸收能量、轧辊振幅不断增大，最终出现轧制过程失稳情况。

国内外研究人员针对此类轧机振动问题做了许多的相关研究。中国发明专利“一种六辊冷轧机的极限轧制速度预测方法”(专利号：CN 109078989 A)提出了一种预测轧机极限轧制速度的方法，该方法通过建立实际等效阻尼系数与轧制速度的关系，然后计算带钢出口张力波动导致的负阻尼效应大小，最后利用阻尼系数与轧制速度的关系得到轧机极限轧制速度。中国发明专利“高速薄板轧机自激振动预测与抑制的工艺方法”(专利号：CN106734194 A)发明了一种根据轧制速度与轧机自激振动诱发条件之间关系的自激振动预测与抑制方法，该方法首先考虑轧件塑性变形和上下辊系非对称性建立轧机结构动力学模型，然后对轧制力解析解经积分运算和泰勒(TAYLOR)公式展开得到轧制力增量，再利用劳斯(ROUTH)稳定性判据求解临界轧制速度，最后根据不同工艺参数下轧制速度的稳定性裕量及前滑值异常变化来预测自激振动，并给出优化措施。中国发明专利“一种冷连轧机组振动抑制方法”(专利号：CN 105522000 A)提出了一种通过控制乳化液喷射量来抑制轧机振动的方法，该方法根据在第4或第5机架上振动传感器采集的信号能量大小来判断轧机是否发生振动，然后由振动信号能量或前滑值决定第4或第5机架入口乳化液喷射量。中国期刊文章“轧机整体耦合建模问题研究”(机械工程学报,2015,51(14):46-53)分析了不同类型振动产生机理及其相互关系，并进一步建立了轧机垂直-水平-扭转以及轧件水平颤振的耦合振动结构模型，在此基础上构建了轧机颤振整体耦合动力学模型。外国期刊文章“Dynamic analysis of cold-rolling process using the finite-element method”(Journal of Manufacturing Science and Engineering,2016,138(4):041002)提出了一种基于有限元方法的四辊轧机简化模型，该方法将轧辊等效为Euler-Bernoulli梁，研究了轧机振动对带钢出口形貌、应力和轧制压力沿宽向分布情况的影响，还分析了带钢宽度对轧机振动的影响。

上述中国专利和期刊中研究的主要不足在于，他们将轧机各部件简化为质点，忽略了弯窜辊和带钢宽度等宽向参数对轧机振动的影响。而外国期刊中研究的主要不足在于，将轧辊简化为Euler-Bernoulli梁，忽略了轧辊垂向振动时的剪切变形，而且轧制力计算未考虑轧辊振动速度的影响，模型精度有所欠缺。

发明内容

针对现有技术的不足，本发明提出一种基于三维模型的六辊冷轧机临界振动速度预测方法，包括：

步骤1：获取冷连轧产线实际生产中一段时间内的生产过程参数；

步骤2：根据六辊冷轧机中各个轧辊的有限元模型和辊间接触刚度，建立六辊冷轧机的整体刚度、质量和阻尼矩阵；

步骤3：将变形区沿轧制方向和宽度方向离散为若干份，计算动态轧制力；

步骤4：根据轧辊和轧件间受力关系以及计算得到的动态轧制力，建立轧机-轧件系统的垂向振动动力学方程，求解得到轧辊位移响应用于预测轧机临界振动速度。

所述生产过程参数包括轧机结构参数、轧制工艺参数、带钢参数以及润滑油参数；

所述轧机结构参数包括：轧辊材质、弹性模量、泊松比、密度、辊身长度和直径、辊颈长度和直径以及牌坊材质、质量和尺寸；

所述轧制工艺参数包括：机架间前后张力、各道次轧制速度、各道次带钢入口速度、各道次带钢出入口厚度、工作辊弯辊力以及中间辊弯辊力和横移量；

所述带钢参数以及润滑油参数包括：带钢牌号和宽度、来料厚度以及润滑油黏度和黏压系数。

所述步骤2包括：

步骤2.1：建立六辊冷轧机中的各个轧辊的有限元模型；

步骤2.2：计算轧辊辊间接触刚度；

步骤2.3：根据轧辊的有限元模型和辊间接触刚度，建立六辊冷轧机的整体刚度、质量和阻尼矩阵。

所述步骤3包括：

步骤3.1：根据带钢出入口厚度、轧辊辊径和轧辊垂向振动速度计算得到变形区动态接触弧长l_d：

其中，R为轧辊压扁半径；y_in和y_out为带钢入口和出口厚度；θ为咬入角变化量；v_y为轧辊垂向振动速度，速度向上为正；

步骤3.2：根据带钢材质和微元体厚度利用变形抗力模型计算得到各微元体的平均变形抗力；

步骤3.3：计算变形区摩擦应力分布；

步骤3.4：建立变形区微元体的力平衡微分方程；

步骤3.5：将步骤3.3所得的摩擦应力分布代入力平衡微分方程，对其沿轧制方向和宽度方向积分得到动态轧制力。

所述步骤4包括：

步骤4.1：建立轧机-轧件系统的垂向振动动力学方程：

其中，x、和/>分别为轧辊梁单元节点的位移、速度和加速度矢量；M为轧机整体质量矩阵；C为轧机整体阻尼矩阵；K_z为轧机整体刚度矩阵；K为支撑辊、中间辊和工作辊梁单元模型组合而成的总刚度矩阵；K_iw为中间辊和工作辊之间的辊间接触刚度矩阵；K_bi为支撑辊和中间辊之间的辊间接触刚度矩阵；F_iw为工作辊和中间辊间的接触压力分布，F_iw＝K_iwx；F_bi为中间辊和支撑辊间的接触压力分布，F_bi＝K_bix；F为轧制力宽向分布；β₁和β₂为比例系数；ξ₁和ξ₂为阻尼比；ω₁和ω₂为频率；

步骤4.2：对动力学方程采用Newmark-Beta法进行求解，可得到特定速度下的轧辊位移响应曲线，位移响应曲线的幅值恒定时对应的速度即为轧机的临界振动速度。

所述步骤2.1包括：

步骤2.1.1：将六辊冷轧机中的各个轧辊简化为梁并离散为若干个单元；

步骤2.1.2：根据插值函数确定各个单元的节点位移矢量的形函数表达式；

步骤2.1.3：利用虚功原理推导单元刚度矩阵并进行组集；

步骤2.1.4：利用虚功原理推导单元质量矩阵并进行组集。

所述步骤2.3包括：

步骤2.3.1：利用压靠实验得到轧机刚度系数，采用有限元分析软件计算牌坊刚度系数；

步骤2.3.2：结合计算得到的辊间接触刚度确定支撑辊刚度系数，最后得到轧机整体刚度矩阵：

K_z＝K+K_iw+K_bi

其中，K_z为轧机整体刚度矩阵；K为支撑辊、中间辊和工作辊梁单元模型组合而成的总刚度矩阵；K_iw为中间辊和工作辊之间的辊间接触刚度矩阵；K_bi为支撑辊和中间辊之间的辊间接触刚度矩阵：

其中，n＝n_w+n_i+n_b，n₁～n₆为行和列的序号，根据以下条件选择：

其中，n_w、n_i和n_b分别为工作辊、中间辊和支撑辊的节点数目；n_s为中辊间横移量对应的节点数目；

步骤2.3.3：利用Rayleigh阻尼公式、轧机整体刚度矩阵和质量矩阵获得轧机整体阻尼矩阵。

本发明的有益效果是：

本发明提出了基于三维模型的六辊冷轧机临界振动速度预测方法，基于三维的六辊冷轧机模型预测临界振动速度，考虑到轧辊应视为短粗梁且需考虑剪切变形的影响，选择Timoshenko梁，同时对节点位移矢量采用Hermite插值；通过对轧件、轧辊和牌坊间的受力分析，可以建立轧机-轧件系统的垂向振动动力学方程，采用Newmark-Beta法进行求解，可得到特定速度下的轧辊位移响应曲线，若位移响应曲线的幅值恒定，则该速度为轧机的临界振动速度；本发明方法所提供的基于三维模型的六辊冷轧机临界振动速度预测方法，首先，通过改进的Timoshenko梁构建三维六辊冷轧机模型并耦合动态轧制力计算模型，提高了计算速度和精度。其次，本发明不仅可以研究二维轧制工艺参数对轧制过程稳定性的影响，还可以分析弯蹿辊等宽向参数对临界轧制速度的影响，与其他方法相比分析参数更加全面合理。再者，采用本发明提供的方法可以在轧制规程制定阶段就预测出轧机的临界振动速度，为工艺参数优化提供理论支撑。

附图说明

图1为本发明中基于三维模型的六辊冷轧机临界振动速度预测方法流程图；

图2为本发明中轧机上辊系简化模型示意图；

图3为本发明中Timoshenko梁单元示意图；

图4为本发明中辊间接触刚度计算示意图；

图5为本发明中三维变形区网格划分示意图；

图6为本发明中工作辊弯辊力对临界轧制速度的影响及对应带钢出口形貌图，其中(a)为工作辊弯辊力对轧机临界轧制速度的影响图，(b)为对应临界振动速度下的带钢出口形貌图；

图7为本发明中中间辊横移量对临界轧制速度的影响及对应带钢出口形貌图，其中(a)为中间辊横移量对轧机临界轧制速度的影响图，(b)为对应临界振动速度下的带钢出口形貌图；

图8为本发明中中间辊弯辊力对临界轧制速度的影响及对应带钢出口形貌图，其中(a)为中间辊弯辊力对轧机临界轧制速度的影响图，(b)为对应临界振动速度下的带钢出口形貌图。

具体实施方式

下面结合附图和具体实施实例对发明做进一步说明。本发明的目的是提供一种基于三维模型的六辊冷轧机临界振动速度预测方法，本方法采用一种改进的Timoshenko梁，只需较少单元数就可分析弯蹿辊等宽向参数对临界轧制速度的影响，为制定更准确合理的轧制规程提供依据，避免因轧制速度过高或工艺参数设定不合理引发的轧机自激振动，以达到冷连轧机组高效和稳定运行的目的。

本实施例中以某厂的1450mm UCM六辊冷连轧机组为例，对其不同弯蹿辊参数下的临界振动速度进行预测，支撑辊、中间辊和工作辊均为平辊。

如图1所示，一种基于三维模型的六辊冷轧机临界振动速度预测方法，包括：

步骤1：获取冷连轧产线实际生产中一段时间内的生产过程参数；所述生产过程参数包括轧机结构参数、轧制工艺参数、带钢参数以及润滑油参数；

所述轧机结构参数包括：轧辊材质、弹性模量、泊松比、密度、辊身长度和直径、辊颈长度和直径以及牌坊材质、质量和尺寸，所需参数可根据轧辊设计图纸获得；

所述轧制工艺参数包括：机架间前后张力、各道次轧制速度、各道次带钢入口速度、各道次带钢出入口厚度、工作辊弯辊力以及中间辊弯辊力和横移量，所需参数可从冷连轧生产线的一级和二级控制系统获得；

所述带钢参数以及润滑油参数包括：带钢牌号和宽度、来料厚度以及润滑油黏度和黏压系数，其中带钢参数可从冷轧生产线的一级和二级控制系统获得，润滑油参数可由实验测定。

本实施例中，轧机结构参数如表1所示：

表1轧机结构参数

本实施例中，轧制工艺参数如表2所示：

表2轧制工艺参数

本实施例中，带钢参数以及润滑油参数如表3所示：

表3带钢参数以及润滑油参数

步骤2：根据六辊冷轧机中各个轧辊的有限元模型和辊间接触刚度，建立六辊冷轧机的整体刚度、质量和阻尼矩阵；包括：

步骤2.1：建立六辊冷轧机中的各个轧辊的有限元模型；如图2所示，图中，WRB表示工作辊弯辊力，IRB表示中间辊弯辊力，IRS表示中间辊横移，ΔS为中间辊横移量。首先将轧辊简化为梁并离散为若干单元，然后根据轧辊参数建立各轧辊的有限元模型并对刚度矩阵和质量矩阵组集；包括：

步骤2.1.1：将六辊冷轧机中的各个轧辊简化为梁并离散为若干个单元，辊身以5mm为一个单元，辊颈等分为10份；

步骤2.1.2：根据插值函数确定各个单元的节点位移矢量的形函数表达式；确定梁单元位移矢量和插值函数并推导形函数表达式，考虑到轧辊应视为短粗梁且需要考虑剪切变形的影响，因此选择Timoshenko梁，由于本研究中仅需考虑轧辊垂向变形，因此节点单元位移矢量仅包含垂向位移及其对应转角即可，如图3所示，则梁单元的节点位移矢量表示为：

δ^e＝[v₁ θ₁ v₂ θ₂]^T

其中，δ^e为节点位移矢量；v₁为节点1的垂向位移；θ₁为节点1的转角；v₂为节点2的垂向位移；θ₂为节点2的转角；

为了提高计算精度，对节点位移矢量采用Hermite插值：

v(x)＝a₀+a₁x+a₂x²+a₃x³

θ(x)＝a₁+2a₂x+3a₃x²-γ

其中，a₀、a₁、a₂、a₃为待定系数；

由材料力学关系可得：

其中，M为弯矩；Q为剪切力；E为轧辊弹性模量；I为轧辊截面惯性矩；G为轧辊剪切模量；S为截面面积；γ为剪切应变；φ为考虑实际剪切应变和剪切应力不是均匀分布而引入的校正因子；

因此可以求得：

则θ(x)的表达式可以简化为：

θ(x)＝a₁+2a₂x+(3x²+6Λ)a₃

将边界条件v₁＝v(0)，θ₁＝θ(0)，v₂＝v(l)，θ₂＝θ(l)代入节点位移矢量表达式，即可得到节点位移矢量的形函数表达式：

其中，为垂向位移的形函数；/>为转角的形函数；l为梁单元长度；

步骤2.1.3：利用虚功原理推导单元刚度矩阵并进行组集；

根据虚功原理，单元应变能U^e由弯曲应变能U_b和剪切应变能U_s组成，表达式如下：

其中，则：

其中，和/>分别为单元弯曲刚度矩阵和单元剪切刚度矩阵；

步骤2.1.4：利用虚功原理推导单元质量矩阵并进行组集；

根据虚功原理，单元动能T^e由平动动能T_t和转动动能T_r组成，表达式如下：

其中，则：

其中，M^e _t和M^e _r分别为单元平动质量矩阵和单元转动质量矩阵；ρ为轧辊密度；I为轧辊截面惯性矩；S为截面面积；

步骤2.2：计算轧辊辊间接触刚度；如图4所示，

根据Hertz接触理论，两轧辊压缩时中心距离的变化量d_c为：

整理后可得：

两侧同时对p求导：

则，单位长度的轧辊辊间刚度系数可用如下式表示：

其中，中的1和2分别代指支撑辊与中间辊之间的辊间接触刚度系数和支撑辊与中间辊之间的辊间接触刚度系数，i和j为行号和列号；p为单位长度轧制力；E为轧辊弹性模量；v为轧辊泊松比；D₁和D₂为两轧辊辊径；

步骤2.3：根据轧辊的有限元模型和辊间接触刚度，建立六辊冷轧机的整体刚度、质量和阻尼矩阵；包括：

步骤2.3.1：利用压靠实验得到轧机刚度系数，采用有限元分析软件计算牌坊刚度系数；

步骤2.3.2：结合计算得到的辊间接触刚度确定支撑辊刚度系数，最后得到轧机整体刚度矩阵；

本实施例中，首先通过对某1450mm UCM冷轧机进行压靠实验，得到轧机总刚度系数为K＝4.4×10⁹N/m。然后利用有限元分析软件对牌坊进行建模并求解其垂直于轧制方向的刚度系数K_s＝4.1×10¹⁰N/m。

在确定轧机总刚度和牌坊刚度系数后，支撑辊刚度系数K_b可根据Hooke定律计算：

其中，K_w和K_im分别为工作辊和中间辊刚度系数，可由步骤2.2沿辊身方向累计求和得到；

由轧机的受力情况可知，轧机整体刚度矩阵K_z为：

K_z＝K+K_iw+K_bi

其中，K为支撑辊、中间辊和工作辊梁单元模型组合而成的总刚度矩阵；K_iw为中间辊和工作辊之间的辊间接触刚度矩阵；K_bi为支撑辊和中间辊之间的辊间接触刚度矩阵：

其中，n＝n_w+n_i+n_b；n₁～n₆为行和列的序号，根据以下条件选择：

其中，n_w、n_i和n_b分别为工作辊、中间辊和支撑辊的节点数目；n_s为中辊间横移量对应的节点数目。

步骤2.3.3：得到轧机的整体质量矩阵，利用Rayleigh阻尼公式、轧机整体刚度矩阵和质量矩阵获得轧机整体阻尼矩阵；

根据Rayleigh阻尼的构造方法，其假设结构的阻尼矩阵C是质量矩阵M和刚度矩阵K_z的线性组合，即：

C＝β₁M+β₂K_z

其中，M为支撑辊、中间辊和工作辊梁单元模型组合而成的整体质量矩阵；ω₁和ω₂为频率，取值分别为100Hz和500Hz；ξ₁和ξ₂为阻尼比，取值为0.03。

步骤3：将变形区沿轧制方向和宽度方向离散为若干份，将变形区如图5所示，图中Lubricate表示润滑油，σ_f表示前张力，σ_b表示后张力。沿轧制方向离散为400份、沿宽度方向离散为200份，计算动态轧制力；包括：

步骤3.1：根据带钢出入口厚度、轧辊辊径和轧辊垂向振动速度计算得到变形区动态接触弧长度；

考虑轧辊垂向振动会造成变形区变化，得到动态接触弧长度l_d的计算公式如下：

其中，R为轧辊压扁半径，单位mm；y_in和y_out为带钢入口和出口厚度，单位mm；θ为咬入角变化量，单位rad；v_y为轧辊垂向振动速度，速度向上为正，单位m/s。

轧辊压扁半径根据Hitchcock公式计算，如下式所示：

其中，R₀为轧辊初始半径；E_w为工作辊弹性模量，取值为2.1×10¹¹Pa；p为单位长度上的轧制力。

步骤3.2：根据带钢材质和微元体厚度利用变形抗力模型计算得到各微元体的平均变形抗力；

本实施例中，带钢的牌号为Q195，平均变形抗力σ_s采用如下公式计算：

其中，ε_Σ为累计变形量，y₀为热轧来料厚度；/>为道次平均厚度，各系数的经验取值分别为A＝498MPa、B＝136MPa、C＝0.2、D＝5。

步骤3.3：计算变形区摩擦应力分布；

由于高速轧制时，轧制界面处于边界润滑和流体润滑兼有的混合润滑状态，因此变形区摩擦应力分布τ可表示为：

其中，τ为变形区总摩擦应力分布；τ_a为粗糙接触产生的摩擦应力；τ_f为流体润滑产生的摩擦应力；k_s为材料的剪切强度；η为润滑油黏度；v_r为轧制速度；v_s为带钢沿轧制方向速度分布；h_t为平均油膜厚度；A_c为实际接触面积比。

由粗糙度分布假设可知，实际接触面积比A_c和平均油膜厚度h_t可分别表示为：

其中，Z＝h/3R_q为无量纲参数；f(δ)为概率密度函数，可表示为：

其中，δ为粗糙度分布，单位μm；R_q为带钢和轧辊的综合表面粗糙度，单位μm。

变形区油膜厚度分布h(x)可以用如下式表示：

其中，v_in为带钢入口速度，单位m/s；h₀为入口油膜厚度，单位mm，可由下式确定：

其中，η₀为大气压下的润滑油黏度；α为Barus公式的黏压系数，单位MPa^-1；l₀为v_y＝0时的接触弧长度。

步骤3.4：建立变形区微元体的力平衡微分方程；

结合变形区微元体的力平衡微分方程及变形区摩擦应力分布计算动态轧制力；

根据微元体的受力情况，力平衡微分方程如下所示：

其中，K_p＝1.155σ_s为材料的平面变形抗力，单位MPa；“+”为后滑区，“-”为前滑区；

步骤3.5：将步骤3.3所得的摩擦应力分布代入力平衡微分方程，对其沿轧制方向和宽度方向积分得到动态轧制力。

步骤4：根据轧辊和轧件间受力关系以及计算得到的动态轧制力，建立轧机-轧件系统的垂向振动动力学方程，求解得到轧辊位移响应用于预测轧机临界振动速度；包括：

步骤4.1：通过对轧件、轧辊和牌坊间的受力分析，建立轧机-轧件系统的垂向振动动力学方程：

其中，x、和/>分别为轧辊梁单元节点的位移、速度和加速度矢量；M为轧机整体质量矩阵；C为轧机整体阻尼矩阵；K_z为轧机整体刚度矩阵；K为支撑辊、中间辊和工作辊梁单元模型组合而成的总刚度矩阵；K_iw为中间辊和工作辊之间的辊间接触刚度矩阵；K_bi为支撑辊和中间辊之间的辊间接触刚度矩阵；F_iw为工作辊和中间辊间的接触压力分布，F_iw＝K_iwx；F_bi为中间辊和支撑辊间的接触压力分布，F_bi＝K_bix；F为轧制力宽向分布；β₁和β₂为比例系数；ξ₁和ξ₂为阻尼比；ω₁和ω₂为固有频率；

步骤4.2：对动力学方程采用Newmark-Beta法进行求解，可得到特定速度下的轧辊位移响应曲线。在设定的最大计算时间t_max下，当t≤t_max时继续下一时刻的动力学计算，当t>t_max时动力学计算终止，位移响应曲线的幅值恒定时对应的速度即为轧机的临界振动速度。

由于实际生产中，轧机振动频发在第四和五机架，因此本实施例对第四机架的临界振动速度进行预测。图6、图7和图8分别为工作辊弯辊力、中间辊横移量和中间辊弯辊力对轧机临界振动速度的影响以及对应临界振动速度下的带钢出口形貌。当轧制速度高于临界振动速度时，轧机发生振动；当轧制速度低于临界振动速度时，轧机处于稳定轧制状态。

本发明首先提供一种基于三维六辊冷轧机模型的临界振动速度预测方法，通过改进的Timoshenko梁构建三维六辊冷轧机模型并耦合动态轧制力计算模型，提高了计算速度和精度。其次，本发明不仅可以研究二维轧制工艺参数对轧制过程稳定性的影响，还可以分析弯蹿辊等宽向参数对临界轧制速度的影响，与其他方法相比分析参数更加全面合理。再次，采用本发明提供的方法可以在轧制规程制定阶段就预测出轧机的临界振动速度，为工艺参数优化提供理论支撑。最后，本发明基于振动机理与模拟仿真，可避免轧机振动导致的设备损坏或安全事故。

Claims (5)

1.一种基于三维模型的六辊冷轧机临界振动速度预测方法，其特征在于，包括：

步骤1：获取冷连轧产线实际生产中一段时间内的生产过程参数；

步骤2：根据六辊冷轧机中各个轧辊的有限元模型和辊间接触刚度，建立六辊冷轧机的整体刚度、质量和阻尼矩阵；

步骤3：将变形区沿轧制方向和宽度方向离散为若干份，计算动态轧制力；

步骤3.1：根据带钢出入口厚度、轧辊辊径和轧辊垂向振动速度计算得到变形区动态接触弧长l_d：

其中，R为轧辊压扁半径；y_in和y_out为带钢入口和出口厚度；θ为咬入角变化量；v_y为轧辊垂向振动速度，速度向上为正；

步骤3.2：根据带钢材质和微元体厚度利用变形抗力模型计算得到各微元体的平均变形抗力；

步骤3.3：计算变形区摩擦应力分布；

步骤3.4：建立变形区微元体的力平衡微分方程；

步骤3.5：将步骤3.3所得的摩擦应力分布代入力平衡微分方程，对其沿轧制方向和宽度方向积分得到动态轧制力；

步骤4：根据轧辊和轧件间受力关系以及计算得到的动态轧制力，建立轧机-轧件系统的垂向振动动力学方程，求解得到轧辊位移响应用于预测轧机临界振动速度；

步骤4.1：建立轧机-轧件系统的垂向振动动力学方程：

其中，x、和/>分别为轧辊梁单元节点的位移、速度和加速度矢量；M为轧机整体质量矩阵；C为轧机整体阻尼矩阵；K_z为轧机整体刚度矩阵；K为支撑辊、中间辊和工作辊梁单元模型组合而成的总刚度矩阵；K_iw为中间辊和工作辊之间的辊间接触刚度矩阵；K_bi为支撑辊和中间辊之间的辊间接触刚度矩阵；F_iw为工作辊和中间辊间的接触压力分布，F_iw＝K_iwx；F_bi为中间辊和支撑辊间的接触压力分布，F_bi＝K_bix；F为轧制力宽向分布；β₁和β₂为比例系数；ξ₁和ξ₂为阻尼比；ω₁和ω₂为频率；

步骤4.2：对动力学方程采用Newmark-Beta法进行求解，可得到特定速度下的轧辊位移响应曲线，位移响应曲线的幅值恒定时对应的速度即为轧机的临界振动速度。

2.根据权利要求1所述的一种基于三维模型的六辊冷轧机临界振动速度预测方法，其特征在于，所述生产过程参数包括轧机结构参数、轧制工艺参数、带钢参数以及润滑油参数；

所述轧机结构参数包括：轧辊材质、弹性模量、泊松比、密度、辊身长度和直径、辊颈长度和直径以及牌坊材质、质量和尺寸；

所述轧制工艺参数包括：机架间前后张力、各道次轧制速度、各道次带钢入口速度、各道次带钢出入口厚度、工作辊弯辊力以及中间辊弯辊力和横移量；

所述带钢参数以及润滑油参数包括：带钢牌号和宽度、来料厚度以及润滑油黏度和黏压系数。

3.根据权利要求1所述的一种基于三维模型的六辊冷轧机临界振动速度预测方法，其特征在于，所述步骤2包括：

步骤2.1：建立六辊冷轧机中的各个轧辊的有限元模型；

步骤2.2：计算轧辊辊间接触刚度；

步骤2.3：根据轧辊的有限元模型和辊间接触刚度，建立六辊冷轧机的整体刚度、质量和阻尼矩阵。

4.根据权利要求3所述的一种基于三维模型的六辊冷轧机临界振动速度预测方法，其特征在于，所述步骤2.1包括：

步骤2.1.1：将六辊冷轧机中的各个轧辊简化为梁并离散为若干个单元；

步骤2.1.2：根据插值函数确定各个单元的节点位移矢量的形函数表达式；

步骤2.1.3：利用虚功原理推导单元刚度矩阵并进行组集：

步骤2.1.4：利用虚功原理推导单元质量矩阵并进行组集。

5.根据权利要求3所述的一种基于三维模型的六辊冷轧机临界振动速度预测方法，其特征在于，所述步骤2.3包括：

步骤2.3.1：利用压靠实验得到轧机刚度系数，采用有限元分析软件计算牌坊刚度系数；

步骤2.3.2：结合计算得到的辊间接触刚度确定支撑辊刚度系数，最后得到轧机整体刚度矩阵：

K_z＝K+K_iw+K_bi

其中，K_z为轧机整体刚度矩阵；K为支撑辊、中间辊和工作辊梁单元模型组合而成的总刚度矩阵；K_iw为中间辊和工作辊之间的辊间接触刚度矩阵；K_bi为支撑辊和中间辊之间的辊间接触刚度矩阵：

其中，n＝n_w+n_i+n_b，n₁～n₆为行和列的序号，根据以下条件选择：

其中，n_w、n_i和n_b分别为工作辊、中间辊和支撑辊的节点数目；n_s为中辊间横移量对应的节点数目；

步骤2.3.3：利用Rayleigh阻尼公式、轧机整体刚度矩阵和质量矩阵获得轧机整体阻尼矩阵。