CN101201871B - 一种预测热轧过程轧制力的刚塑性有限元方法 - Google Patents

Abstract

一种热轧过程预测轧制力的刚塑性有限元方法，根据刚塑性材料变分原理，采用二维刚塑性有限元法建立刚塑性材料能量泛函，求解满足能量泛函得到极小值时的速度场，利用初等方法设定初始速度场，利用黄金分割法进行一维搜索获得修正Newton法的阻尼因子，利用一维变带宽存储法求解线性方程组，根据得到的速度场求解应力场，进而根据轧制条件求解轧制力。这种求解方法未知数数量较少，计算速度较快，与实测结果比较表明计算结果具有较高的计算精度。

Description

一种预测热轧过程轧制力的刚塑性有限元方法

技术领域

本发明属于轧制技术领域，特别涉及一种预测热轧过程轧制力的刚塑性有限元方法。

背景技术

目前，轧制过程预测轧制力的方法多采用数学模型。由于轧制过程中轧件变形参数、力能参数、材料特性等固有的多重非线性关系，长期以来都没有形成精确的数学关系式加以描述和分析。有限元方法能够分析轧制过程详尽的场变量以及轧制力和前滑等参数积分量。但是通用的有限元软件由于其固有的求解方程系统复杂，专业应用性不强以及采用弹塑性或其它不同的材料模型等使得计算效率低，轧制力计算精度不高。

发明内容

为了克服上述方法预测精度低和计算效率低的缺点，本发明提供一种预测热轧过程轧制力的刚塑性有限元方法，其目的是提高轧制力预测精度，提高计算效率。

本发明的技术解决方案包括如下步骤：

(1)通过刚塑性材料变分原理，采用二维刚塑性可压缩有限元法求解板材轧制问题，建立可压缩法的能量泛函为：

$\mathrm{\Φ}=\frac{1}{m+1}{\∫\∫}_A\stackrel{\‾}{\mathrm{\σ}}\stackrel{\·}{\stackrel{\‾}{\mathrm{\ϵ}}}\mathrm{dA}+{\∫}_{\mathrm{Lf}}{\mathrm{\τ}}_f\mathrm{\Δ}{v}_f\mathrm{dl}\±{\∫}_{\mathrm{Lv}}{T}_1{v}_s\mathrm{dl}={\mathrm{\φ}}^P+{\mathrm{\φ}}^f+{\mathrm{\φ}}^t---\left(1\right)$

其中：

φ^p-轧件内部的塑性变形功率

-等效应力

-等效变形速度

m-速度敏感指数

φ^f-轧件与轧辊之间的摩擦功率

ΔV_f-轧件与轧辊之间的相对滑动速度

τ_f-摩擦剪应力

φ^t-外张力功率

T-张应力，‘+’为前张力，‘-’为后张力

v_s-相应表面处的位移速度。

对刚塑性可压缩有限元法来说，刚塑性可压缩材料等效变形速度表达式为：

$\stackrel{\·}{\stackrel{\‾}{\mathrm{\ϵ}}}=\frac{4}{9}({\stackrel{\·}{\mathrm{\ϵ}}}_x^2+{\stackrel{\·}{\mathrm{\ϵ}}}_y^2-{\stackrel{\·}{\mathrm{\ϵ}}}_x{\stackrel{\·}{\mathrm{\ϵ}}}_y)+\frac{1}{3}{\stackrel{\·}{\mathrm{\γ}}}_{\mathrm{xy}}+\frac{1}{g}{\stackrel{\·}{\mathrm{\ϵ}}}_v^2---\left(2\right)$

(2)对能量泛函求极值，获得速度场修正量求解的线性方程组和速度场求解的迭代方程

把未知数序列(v₁…v_n)记为矢量v，以v_k表示第k迭代步中得出的近似解，将泛函Φ＝f(v)在v＝v_k邻域内以泰勒级数展开取前三项：

$\mathrm{\Φ}=f\left(v\right)\≈f\left(v_k\right)+\▿f\left(v_k\right)(v-v_k)+\frac{1}{2}{\▿}^{2}f\left(v_k\right)\·{(v-v_k)}^2---\left(3\right)$

其中：v-v_k＝Δv_k            (4)

那么，泛函Φ即为Δv_k的二次函数：

