轧制力是进行轧机强度校核与轧制工艺设定及优化的依据。构建预测精准的轧制力模型对于提高板材的尺寸精度和质量有重要的意义。其中，速度场的设定是进行轧制过程轧件运动及变形分析的基础，在轧制力模型的构建中处于基础性地位。关于轧制速度场的研究可以追溯到1975年Oh和Kobayashi^[1]的工作。他们提出了流函数速度场，并基于该速度场建立了一个轧制力模型。考虑到速度场设定对轧制力模型精度的影响，Kato等^[2]随后在1980年针对棒材轧制过程中不均匀变形的特点提出了一个加权速度场。这两种速度场的提出很好地促进了轧制力解析建模的发展。

数值模拟是求解轧制力的另一种有效方法。它能够准确地模拟出轧制过程中的力能参数及其变化规律。Kobayashi等^[3]在1989年提出使用有限元法模拟金属变形的设想，但碍于当时的计算机水平，未能给出相应的数值结果。Hwang等^[4]于1992年提出了一种用于分析热轧带钢的刚塑性有限元法，并在板材的模拟中获得了成功。Heislitz等^[5]使用了有限元法对轧制的应力分布以及成品的几何形状进行了模拟，并用实验结果对模拟结果的正确性进行了验证。Mei等^[6-7]为了提高有限元模拟的精度并减少迭代次数，使用了多种函数对轧制初始速度场进行了对比优化，提出了设定速度场的GF-RM法。王宝明等^[8-9]使用有限元模拟软件MSC.MARC对全浮动芯棒连轧过程进行了模拟。结果表明，开轧温度和压下量是调整轧制工艺时需要考虑的两个主要因素，而张力变化对连轧机组轧制力的影响主要在第二机架上。孙建亮等^[10]为了解决轧制力偏差过大的问题，使用有限元软件ANSYS/LS-DANA对厚板轧制进行了模拟。通过与现场实测数据作对比，确定了轧辊交叉、非对称弯辊和窜辊是造成轧机产生轧制力偏差的主要原因，并提出了降低偏差的相应措施。彭林等^[11]为了研究非线性塑性大变形过程，使用有限元软件对H型钢进行了模拟研究。模拟结果直观地体现出了金属流动的演化过程，并且揭示了造成轧制后出现“舌头”缺陷的原因。尽管数值模拟能够明确模拟出轧制过程中各个参数的具体数值，但所需的计算量仍然较大。另一方面，数值模拟得到的结果是针对某一具体工艺的离散数值解，不能反映轧制过程中各个参数间的函数制约关系，不利于进行现场的工艺设计与优化。

理论解析作为另一种常用的轧制力求解方法弥补了这一不足。它能明确给出各个参数之间的关系，且具有物理意义明确的特点。Sims^[12]较早基于工程法对轧制过程进行了一定的简化，得到了轧制力与轧制力矩的解析解。Alexander等^[13]则使用了滑移线法描述了金属轧制变形过程，并建立了相应的模型。Moon等^[14]则考虑了轧制过程中几何因子与形状参数的影响，建立了一个预测轧制力的近似模型。Pan等^[15]对传统的三角形速度场进行了分析。他们发现，在多三角形速度场中，可以通过描述三角形的数量来预估变形功率的上下限。赵德文等^[16]基于柱坐标系速度场对孔型轧制进行了理论解析。结果证明，初始速度场设定的合理性将会对后续的轧制力预测精度有较大影响。近些年来不少研究者们也对轧制速度场的设定进行了相关研究，如双曲正弦速度场^[17]、抛物线速度场^[18]以及整体加权速度场^[19]等。这些研究为提高轧制力模型的精度做出了许多贡献。但是，这些速度场还不能较好地反应厚板轧制变形的特点。实际应用中，已有模型由于作了较多简化，因此预测精度有限。另一方面，已有模型多对非线性Mises比塑性功率进行近似替代处理，也不可避免地带来了一定的计算误差。

针对以上问题，根据轧件变形时金属流动的特点提出一个余弦速度场，并基于该速度场进行相应的能量分析。提出使用分矢量内积加和法对Mises非线性轧制功率进行等价转化，解决了非线性轧制功率积分困难的问题。同时，基于提出的速度场得到了摩擦功率和剪切功率表达式，并基于刚塑性变分原理导出了轧制力的解析模型。最后，将得到的轧制力理论结果与实验数据进行了对比，验证了该轧制力模型的精度，并分析了各轧制参数与轧制力之间的关系。

1 轧制参数描述 1.1 变形参数

考虑到变形区的对称性，只取整体变形区的1/4进行分析，见图 1。

图中，h₀为轧件入口厚度，h₁为出口厚度，R为轧辊半径，O点位于轧件变形区入口位置，v₀为入口速度，v₁为出口速度。θ为接触角，α为中性角，x为变形区任意处到变形区入口的水平距离，h_x为变形区内任意处轧件厚度。中性点处轧件切向速度与此处轧辊的切向速度相等，切向速度不连续量与摩擦功均为零。根据图中的几何关系，接触弧方程、参数方程以及其一阶、二阶导数可表示为