$\mathrm{\Φ}\≈f\left(v_k\right)+\▿f\left(v_k\right)\·\mathrm{\Δ}{v}_k+\frac{1}{2}{\▿}^{2}f\left(v_k\right)\·\mathrm{\Δ}{v}_k---\left(5\right)$

将式(3)泛函对Δv_k求一阶偏导并置零，它的极值可以通过求解线性方程组(6)获得：

${\▿}^{2}f\left(v_k\right)\·\mathrm{\Δ}{v}_k=-\▿f\left(v_k\right)---\left(6\right)$

由式(6)解出的Δv_k为速度修正量，可使Φ接近其极值。由式(4)取k+1迭代步中

v_k+1＝v_k+Δv_k            (7)

得到新的速度场v_k+1更加接近真实解，求解过程反复迭代，直到满足收敛条件，Δv→0，此时v即为最终解

或者在由式(4)取K+1迭代步中，采用修正牛顿法，修正牛顿法是在速度修正量前加一个阻尼因子α，使迭代过程变为：

v_k+1＝v_k+αΔv_k                        (8)

(3)利用刚塑性有限元法，利用能量泛函和牛顿迭代法求解板材轧制力按以下步骤进行：

①采集轧制条件原始数据、单元划分数据和收敛条件

轧制条件原始数据包括：

轧辊直径(mm)，轧制速度(mm/s)，入口厚度(mm)，出口厚度(mm)，板坯宽度(mm)，轧制温度(℃)，前张力和后张力(MPa)，摩擦因子，摩擦系数，可压缩因子。

单元划分数据包括：轧制方向单元数目，厚度方向单元数目

收敛条件包括：速度收敛条件和能量泛函收敛条件

②进行单元节点划分、编号和调查，计算节点坐标，设定速度边界条件。

根据有限元法计算需要和单元划分数据对轧件进行四边形线性等参单元均匀划分、编号和调查，并计算节点坐标，如图1所示。图中α为板材轧制的咬入角，β为接触弧BD上任一点的接触角，H0，H1分别为轧制入口和轧制出口半厚度，R为轧辊半径，i为单元，j为节点，X方向为轧制方向，Y方向为厚度方向。

单元和节点编号规则为从A点沿厚度方向和轧制方向依次增加，如图1所示。

节点坐标计算方法为：以A点为O点，对每一个节点进行坐标计算。轧制方向上每个单元的宽度相同，厚度方向上每相同列的单元高度相同。

为了简化计算，制定已知速度边界条件为：线段AB上每个节点的轧向速度v_x相同，厚向速度v_y为零；线段CD上每个节点的轧向速度v_x相同，厚向速度v_y为零；线段AC上的厚向速度v_y为零；接触弧BD上节点的轧向节点速度v_x与厚向速度v_y满足：v_y＝-v_xtanβ。

③利用初等方法进行初速度场设定，初等方法主要是结合板带轧制条件(摩擦条件和咬入角)和工程法设定近似初始速度场。

初等方法设定初速度场主要假设如下：

(1)同一垂直截面上节点速度在轧制方向上的速度分量v_x相同(平断面假设)

(2)沿厚度方向上速度分量v_y呈线性分布

(3)沿轧制方向上任意一垂直截面的金属秒流量体积相等。

首先确定一垂直面的秒流量，由于秒流量与中性面有关，以中性面秒流量为基准设定初始速度场。中性面的中性角依据工程法确定：

$\mathrm{\γ}=\frac{\mathrm{\α}}{2}(1-\frac{\mathrm{\α}}{2t_f})---\left(9\right)$

其中：r-中性角；t_f-为摩擦系数。

中性面处的秒体积为：

V_N＝h_N·v_R·cosγ            (10)

其中：V_N表示中性面处秒体积，v_R轧制速度，h_N表示中性面处轧件半厚度。由假设条件可以获得第i列所有节点的轧向速度：

h_i·v_xi·cosβ＝h_N·v_R·cosγ

$\⇒v_{\mathrm{xi}}=v_R\·\frac{h_N}{h_i}\·\frac{\cos \mathrm{\γ}}{\cos \mathrm{\β}}---\left(11\right)$

然后，由假设条件和速度边界条件，可以获得内节点j处的轧向速度v_yj

$v_{\mathrm{yj}}=2v_{\mathrm{xi}}\·(-\tan \mathrm{\β})\·\frac{y_j}{h_i}---\left(12\right)$

④计算能量泛函Φ，求解能量泛函一阶和二阶偏导数。

利用有限元法中局部坐标和整体坐标映射基本原理，求解等参四边形单元的形函数，B矩阵和雅克比矩阵J及其行列式|J|，然后利用高斯积分法，对能量泛函式(1)进行计算。经过高斯积分后塑性功率和摩擦功率的单元能耗率泛函如下：

${\mathrm{\φ}}_e=\underset{j=1}{\overset{4}{\mathrm{\Σ}}}\stackrel{\‾}{\mathrm{\σ}}\stackrel{\·}{\stackrel{\‾}{\mathrm{\ϵ}}}\left|J\right|+\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{m}_f{\mathrm{\τ}}_f\·\mathrm{\Δ}{v}_f\·\frac{x_2-x_1}{2}---\left(13\right)$

为了充分考虑温度和材料化学成分对轧制过程影响，从而提高计算精度，式(13)中，变等效应力

采用志田模型，τ_f-摩擦剪应力，

m为摩擦因子。

张力功率泛函为：

φ_t＝t₀v_x1-t₁v_xnp                (14)

其中：e表示每个单元；j表示节点数；m_f-摩擦系数；x₂，x₁为接触面上相邻节点x坐标；v_x1，v_xnp分别表示入口端和出口端轧向速度，N为厚向单元数；t₀，t₁分别表示前张力和后张力。

对方程式(13)和(14)分别求一阶和二阶偏导数，从而获得与速度场相关的线性方程组(6)的系数矩阵和常数矩阵。

⑤对形成的线性方程组进行速度修正量求解，在速度修正量求解中，采用一维变带宽存储法求解大型线性方程组，从而节约存储空间，提高计算效率。一维变带宽存储方法可以说明如下：对一般的线性方程组[A]X＝[B]，如果系数矩阵为[A]对称稀疏矩阵(零元素较多，)，见式(15)。可以从第一列非零元素开始进行存储，直到对角线元素为止

存储为一维矩阵为：

[A]＝[a₁₁ a₁₂ a₂₂ a₂₃ a₃₃ a₁₄ a₂₄ 0 a₄₄ a₃₅ 0 a₅₅… a_n-1，n  a_nn]

其中：a-为矩阵元素，下表第一项表示行数，第二项表示列数。

⑥采用黄金分割法进行阻尼因子一维线性搜索

为了提高计算精度和计算速度，在利用修正牛顿法求解时，阻尼因子通过黄金分割法进行一维线性搜索获得，黄金分割法被证明是一种搜索α因子的高效率算法。

⑦根据修正牛顿迭代法，利用式(8)求解速度场

⑧由速度场求出变形速度

获得速度修正量和阻尼因子后，可以由式(8)从前一迭代步得到下一迭代步的速度场，即获得每个节点的变形速度。由变形速度和节点速度关系式计算变形速度。

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\ϵ}}}_x\\ {\stackrel{\·}{\mathrm{\ϵ}}}_y\\ {\stackrel{\·}{\mathrm{\γ}}}_{\mathrm{xy}}\end{array}\right\}=\left[B\right]\left\{v\right\}---\left(16\right)$

其中：

方向变形速度；

方向变形速度；，-剪切变形速度；[B]-为B矩阵；{v}节点速度矢量。

⑨根据收敛条件进行收敛判定

根据所得变形速度场，由式(1)求出第k迭代步的能量泛函，然后和k-1步的能量泛函组成能量泛函收敛条件：

$\frac{{\mathrm{\&Delta;\&phi;}}_k}{{\mathrm{\&phi;}}_k}=\frac{|{\mathrm{\&phi;}}_k-{\mathrm{\&phi;}}_{k-1}|}{{\mathrm{\&phi;}}_k}<{\mathrm{\&epsiv;}}_{\mathrm{\&phi;}}---\left(17\right)$

由第k迭代步的速度场修正量和速度场，计算速度收敛条件：

$\frac{\left|\right|\mathrm{\&Delta;}{v}_k\left|\right|}{\left|\right|v_k\left|\right|}=\frac{{\left(\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{v}_j^2\right)}^{\frac{1}{2}}}{{\left(\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_j^2\right)}^{\frac{1}{2}}}<{\mathrm{\&epsiv;}}_v---\left(18\right)$

式中：

k为迭代步数

Δφ_k为能量泛函变化率

φ_k为能量泛函

ε_φ为能量泛函收敛条件

j为节点编号

n为总节点数

ε_v为速度收敛条件

Δv_j为节点速度变化率