$ \left\{\begin{array}{l} z=h_{x}=R+h_{1}-\left[R^{2}-(l-x)^{2}\right]^{1 / 2} \\ z=h_{\alpha}=R+h_{1}-R \cos \alpha \\ l-x=R \sin \alpha, \mathrm{d} x=-R \cos \alpha \mathrm{d} \alpha \\ {h^{\prime}}_{x}=-\tan \alpha, {h^{\prime \prime}}_{x}=\left(R \cos ^{3} \alpha\right)^{-1} \end{array}\right. $

(1)

由图(1)可知边界条件

$ \left\{\begin{array}{l} x=0, \alpha=\theta ; h_{x}=h_{\alpha}=h_{\theta}=h_{0}, {h^{\prime}}_{x}=-\tan \theta \\ x=l, \alpha=0 ; h_{x}=h_{\alpha}=h_{1}, {h^{\prime}}_{x}=0 \end{array}\right. $

(2)

轧件在变形过程中不仅存在着垂直方向的变形，在宽度方向同样存在着变形。但是轧件的宽厚比远大于10，宽展量很小^[20]。图 2为轧件的宽展俯视图，其中的宽展假定为抛物线形式，其数学表达式、一阶导数以及平均宽度表达见式(3)~(5)。

$ b_{x}=b_{1}-\frac{\Delta b}{l^{2}}(l-x)^{2} $

(3)

$ {b^{\prime}}_{x}=\frac{2 \Delta b}{l^{2}}(l-x) $

(4)

$ \bar{b}_{x}=\frac{b_{0}+2 b_{1}}{3} $

(5)

1.2 余弦速度场

本文假定轧制时从入口到出口处发生的金属流动符合余弦函数变化规律，写作v_x=a－bcos cx的形式，其表达式以及一、二阶导数为：

$ v_{x}=\frac{v_{0}}{2}\left[(1+r)-(r-1) \cos \left(\frac{{\rm{ \mathsf{ π} }}}{l} x\right)\right] $

(6)

$ {v^{\prime}}_{x}=\frac{v_{0} {\rm{ \mathsf{ π} }}}{2 l}(r-1) \sin \left(\frac{{\rm{ \mathsf{ π} }}}{l} x\right) $

(7)

$ {v^{\prime \prime}}_{x}=\frac{v_{0} {\rm{ \mathsf{ π} }}^{2}}{2 l^{2}}(r-1) \cos \left(\frac{{\rm{ \mathsf{ π} }}}{l} x\right) $

(8)

式中r=(h₀b₀)/(h₁b₁)称为截面压缩比。

根据速度协调条件，可以导出v_y、v_z如下：

$ v_{y}=v_{x} \frac{{h^{\prime}}_{x}}{h_{x}} y $

(9)

$ v_{z}=v_{x} \frac{{h^{\prime}}_{x}}{h_{x}} z $

(10)

于是，可得到本文的余弦速度场：

$ \left\{\begin{array}{l} v_{x}=\frac{v_{0}}{2}\left[(1+r)-(r-1) \cos \left(\frac{{\rm{ \mathsf{ π} }}}{l} x\right)\right] \\ v_{y}=v_{x} \frac{{h^{\prime}}_{x}}{h_{x}} y \\ v_{z}=v_{x} \frac{{h^{\prime}}_{x}}{h_{x}} z \end{array}\right. $

(11)

式中v_x、v_y、v_z分别是轧制方向、宽度方向以及压下方向的速度分量。根据几何方程^[21]，速度场对应的应变速率场为：

$ \left\{\begin{array}{l} \dot{\varepsilon}_{x}=-v_{x}\left(\frac{{b^{\prime}}_{x}}{b_{x}}+\frac{{h^{\prime}}_{x}}{h_{x}}\right) \\ \dot{\varepsilon}_{y}=v_{x} \frac{{b^{\prime}}_{x}}{b_{x}} \\ \dot{\varepsilon}_{z}=v_{x} \frac{{h^{\prime}}_{x}}{h_{x}} \\ \dot{\varepsilon}_{x y}=-\frac{y}{2}\left({v^{\prime \prime}}_{x}+\frac{{v^{\prime}}_{x} {h^{\prime}}_{x}+v_{x} {h^{\prime \prime}}_{x}}{h_{x}}-\frac{v_{x} h_{x}^{\prime 2}}{h_{x}^{2}}\right) \\ \dot{\varepsilon}_{x z}=\frac{z}{2 h_{x}}\left({v^{\prime}}_{x} {h^{\prime}}_{x}+v_{x} {h^{\prime \prime}}_{x}-\frac{v_{x} h_{x}^{\prime 2}}{h_{x}}\right) \\ \dot{\varepsilon}_{y z}=0 \end{array}\right. $

(12)

该速度场能够严格满足体积不变条件、出入口速度边界条件和几何方程，因而可以充分满足运动许可条件。以下部分将以该速度场为基础进行轧制力的能量分析。

2 轧制力的能量解析 2.1 内部变形功率

内部变形功率N_d可由变形材料的等效应力和等效应变速率确定，其计算式为

$ \begin{gathered} N_{\mathrm{d}}=\iiint\limits_{V} \bar{\sigma} \bar{\varepsilon} \mathrm{d} V= \\ 4 \iiint\limits_{V} \sigma_{\mathrm{s}} \sqrt{\frac{2}{3} \dot{\varepsilon}_{i j} \dot{\varepsilon}_{i j}} \mathrm{~d} V \end{gathered} $

(13)

式中: σ为轧件等效应力，$ \bar {\dot {\varepsilon}} $为轧件的等效应变速率，σ_s为变形抗力，ε_ij为应变速率分量。

由于式(13)基于的是非线性的Mises屈服准则，会给后续的积分计算带来困难，因此使用分矢量内积加和法进行计算，转化过程如下:

$ \begin{aligned} N_{\mathrm{d}}=& 4 \sqrt{\frac{2}{3}} \sigma_{\mathrm{s}} \int_{0}^{h_{x}} \int_{0}^{b_{x}} \int_{0}^{l} \dot{\boldsymbol{\varepsilon}} \cdot \dot{\boldsymbol{\varepsilon}}_{0} \mathrm{d} x \mathrm{d} y \mathrm{d} z=\\ & 4 \sqrt{\frac{2}{3}} \sigma_{\mathrm{s}}\left(\iiint\limits_{V} \eta_{x} l_{x} \mathrm{~d} V+\iiint\limits_{V} \eta_{y} l_{y} \mathrm{~d} V+\iiint\limits_{V} \eta_{z} l_{z} \mathrm{~d} V\right)=\\ & 4 \sqrt{\frac{2}{3}} \sigma_{\mathrm{s}}\left(I_{x}+I_{y}+I_{z}\right) \end{aligned} $

(14)

式中，$\mathit{\boldsymbol{\dot \varepsilon = }}{\eta _x}\mathit{\boldsymbol{i}} + {\eta _y}\mathit{\boldsymbol{j}} + {\eta _z}\mathit{\boldsymbol{k}}$为轧件变形过程中的应变速率矢量，${\mathit{\boldsymbol{\dot \varepsilon }}_0}\mathit{\boldsymbol{ = }}{l_x}\mathit{\boldsymbol{i}} + {l_y}\mathit{\boldsymbol{j}} + {l_z}\mathit{\boldsymbol{k}}$为单位向量，l_x=cos α，l_y=cos β，l_z=cos γ分别为单位向量在各个轴上夹角的余弦值。

以下为轧制过程中各个方向的矢量模长以及它们与主轴的夹角余弦。

$ \left\{\begin{array}{l} \eta_{x}=\left|\dot{\boldsymbol{\varepsilon}}_{x}\right|=v_{x}\left(\frac{{b^{\prime}}_{x}}{b_{x}}+\frac{{h^{\prime}}_{x}}{h_{x}}\right) \\ \eta_{y}=\left|\dot{\boldsymbol{\varepsilon}}_{y}\right|+\left|\dot{\boldsymbol{\varepsilon}}_{x y}\right|=v_{x} \frac{{b^{\prime}}_{x}}{b_{x}}+\frac{y}{\sqrt{2}} v_{x}\left[\frac{{b^{\prime \prime}}_{x}}{b_{x}}-\frac{{h^{\prime}}_{x}}{h_{x}} \cdot \frac{{b^{\prime}}_{x}}{b_{x}}-2\left(\frac{{b^{\prime}}_{x}}{b_{x}}\right)^{2}\right] \\ \eta_{z}=\left|\dot{\boldsymbol{\varepsilon}}_{z}\right|+\left|\dot{\boldsymbol{\varepsilon}}_{x z}\right|=v_{x} \frac{{h^{\prime}}_{x}}{h_{x}}+\frac{z}{\sqrt{2}} v_{x}\left[\frac{{h^{\prime \prime}}_{x}}{h_{x}}-\frac{{h^{\prime}}_{x}}{h_{x}} \cdot \frac{{b^{\prime}}_{x}}{b_{x}}-2\left(\frac{{h^{\prime}}_{x}}{h_{x}}\right)^{2}\right] \end{array}\right. $

(15)

$ \left\{\begin{array}{l} l_{x}=\cos \alpha=\sqrt{1+\left(\frac{\mathrm{d} y}{\mathrm{~d} x}\right)^{2}+\left(\frac{\mathrm{d} z}{\mathrm{~d} x}\right)^{2}} \\ l_{y}=\cos \beta=\sqrt{1+\left(\frac{\mathrm{d} x}{\mathrm{~d} y}\right)^{2}+\left(\frac{\mathrm{d} z}{\mathrm{~d} y}\right)^{2}} \\ l_{z}=\cos \gamma=\sqrt{1+\left(\frac{\mathrm{d} x}{\mathrm{~d} z}\right)^{2}+\left(\frac{\mathrm{d} y}{\mathrm{~d} z}\right)^{2}} \end{array}\right. $

(16)

注意到式(15)为x的单值函数，使用积分中值定理后可得

$ \left\{\begin{array}{l} \frac{\bar{h}_{x}^{\prime}}{h_{x}}=\frac{1}{l} \int_{0}^{l} \frac{h^{\prime}{ }_{x}}{h_{x}} \mathrm{~d} x=-\frac{\varepsilon_{3}}{l} ; \frac{\bar{b}_{x}^{\prime \prime}}{b_{x}}=0 \\ \overline{h_{x}}=\frac{h_{0}+2 h_{1}}{3} ; \overline{h}^{\prime \prime}_{x}=\frac{1}{l} \int_{0}^{l} h^{\prime \prime}{ }_{x} \mathrm{~d} x=\frac{1}{R-\Delta h} \\ \frac{\overline{h}^{\prime \prime}_{x}}{\overline{h}_{x}}=\frac{3 \tan \theta}{l\left(h_{0}+2 h_{1}\right)}=\frac{3}{(R-\Delta h)\left(h_{0}+2 h_{1}\right)} \\ \bar{v}_{x}=\frac{1}{l} \int_{0}^{l} v_{x} \mathrm{~d} x=v_{0}+\frac{\Delta v}{2} ; \bar{v}_{x}^{\prime \prime}=\frac{1}{l} \int_{0}^{l} {v^{\prime \prime}}_{x} \mathrm{d} x=0 \end{array}\right. $

(17)

其中: l为接触弧的水平投影长度，$l = \sqrt {2R\Delta h} $, ${\bar h_x}$、${\bar b_x}$、${\bar v_x}$分别为按积分中值定理计算所得的变形区平均厚度、平均宽度和轧件沿x轴方向平均速度；ε₂=ln(b₁/b₀)和ε₃=ln(h₀/h₁)为轧件宽、厚度方向的主应变。

根据变形区几何关系可以确定：

$ \left\{\begin{array}{l} \frac{\mathrm{d} y}{\mathrm{~d} x}={b^{\prime}}_{x} \approx \frac{\Delta b}{l}=\bar{b}_{x}^{\prime} \\ \frac{\mathrm{d} z}{\mathrm{~d} x}={h^{\prime}}_{x}=-\frac{\Delta h}{l}=\bar{h}_{x}^{\prime} \\ \frac{\mathrm{d} z}{\mathrm{~d} y}=0 \end{array}\right. $

(18)

由此可将各方向的应变速率与主轴的余弦值表示成

$ \left\{\begin{array}{l} l_{x}=\frac{l}{\sqrt{l^{2}+\Delta b^{2}+\Delta h^{2}}} \\ l_{y}=\frac{\Delta b}{\sqrt{\Delta b^{2}+l^{2}}} \\ l_{z}=\frac{\Delta h}{\sqrt{\Delta h+l^{2}}} \end{array}\right. $

(19)

将式(15)~(17)以及积分中值得到的各参数简化关系代入式(14)后可得到:

$ \begin{aligned} I_{x}=& \int_{0}^{h_{x}} \int_{0}^{b_{x}} \int_{0}^{l} \frac{v_{x}\left(\frac{{b^{\prime}}_{x}}{b_{x}}+\frac{{h^{\prime}}_{x}}{h_{x}}\right)}{\sqrt{1+\left({b^{\prime}}_{x}\right)^{2}+\left({h^{\prime}}_{x}\right)^{2}}} \mathrm{d} x \mathrm{d} y \mathrm{d} z=\\ & \frac{l U \ln \frac{b_{1} h_{0}}{b_{0} h_{1}}}{\sqrt{l^{2}+\Delta b^{2}+\Delta h^{2}}}=U f_{1} \end{aligned} $

(20)

$ \begin{aligned} I_{y} =&\int_{0}^{h_{x}} \int_{0}^{b_{x}} \int_{0}^{l} \frac{v_{x} \frac{{b^{\prime}}_{x}}{b_{x}}+\frac{y}{\sqrt{2}} v_{x}\left[\frac{{b^{\prime \prime}}_{x}}{b_{x}}-\frac{{h^{\prime}}_{x}}{h_{x}} \cdot \frac{{b^{\prime}}_{x}}{b_{x}}-2\left(\frac{{b^{\prime}}_{x}}{b_{x}}\right)^{2}\right]}{\sqrt{1+\left(1 / {b^{\prime}}_{x}\right)^{2}}} \mathrm{d} x \mathrm{d} y \mathrm{d} z=\\ & \frac{\Delta b}{\sqrt{\Delta b^{2}+l^{2}}} \int_{0}^{h_{x}} \int_{0}^{b_{x}} \int_{0}^{l} v_{x} \frac{{b^{\prime}}_{x}}{b_{x}}+\\ & \frac{y}{\sqrt{2}} v_{x}\left[\frac{{b^{\prime \prime}}_{x}}{b_{x}}-\frac{{h^{\prime}}_{x}}{h_{x}} \cdot \frac{{b^{\prime}}_{x}}{b_{x}}-2\left(\frac{{b^{\prime}}_{x}}{b_{x}}\right)^{2}\right] \mathrm{d} x \mathrm{d} y \mathrm{d} z=\\ & \frac{\Delta b U}{\sqrt{\Delta b^{2}+l^{2}}}\left[\ln \frac{b_{1}}{b_{0}}+\frac{\left(\varepsilon_{2} \varepsilon_{3}-2 \varepsilon_{3}^{2}\right) \bar{b}_{x}}{2 \sqrt{2} l}\right]=U f_{2} \end{aligned} $

(21)

注意到y与z无关，dx/dz=1/h′_x=－l/2Δh，有

$ \begin{aligned} I_{z}=& \int_{0}^{h_{x}} \int_{0}^{b_{x}} \int_{0}^{l}\frac{v_{x} \frac{{h^{\prime}}_{x}}{h_{x}}+\frac{z}{\sqrt{2}} v_{x}\left[\frac{{h^{\prime \prime}}_{x}}{h_{x}}-\frac{{h^{\prime}}_{x}}{h_{x}} \cdot \frac{{b^{\prime}}_{x}}{b_{x}}-2\left(\frac{{h^{\prime}}_{x}}{h_{x}}\right)^{2}\right]}{\sqrt{1+\left(1 / {h^{\prime}}_{x}\right)^{2}}} \mathrm{d} x \mathrm{d} y \mathrm{d} z=\\ & \frac{2 \Delta h U}{\sqrt{4 \Delta h^{2}+l^{2}}}\left\{\begin{array}{l} \ln \frac{h_{0}}{h_{1}}+\frac{1}{2 \sqrt{2} l}\left[\frac{l^{2}}{(R-\Delta h)}+\right. \\ \left.\left(\varepsilon_{2} \varepsilon_{3}-2 \varepsilon_{3}^{2}\right) \bar{h}_{x}\right] \end{array}\right\}=U f_{3} \end{aligned} $

(22)

将逐项积分结果I_x、I_y、I_z代入式(14)后整理得

$\begin{aligned} &N_{\mathrm{d}}=4 \sqrt{\frac{2}{3}} \sigma_{\mathrm{s}}\left(U f_{1}+U f_{2}+U f_{3}\right)=\\ &4 \sqrt{\frac{2}{3}} \sigma_{\mathrm{s}} U \left(\begin{aligned} &\frac{l\ln \frac{b_{1} h_{0}}{b_{0} h_{1}}}{\sqrt{l^{2}+\Delta b^{2}+\Delta h_{1}^{2}}}+\frac{\Delta b}{\sqrt{\Delta b^{2}+l^{2}}}\left[\begin{aligned}&\ln \frac{b_{1}}{b_{0}}+\\ &\frac{\left(\varepsilon_{2} \varepsilon_{3}-2 \varepsilon_{3}^{2}\right) \bar{b}_{x}}{2 \sqrt{2} l}\end{aligned}\right]+ \\ &\frac{2 \Delta h}{\sqrt{4 \Delta h^{2}+l^{2}}}\left\{\ln \frac{h_{0}}{h_{1}}+\frac{1}{2 \sqrt{2} l}\left[\begin{array}{l} \frac{l^{2}}{(R-\Delta h)}+ \\ \left(\varepsilon_{2} \varepsilon_{3}-2 \varepsilon_{3}^{2}\right) \bar{h}_{x} \end{array}\right]\right\} \end{aligned}\right) \end{aligned} $

(23)

2.2 摩擦功率

摩擦功率N_f可由摩擦剪应力与速度不连续量的乘积求得，计算式为

$ N_{\mathrm{f}}=4 \int\limits_{S} \tau_{\mathrm{f}} \Delta v_{\mathrm{f}} \mathrm{d} S $

(24)

式中Δv_f为速度不连续量，τ_f为摩擦剪应力，其中Δv_f和dS的计算式为

$ \left\{\begin{array}{l} \Delta v_{\mathrm{f}}=v_{\mathrm{R}}-v_{x} \sqrt{1+h_{x}^{\prime 2}}=v_{\mathrm{R}}-v_{x} \sec \alpha \\ \mathrm{d} S=\sqrt{1+h^{\prime 2}} \mathrm{~d} x \mathrm{d} y=\sec \alpha \mathrm{d} x \mathrm{d} y \end{array}\right. $

(25)

同样，将轧件表面速度不连续量与摩擦剪应力表示为矢量的形式

$ \left\{\begin{array}{l} \Delta v_{\mathrm{f}}=\Delta v_{x} \boldsymbol{i}+\Delta v_{y} \boldsymbol{j}+\Delta v_{z} \boldsymbol{k} \\ \tau_{\mathrm{f}}=\tau_{\mathrm{f}_{x}} \boldsymbol{i}+\tau_{\mathrm{f}_{y}} \boldsymbol{j}+\tau_{\mathrm{f}_{z}} \boldsymbol{k} \end{array}\right. $

(26)

将上式代入式(24)后，即可把摩擦功率写成矢量内积的形式

$ \begin{aligned} &N_{\mathrm{f}}=4 \int_{0}^{l} \int_{0}^{b_{x}} \tau_{\mathrm{f}} \cdot \Delta v_{\mathrm{f}} \mathrm{d} S= \\ &4 \int_{0}^{l} \int_{0}^{b_{x}}\left(\tau_{\mathrm{f}_{x}} \Delta v_{x}+\tau_{\mathrm{f}_{y}} \Delta v_{y}+\tau_{\mathrm{f}_{z}} \Delta v_{z}\right) \sqrt{1+h_x^{\prime 2}} \mathrm{d} x \mathrm{d} y= \\ &4 m k \int_{0}^{l} \int_{0}^{b_{x}}\left(\Delta v_{x} \cos \alpha+\Delta v_{y} \cos \beta+\Delta v_{z} \cos \gamma\right) \sec \alpha \mathrm{d} x \mathrm{d} y= \\ &4 m k\left(\begin{array}{l} \int_{0}^{l} \int_{0}^{b_{x}} \Delta v_{x} \cos \alpha \sec \alpha \mathrm{d} x \mathrm{d} y+\int_{0}^{l} \int_{0}^{b_{x}} \Delta v_{y} \cos \beta \sec \alpha \mathrm{d} x \mathrm{d} y+ \\ \int_{0}^{l} \int_{0}^{b_{x}} \Delta v_{z} \cos \gamma \sec \alpha \mathrm{d} x \mathrm{d} y \end{array}\right) \end{aligned} $

(27)

根据图 1，轧辊表面各方向的余弦值如下

$ \begin{gathered} \cos \alpha=\pm \frac{\sqrt{R^{2}-(l-x)^{2}}}{R} ; \cos \beta=0 \\ \cos \gamma=\sin \alpha=\pm \frac{(l-x)}{R} \end{gathered} $

(28)

将式(28)代入式(27)后进行积分计算，并注意以中性点位置为分界点，中性点两侧的积分上下限要交换，分段进行积分后可得到摩擦功率的表达式

$ N_{\mathrm{f}}=4 m k\left[R v_{\mathrm{R}} b_{1}\left(\theta-2 \alpha_{n}\right)+\frac{R U}{h_{\mathrm{m}}} \ln \frac{\tan ^{2}\left(\frac{\alpha_{n}}{2}+\frac{{\rm{ \mathsf{ π} }}}{4}\right)}{\tan \left(\frac{\theta}{2}+\frac{{\rm{ \mathsf{ π} }}}{4}\right)}\right] $

(29)

2.3 剪切功率

轧制时在变形区入口与出口处的剪切功耗记为N_s0与N_s1，以下分别进行计算:

$ \begin{aligned} N_{\mathrm{s}_{0}}=& 4 \int k\left|\Delta \bar{v}_{t}\right| \mathrm{d} S=\\ & 4 k \int_{0}^{b_{0}} \int_{0}^{h_{0}}\left(\bar{v}_{y} \sqrt{1+\left(\bar{v}_{z} / \bar{v}_{y}\right)^{2}}\right) \mathrm{d} z \mathrm{d} y=\\ & \frac{k \tan \theta \Delta b U b_{0}}{b_{1} h_{0}} \sqrt{1+\frac{4 b_{1}^{2} h_{0}^{2}}{\Delta b^{2} b_{0}^{2}}} \end{aligned} $

(30)

式中存在着如下关系:

$ \left.\Delta \bar{v}_{t}\right|_{x=0}=\left.\sqrt{\Delta \bar{v}_{y}^{-2}+\Delta \bar{v}_{z}^{2}}\right|_{x=0}=\left.\bar{v}_{y} \sqrt{1+\left(\bar{v}_{z} / \bar{v}_{y}\right)^{2}}\right|_{x=0} $

(31)

$ \left.\bar{v}_{z}\right|_{x=0}=\frac{1}{h_{0}} \int_{0}^{h_{0}} v_{z} \mathrm{~d} z=-\frac{\tan \theta v_{0}}{2} $

(32)

$ \left.\bar{v}_{y}\right|_{x=0}=\frac{1}{b_{0}} \int_{0}^{b_{0}} \left.v_{y}\right|_{x=0} \mathrm{~d} y=\frac{\Delta b v_{0} b_{0} \tan \theta}{4 b_{1} h_{0}} $

(33)

在变形区出口处，由于h′_x=b′_x=0，因此出口处的金属流动存在着以下规律

$ \left.\bar{v}_{y}\right|_{x=l}=\left.\bar{v}_{z}\right|_{x=l}=\left.\Delta v_{y}\right|_{x=l}=\left.\Delta v_{z}\right|_{x=l}=0 $

(34)

因此，在轧制出口处剪切功率N_s1=0，入口处的剪切功率消耗为轧制过程中的总剪切功率，可表示为

$ N_{\mathrm{s}}=N_{\mathrm{s}_{0}}=\frac{k \tan \theta \Delta b U b_{0}}{b_{1} h_{0}} \sqrt{1+\frac{4 b_{1}^{2} h_{0}^{2}}{\Delta b^{2} b_{0}^{2}}} $

(35)

2.4 总功泛函最小化

根据式(23)、(29)和(35)可将轧制的总功率泛函Φ表达成

$ \begin{aligned} \varPhi=& N_{\mathrm{d}}+N_{\mathrm{f}}+N_{\mathrm{s}}=4 \sqrt{\frac{2}{3}} \sigma_{\mathrm{s}}\left(U f_{1}+U f_{2}+U f_{3}\right)+\\ & 4 m k\left[R v_{\mathrm{R}} b_{1}\left(\theta-2 \alpha_{n}\right)+\frac{R U}{h_{\mathrm{m}}} \ln \frac{\tan ^{2}\left(\frac{\alpha_{n}}{2}+\frac{{\rm{ \mathsf{ π} }}}{4}\right)}{\tan \left(\frac{\theta}{2}+\frac{{\rm{ \mathsf{ π} }}}{4}\right)}\right]+\\ & \frac{k \tan \theta \Delta b U b_{0}}{b_{1} h_{0}} \sqrt{1+\frac{4 b_{1}^{2} h_{0}^{2}}{\Delta b^{2} b_{0}^{2}}} \end{aligned} $

(36)

式中U为轧制变形区内任意变形区的秒体积流量。其表达式与一阶偏导式可写为:

$ \begin{aligned} U=&v_{x} h_{x} b_{x}=\left\{v_{0}+\frac{\Delta v}{2}-\frac{\Delta v}{2} \cos \left[\frac{{\rm{ \mathsf{ π} }}}{l}\left(l-R \sin \alpha_{n}\right)\right]\right\} \\ &\left(R+h_{1}-R \cos \alpha_{n}\right)\left\{b_{1}-\frac{\Delta b}{l^{2}}\left[l-\left(l-R \sin \alpha_{n}\right)\right]^{2}\right\} \end{aligned} $

(37)

$ \begin{aligned} \frac{\mathrm{d} U}{\mathrm{~d} \alpha_{n}}=&\left(R+h_{1}-R \cos \alpha_{n}\right) \\ &\left\langle\left\{-\frac{\Delta v {\rm{ \mathsf{ π} }}}{2 l} \cos \alpha_{n} \sin \left[\frac{{\rm{ \mathsf{ π} }}}{l}\left(l-R \sin \alpha_{n}\right)\right]\right\}\right.\\ &\left\{b_{1}-\frac{\Delta b}{l^{2}}\left[l-\left(l-R \sin \alpha_{n}\right)\right]^{2}\right\}-\\ &\left.\frac{2 \Delta b R}{l^{2}} \cos \alpha_{n}\left(l-R \sin \alpha_{n}\right)\right\rangle+\\ & R \sin \alpha_{n}\left\{b_{1}-\frac{\Delta b}{l^{2}}\left[l-\left(l-R \sin \alpha_{n}\right)\right]^{2}\right\}=N \end{aligned} $

(38)

式(23)、(29)和(35)分别对中性角α_n求导可得:

$ \frac{\mathrm{d} N_{\mathrm{d}}}{\mathrm{d} \alpha_{n}}=4 \sqrt{\frac{2}{3}} \sigma_{\mathrm{s}} N\left(f_{1}+f_{2}+f_{3}\right) $

(39)

$ \frac{\mathrm{d} N_{\mathrm{f}}}{\mathrm{d} \alpha_{n}}=4 m k\left\{\begin{array}{l} 4 v_{1}\left(\sec \alpha_{n}-\cos 2 \alpha_{n}\right)+v_{\mathrm{R}}\left(\cos 2 \alpha_{n}-1\right)+ \\ \frac{2 \Delta v R^{2}}{l^{2}}\left[\begin{array}{l} \cos \alpha_{n}\left(1+\cos \alpha_{n}-\frac{1}{\cos \alpha_{n}+1}\right)- \\ \tan \alpha_{n} \end{array}\right] \end{array}\right\} $

(40)

$ \frac{\mathrm{d} N_{\mathrm{s}}}{\mathrm{d} \alpha_{n}}=\frac{k \tan \theta \Delta b N b_{0}}{b_{1} h_{0}} \sqrt{1+\frac{4 b_{1}^{2} h_{0}^{2}}{\Delta b^{2} b_{0}^{2}}} $

(41)

对总功率泛函变分有

$ \frac{\mathrm{d} \varPhi}{\mathrm{d} \alpha_{n}}=\frac{\partial N_{\mathrm{d}}}{\partial \alpha_{n}}+\frac{\partial N_{\mathrm{f}}}{\partial \alpha_{n}}+\frac{\partial N_{\mathrm{s}}}{\partial \alpha_{n}} $

(42)

基于式(42)可得到摩擦因子m的表达式

$ m = \frac{{\sqrt {\frac{2}{3}} {\sigma _{\rm{s}}}N\left( {{f_1} + {f_2} + {f_3}} \right) + \frac{{k\tan \theta \Delta bN{b_0}}}{{4{b_1}{h_0}}}\sqrt {1 + \frac{{4b_1^2h_0^2}}{{\Delta {b^2}b_0^2}}} }}{{k\left\{ {\begin{array}{*{20}{l}} {4{v_1}\left( {\sec {\alpha _n} - \cos 2{\alpha _n}} \right) + {v_{\rm{R}}}\left( {\cos 2{\alpha _n} - 1} \right) + }\\ {\frac{{2\Delta v{R^2}}}{{{l^2}}}\left[ {\cos {\alpha _n}\left( {1 + \cos {\alpha _n} - \frac{1}{{\cos {\alpha _n} + 1}}} \right) - \tan {\alpha _n}} \right]} \end{array}} \right\}}} $

(43)

因此，轧制力矩、轧制力以及应力状态系数的解析模型可以根据以下关系进行计算:

$ M=\frac{R}{2 v_{\mathrm{R}}} \varPhi_{\min } ; F=\frac{M}{\chi \sqrt{2 R \Delta h}} ; n_{\sigma}=\frac{p}{2 k}=\frac{F}{4 b l k} $

(44)

3 实验验证

以下用国内某厂实测轧制数据来验证轧制力模型的预测精度。其中，轧件的尺寸为250 mm×2 100 mm×3 500 mm(高×长×宽)，轧辊直径为1 200 mm。在本文中，每道次的实际力臂系数χ分别为0.54、0.54、0.53、0.52和0.51。经过第一道次的整形轧制后进入展宽轧制阶段，此时轧件的初始厚度为235.47 mm。材料为Q345R钢，其变形抗力模型为^[22]:

$ \sigma_{\mathrm{s}}=6、310.7 \bar{\varepsilon}^{-0.407} \bar{\dot{\varepsilon}} \exp \left(-2.62 \times 10^{-3} T-0.669 \bar{\varepsilon}\right) $

(45)

$ T=t+273 $

(46)

式中：ε为等效应变, $\bar {\dot {\varepsilon}} $为等效应变速率, t为轧制温度, T为开尔文温度。

各道次的解析轧制力可由式(44)计算得到，相关的参数以及对比结果见表 1、2。

由表 2易见，轧制力矩计算值与实测值的最大误差均小于6.83%，解析轧制力与实测轧制力的最大误差为7.55%，平均误差为4.77%。各个道次的解析轧制力、轧制力矩与实测值之间的误差均满足工程应用上最大预测误差小于15%的要求，说明本文建立的轧制力模型具有较好预测精度。

将本文所建立的轧制力预测结果与Sims模型^[12]、Tselikov模型^[20]以及实测数据进行对比，结果如图 3所示。由图可见，Sims模型的预测结果偏大，而Tselikov模型的预测结果偏小，两个经典模型的最大预测误差均超过15%。而本文所建立的模型的预测结果与实测轧制力吻合较好，最大误差为6.83%。

4 分析与讨论

图 4为不同压下率下各轧制功率的变化关系图。由图可知，在各压下率下，内部变形功率总是轧制过程中消耗能量的最主要部分，且它在总轧制功率中的占比随着压下率的增加而增加，摩擦功率和剪切功率几乎在同一个数量级且摩擦功率数值最小。摩擦功率和剪切功率虽然都与压下率成正相关，但其数值均小于内部变形功率。当摩擦因子为0.6时，剪切功率与摩擦功率在压下率为0.35时相等，而当摩擦因子为0.8时，二者在压下率为0.3时相等。此外，随着摩擦因子的增加，摩擦功率的增长速率大于剪切功率的增长速率。因此，随着摩擦因子的增加，摩擦功率等于剪切功率时的所需的压下率会下降。

图 5为不同压下率下变形区长度l与中性点位置x_n/l的关系。随着压下率的增大，变形区的总长度虽然变长了，但中性点在变形区内的相对位置并未向出口移动，反而是向变形区中点(图中的中点线)发生了偏移。

图 6为摩擦因子m、压下率ε与轧制力F之间的关系。由图可见，轧制力与摩擦因子以及压下率都成线性正相关。但在压下率较小时，摩擦因子对轧制力变化的影响较小，而在压下率较大时，摩擦因子对轧制力的影响更明显。

图 7为压下率ε、摩擦因子m与中性点位置x_n/l的关系。图中, 摩擦因子对中性点的影响较为显著, 而压下率的影响不显著。可以看出, 随着摩擦因子的增加, 中性点向入口移动。

图 8为形状因子l/2h_m、径厚比R/2h₀与应力状态系数n_σ的关系。图中，应力状态系数与轧件径厚比成负相关关系，而与形状因子成正相关关系。当径厚比较大且形状因子较小时，所需的轧制力较小。

5 结论

1) 本文根据轧件变形区的金属流动特点提出了一个余弦速度场。经证明，该速度场能够严格满足体积不变条件、出入口速度边界条件以及几何方程，这对拟建立的轧制力模型具有一定的精度保障。

2) 针对Mises比塑性功率积分困难的问题，提出了一个新的解决方法，即分矢量内积加和法。经此方法得到了轧制的内部变形功率。同时，也基于提出的速度场得到了摩擦功率和剪切功率的表达式，并建立了轧制力的解析模型。本文的理论轧制力矩、轧制力与各道次实测值的误差均在6.83%和7.55%以内，具有较高精度，表明本文提出的速度场与解析方法是可行的。

3) 参数规律分析表明：在各个压下率下，内部变形功率的消耗都远大于其他两种功率的消耗；随着摩擦因子的增加，摩擦功率的增长速率大于剪切功率的增长速率；轧制中性点的位置并不在变形区中点，而是在中点偏向出口方向的位置上，但随着压下率的增大，中性点会向变形区中点位置移动。