v_j为节点速度

⑨计算应力场

当收敛条件满足后，认为速度场为真实速度场，然后由速度场通过式(16)求解变形速度，通过变形速度由式(19)求解应力场：

${\mathrm{\&sigma;}}_x=\frac{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}}[\frac{2}{3}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_x+(\frac{1}{g}-\frac{2}{9}){\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_v]---\left(19a\right)$

${\mathrm{\&sigma;}}_y=\frac{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}}[\frac{2}{3}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_y+(\frac{1}{g}-\frac{2}{9}){\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_v]---\left(19b\right)$

${\mathrm{\&tau;}}_{\mathrm{xy}}=\frac{1}{3}\&CenterDot;\frac{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}}\&CenterDot;{\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}_{\mathrm{xy}}---\left(19c\right)$

式中：

-等效应力；

σ_x-x方向应力；

σ_y-y方向应力；

τ_xy-剪应力。

⑩计算轧制力

根据求解所得的应力场，由压下方向上的正应力沿接触表面积分得到轧制力，其均值为平均单位压力，轧制力和轧制压力表达式为：

$F=b\underset{0}{\overset{l}{\&Integral;}}{\mathrm{\&sigma;}}_y\mathrm{dx}---\left(20\right)$

其中：F为轧制力(KN)，l为轧制过程接触弧长(m)，b为板带宽度(m)。

本发明的计算流程如图2所示。

本发明的最大效果是：能够较为准确预测热轧过程轧制力，为轧制过程其它参数设定和优化提供必要的信息；本方法应用性强，计算时间短，提高了计算效率。本发明适用于热连轧过程的粗轧机和精轧机。

附图说明

图1本发明的有限元模型图；

图2本发明的计算流程图；

图中：α为板材轧制的咬入角，β为接触弧BD上任一点的接触角，H0，为轧制入口半厚度，H1为轧制出口半厚度，R为轧辊半径，i为单元，j为节点，X方向为轧制方向，Y方向为厚度方向。

具体实施方式

采集轧制条件原始数据，如表1所示，单元节点划分为18×6，摩擦系数为0.5，可压缩因子为0.01。

表1

按本发明方法步骤，将采集的轧制条件原始数据输入计算机，依图2所示，顺序计算。

表2为计算结果，从表2可以看出，轧制力计算结果与实测值吻合良好，计算误差小于10％，轧制一道次计算时间少于1000ms，计算精度和计算时间能够满足现场需求。

表2

编号	实测轧制力(KN)	计算轧制力(KN)	计算误差(％)	计算时间(ms)
					1	17020	17805	4.6	328
2	26050	25179	-3.3	525
					3	26450	26035	-1.6	515
4	24410	25170	3.1	203
					5	20200	21761	7.7	375
6	23820	22845	-4.1	265

编号	实测轧制力(KN)	计算轧制力(KN)	计算误差(％)	计算时间(ms)
					7	25100	24290	-3.2	196
8	27210	28237	3.8	359
					9	19520	20837	6.7	375
10	15850	16780	5.9	359
					11	21440	23341	8.9	328
12	16680	17782	6.6	265

本发明方法针对板带热轧过程可以获得轧制过程的轧制力，计算值与跟踪测量值吻合良好。如表2所示，提高了热轧过程轧制力的预测精度和效率，能够更好的优化热轧过程参数。

Claims (1)

1.一种预测热轧过程轧制力的刚塑性有限元方法，其特征是利用刚塑性可压缩有限元法，利用能量泛函和牛顿迭代方法求解板材轧制力，包括以下步骤：

①采集轧制条件原始数据、单元划分数据和收敛条件；其中轧制条件原始数据包括：轧辊直径、轧制速度、入口厚度、出口厚度、板坯宽度、轧制温度，前张力和后张力，摩擦因子，摩擦系数和可压缩因子；单元划分数据包括：轧制方向单元数目和厚度方向单元数目，收敛条件包括：速度收敛条件和能量泛函收敛条件；

②进行单元节点划分、编号和调查，计算节点坐标，设定速度边界条件；

③利用初等方法进行初速度场设定；初等方法设定初速度场主要假设如下：

(1)同一垂直截面上节点速度在轧制方向上的速度分量v_x相同；

(2)沿厚度方向上速度分量v_y呈线性分布；

(3)沿轧制方向上任意一垂直截面的金属秒流量体积相等；

④通过刚塑性材料变分原理，采用二维刚塑性可压缩有限元法，建立刚塑性材料能量泛函，计算能量泛函Φ，求解能量泛函一阶和二阶偏导数；

能量泛函为：

$\mathrm{\&Phi;}=\frac{1}{m+1}{\&Integral;\&Integral;}_A\stackrel{\&OverBar;}{\mathrm{\&sigma;}}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}\mathrm{dA}+{\&Integral;}_{\mathrm{Lf}}{\mathrm{\&tau;}}_f\mathrm{\&Delta;}{v}_f\mathrm{dl}\&PlusMinus;{\&Integral;}_{\mathrm{Lv}}{T}_1{v}_s\mathrm{dl}={\mathrm{\&phi;}}^P+{\mathrm{\&phi;}}^f+{\mathrm{\&phi;}}^t---\left(1\right)$

式中：

φ^p-轧件内部的塑性变形功率；

-等效应力；

-等效变形速度；m-速度敏感指数；

φ^f-轧件与轧辊之间的摩擦功率；Δv_f-轧件与轧辊之间的相对滑动速度；

τ_f-摩擦剪应力；φ^t-外张力功率；T₁-张应力，v_s-相应表面处的位移速度；

⑤采用一维变带宽存储法求解大型线性方程组；

⑥采用黄金分割法进行阻尼因子一维线性搜索；

⑦对能量泛函求极值，利用牛顿法反复迭代获得真实速度场，由速度场求出变形速度：

把未知数序列(v₁…v_n)记为矢量v，以v_k表示第k迭代步中得出的近似解，将泛函Φ＝f(v)在v＝v_k邻域内以泰勒级数展开取前三项：

$\mathrm{\&Phi;}=f\left(v\right)\&ap;f\left(v_k\right)+\&dtri;f\left(v_k\right)(v-v_k)+\frac{1}{2}{\&dtri;}^{2}f\left(v_k\right)\&CenterDot;{(v-v_k)}^2---\left(3\right)$

其中：

v-v_k＝Δv_k    (4)

那么，泛函Φ即为Δv_k的二次函数：

$\mathrm{\&Phi;}\&ap;f\left(v_k\right)+\&dtri;f\left(v_k\right)\&CenterDot;\mathrm{\&Delta;}{v}_k+\frac{1}{2}{\&dtri;}^{2}f\left(v_k\right)\&CenterDot;\mathrm{\&Delta;}{v}_k---\left(5\right)$

将式(3)泛函对Δv_k求一阶偏导并置零，它的极值可以通过求解线性方程组(6)获得：

${\&dtri;}^{2}f\left(v_k\right)\&CenterDot;\mathrm{\&Delta;}{v}_k=-\&dtri;f\left(v_k\right)---\left(6\right)$

由式(6)解出的Δv_k为速度修正量，可使Φ接近其极值，由式(4)取k+1迭代步中

v_k+1＝v_k+Δv_k    (7)

求解过程反复迭代，直到满足收敛条件，Δv→0，此时v即为最终解；

或者在由式(4)取K+1迭代步中，采用修正牛顿法，修正牛顿法是在速度修正量前加一个阻尼因子α，使迭代过程变为：

v_k+1＝v_k+αΔv_k

⑧进行收敛判定，计算速度收敛条件；

⑨计算应力场；

通过变形速度依下式求解应力场：

${\mathrm{\&sigma;}}_x=\frac{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}}[\frac{2}{3}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_x+(\frac{1}{g}-\frac{2}{9}){\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_v]$

${\mathrm{\&sigma;}}_y=\frac{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}}[\frac{2}{3}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_y+(\frac{1}{g}-\frac{2}{9}){\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_v]$

${\mathrm{\&tau;}}_{\mathrm{xy}}=\frac{1}{3}\&CenterDot;\frac{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}}\&CenterDot;{\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}_{\mathrm{xy}}$

式中：x方向变形速度，

y方向变形速度，剪切变形速度，

为速度收敛条件，

-等效应力，

-等效变形速度；σ_x-x方向应力；σ_y-y方向应力；τ_xy-剪应力；

⑩计算轧制力

根据求解所得的应力场，计算轧制力，依下式进行：

$F=b\underset{0}{\overset{l}{\&Integral;}}{\mathrm{\&sigma;}}_y\mathrm{dx}$

式中：F为轧制力，l为轧制过程接触弧长，b为板带宽度。

Images