用多极边界元法计算辊式矫直机矫直力的方法

摘要

一种用多极边界元法计算辊式矫直机矫直力的方法，属于辊式矫直的技术领域。特征是首先计算考虑中性层偏移的矫直力；其次计算考虑包辛格效应的矫直力；再次利用自主开发的三维弹塑性多物体接触多极边界元法程序包计算矫直力；最后综合考虑三方面因素，利用各自在矫直力中所占权重最终求得辊式矫直过程中矫直力大小。该发明的优点是得到更为精确的矫直力计算公式，为矫直机设计提供理论支持，矫直力计算更为精确，有效提高矫直效果，进而提高产品质量。

说明书

技术领域

本发明属于辊式矫直的技术领域，具体涉及一种利用多极边界元法计算辊式矫直机矫直力的方法。 

背景技术

在辊式矫直中矫直力的精确度决定了矫直效果的好坏，板材的产品质量，最终决定了企业在板材市场的竞争能力。目前矫直力的数值计算方法主要是有限元法，但由于影响矫直力的因素有很多，并且非线性程度高，导致有限元计算结果并不理想。进而影响了矫直力模型的建立，矫直效果也受到一定影响，不利于提高产品质量和提升企业的竞争力。 

发明内容

本发明目的是提供一种用多极边界元法计算辊式矫直机矫直力的方法，可以有效提高矫直力的计算精度，保证矫直效果。 

本发明是这样实现的： 

<1>考虑中性层偏移的矫直力计算公式推导 

对矫直件进行三维弹塑性接触摩擦分析，利用拉伸区和压缩区微元体的力学平衡关系确定中性层偏移计算公式。以几何中心层为界分为拉伸区和压缩区两部分,在拉伸区和压缩区分别取微元体dα×dρ，反弯曲率半径为ρ_w，公式推导结果如下： 

拉伸区微元体三向主应力分别为： 

${\sigma }_x=\frac{2}{\sqrt{3}}{\sigma }_s(1+\ln \frac{{\rho }_{\mathrm{la}}}{{\rho }_{\max }})$

${\sigma }_y=\frac{2}{\sqrt{3}}{\sigma }_s\ln \frac{{\rho }_{\mathrm{la}}}{{\rho }_{\max }}$

${\sigma }_z=\frac{{\sigma }_x+{\sigma }_y}{2}=\frac{1}{\sqrt{3}}{\sigma }_s(1+2\ln \frac{{\rho }_{\mathrm{la}}}{{\rho }_{\max }})$

ρ_la——拉伸区微元体恰好达到屈服极限时的曲率半径， 

ρ_max——压弯后矫直件上表面曲率半径， 

σ_s——矫直件的屈服应力。 

压缩区微元体三向主应力分别为： 

${{\sigma }^{\prime }}_x=\frac{2}{\sqrt{3}}{\sigma }_s(1+\ln \frac{{\rho }_{\min }}{{\rho }_{\mathrm{ya}}})-\bar{P}$

${{\sigma }^{\prime }}_y=\frac{2}{\sqrt{3}}{\sigma }_s\ln \frac{{\rho }_{\min }}{{\rho }_{\mathrm{ya}}}-\bar{P}$

${{\sigma }^{\prime }}_z=\frac{{{\sigma }^{\prime }}_x+{{\sigma }^{\prime }}_y}{2}=\frac{1}{\sqrt{3}}{\sigma }_s(1+2\ln \frac{{\rho }_{\min }}{{\rho }_{\mathrm{ya}}})-\frac{1}{2}\bar{P}$

ρ_min——压弯后矫直件下表面曲率半径， 

——外部受力， 

ρ_ya——压缩区微元体恰好达到屈服极限时的曲率半径。 

中性层向下偏移量为：$e=\frac{(m-1)\mathrm{\xi h}}{2}+\frac{(1-m^2){\xi }^{2}h}{4}---\left(1\right)$

2h——矫直件厚度， 

m——偏移系数，y_s'为矫直件压缩区弹性变形厚度，y_s为矫直件拉伸区弹性变形厚度 

ξ——弹区比，

在考虑中性层偏移的情况下，矫直力计算公式为： 

$F_1=\frac{4M}{2b}=2{\sigma }_s{h}^2[-\frac{\mathrm{m\xi h}-e}{\mathrm{\xi h}+e}(m^2{\xi }^2-1)+\frac{2}{3}(1+m^3){\xi }^2+(1-{\xi }^2)]---\left(2\right)$

M——矫直弯矩；$M={\mathrm{b\sigma }}_s{h}^2[-\frac{\mathrm{m\xi h}-e}{\mathrm{\xi h}+e}(m^2{\xi }^2-1)+\frac{2}{3}(1+m^3){\xi }^2+(1-{\xi }^2)];$

b——矫直件断面宽度。 

<2>考虑包辛格效应的矫直力计算公式推导 

在直角坐标系下加载和卸载过程进行分析，分别给出弹性区和弹塑性区应力和矫直力的表达式，具体公式推导如下： 

①加载过程： 

弹性加载区：${\sigma }_y=\frac{{2E}_1}{{3r}_c}(y+{2y}_{s^{\prime }})y+\frac{2}{\sqrt{3}}{\sigma }_{s^{\prime }}y---\left(3\right)$

弹塑性加载区：${\sigma }_y=\frac{2}{\sqrt{3}}{\sigma }_{s^{\prime }}(y+h)+\frac{2}{3}\frac{A_2}{r_c}[(y+{2y}_{s^{\prime }})y+h({2y}_{s^{\prime }}-h)]+\bar{P}---\left(4\right)$

$\bar{P}=\frac{2}{3}\frac{A_2}{r_c}{(h-y_{s^{\prime }})}^2-\frac{2}{\sqrt{3}}{\sigma }_{s^{\prime }}h-\frac{{2E}_1}{{3r}_c}{y}_{s^{\prime }}^2---\left(5\right)$

σ_y——加载时y向应力， 

E₁——弹性模量， 

r_c——原始曲率半径， 

y_s'——压缩区弹塑性界面到中心层的厚度， 

σ_s'——加载屈服应力， 

A₂——加工硬化系数， 

——矫直力， 

2h——矫直件厚度。 

②卸载过程： 

弹性加载区：${\sigma *}_y=\frac{{2E}_1}{{3r}_c}(y+{2y}_{s^{\prime }})y+\frac{2}{\sqrt{3}}{\sigma *}_{d}y---\left(6\right)$

弹塑性加载区：${\sigma *}_y=\frac{{2E}_1}{{3r}_c}[(y+{2y}_{s^{\prime }})y-h(h-{2y}_{s^{\prime }})]+\frac{2}{\sqrt{3}}{\sigma *}_d(y+h)+\bar{P}---\left(7\right)$

$\bar{P}=\frac{{2E}_1}{{3r}_c}h(h-{2y}_{s^{\prime }})-\frac{2}{\sqrt{3}}{\sigma *}_{d}h---\left(8\right)$

σ*_y——卸载时y向应力， 

σ*_d——卸载屈服应力。 

考虑包辛格效应的矫直力模型为： 

①初始矫直力：${\bar{F}}_0={\bar{P}}_0=\frac{2}{3}\frac{A_2}{r_{0c}}{(h-y_{0s})}^2-\frac{2}{\sqrt{3}}{\sigma }_{s^{\prime }}h-\frac{{2E}_1}{{3r}_{0c}}{y}_{0s}^2---\left(9\right)$

r_0c——矫直件的初始曲率半径， 

y_0s——初始状态弹性厚度； 

②弹性区矫直力： 

$\left(\begin{array}{c}{\bar{F}}_2^i=\frac{{2A}_2}{3}(\frac{{(h+y_{\mathrm{is}})}^2}{r_{\mathrm{ic}}}-\frac{{(h+y_{(i-1)s})}^2}{r_{(i-1)c}})-\frac{{2E}_1}{3}(\frac{y_{\mathrm{is}}^2}{r_{\mathrm{ic}}}+\frac{{(h+y_{(i-1)s})}^2}{r_{(i-1)c}})\\ -\frac{2}{\sqrt{3}}{f}_1\left({ϵ}_s\right)E_1\frac{y}{r_{(i-1)c}}h\end{array}\right)---\left(10\right)$

y_is——经过第i辊时钢板压缩区的弹性厚度，当压缩区在钢板中心层下方时取负，压缩区在钢板中心层上方时取正。 

r_ic——经过第i辊时矫直件的曲率半径， 

f₁(ε_s)——杨氏模量减少量； 

③弹塑性区矫直力： 

$\left(\begin{array}{c}{\bar{F}}_2^i=\frac{{2A}_2}{3}(\frac{{(h+y_{\mathrm{is}})}^2}{r_{\mathrm{ic}}}-\frac{{(h+y_{(i-1)s})}^2}{r_{(i-1)c}})-\frac{{2E}_1}{3}(\frac{y_{\mathrm{is}}^2}{r_{\mathrm{ic}}}+\frac{{(h+y_{(i-1)s})}^2}{r_{(i-1)c}})\\ -\frac{2}{\sqrt{3}}[f_2\left({ϵ}_s\right){\sigma }_{s^{\prime }}+A_1+A_2\frac{y}{r_{(i-1)c}}]h\end{array}\right)---\left(11\right)$

f₂(ε_s)——反向屈服应力减少量。 

<3>利用自主开发的三维弹塑性多物体接触多极边界元法程序包计算矫直力 

利用fortran语言开发三维弹塑性多物体接触多极边界元法程序包，用于模拟辊式矫直过程。在矫直过程中，将轧件变形视为板弯曲问题的反问题。先建立板弯曲问题的多极边界元法的边界积分方程，利用板弯曲问题的反问题求 解矫直力。 

设板的边界为Γ，内部区域为Ω。用w表示板中面的挠度，则板弯曲问题的控制方程可用w做变量表示为： 

${-D▿}^{4}w+F=0$

——板的弯曲刚度， 

E——材料的弹性模量， 

ν——泊松比， 

F——垂直作用于板面的压力。 

板的边界条件为： 

在Γ₁上：$\left(\begin{array}{c}w-\bar{w}=0\\ {\beta }_n-{\bar{\beta }}_n=0\\ {\beta }_s-{\bar{\beta }}_s=0\end{array}\right)$在Γ₂上：$\left(\begin{array}{c}q-\bar{q}=0\\ m_n-{\bar{m}}_n=0\\ m_s-{\bar{m}}_s=0\end{array}\right)$

β_n，β_s——板边界上任意点处绕法线和切线的转角， 

m_n，m_s——板截面上的弯矩和扭矩， 

q——剪力。 

$m_n=n_1{m}_1+n_2{m}_2=m_{11}{n}_1^2+{2m}_{12}{n}_1{n}_2+m_{22}{n}_2^2$

$m_s=-n_2{m}_1+n_1{m}_2=-(m_{11}-m_{22})n_1{n}_2+m_{12}(n_1^2-n_2^2)$

q=n₁q₁+n₂q₂，${\beta }_s=-\frac{\partial w}{\partial s},{\beta }_n=-\frac{\partial w}{\partial n}.$

边界积分方程可写为： 

$F_3{\int }_{\Omega }{w}^{*}\mathrm{d\Omega }=C_i{w}_i+{\int }_{\Gamma }(m_n^{*}{\beta }_n+t^{*}w)\mathrm{d\Gamma }+\underset{k=1}{\overset{L}{\Sigma }}{t}_k^{c^{*}}{w}_k-{\int }_{\Gamma }(m_n{\beta }_n^{*}+{\mathrm{tw}}^{*})\mathrm{d\Gamma }-\underset{k=1}{\overset{L}{\Sigma }}{t}_k^c{w}_k^{*}---\left(12\right)$

w^*——单位集中力沿x₃轴方向作用于无穷大薄板的ξ点时，板的几何中心 层上任意一点产生的挠度，

$q_1^{*}=\frac{{\partial m}_{11}^{*}}{{\partial x}_1}+\frac{{\partial m}_{12}^{*}}{{\partial x}_2},q_2^{*}=\frac{{\partial m}_{12}^{*}}{{\partial x}_1}+\frac{{\partial m}_{22}^{*}}{{\partial x}_2},t=q+\frac{{\partial m}_s}{\partial s},$

$\underset{k=1}{\overset{L}{\Sigma }}{t}_k^{c^{*}}{w}_k={\int }_{\Gamma }(\frac{{\partial m}_s^{*}}{\partial s}w+m_s^{*}\frac{\partial w}{\partial s})\mathrm{d\Gamma },\underset{k=1}{\overset{L}{\Sigma }}{t}_k^c{w}_k^{*}={\int }_{\Gamma }(\frac{{\partial m}_s}{\partial s}{w}^{*}+m_s\frac{{\partial w}^{*}}{\partial s})\mathrm{d\Gamma }.$

如图3所示，计算流程如下所述: 

①将矫直辊和矫直件划分单元,将单元和节点信息以文本形式读入，并确定边界条件； 

②输入增量步数据并根据接触条件进行多物体接触检查； 

③定义摩擦接触约束； 

④求解单元最大长度，根据单元长度分别采用指数展开和球谐函数展开求多极矩系数； 

⑤利用GMRES(m)求解方程组（12），并检验是否满足精度，若满足则执行步骤（6），否则返回步骤（3）改变单元长度重新计算； 

⑥计算表面应力，判断摩擦状态是否改变。若摩擦状态改变则返回步骤（4）重新计算，否则执行步骤（7）； 

⑦计算内点应力，判断是否发生穿透。若发生穿透则返回步骤（3），细分增量步，否则执行步骤（8）； 

⑧计算是否为最后增量步，是则结束计算，否则返回步骤（2）继续迭代计算。 

<4>综合考虑多极边界元法计算的数值结果和第一部分推导的关于中性层偏移量和包辛格效应的计算公式，最终得到矫直力计算公式： 

F=α₁F₁+α₂F₂+α₃F₃    （13） 

α₁，α₂，α₃——权重系数，取决于具体工况； 

F₁——由中性层偏移决定的矫直力； 

F₂——由包辛格效应决定的矫直力； 

F₃——由多极边界元法计算得到的矫直力。 

本发明的优点是：建立更为精确的矫直力模型，提高矫直力计算精度和矫直效率，提高产品质量。 

附图说明

图1为矫直件微元体受力分析示意图 

图2为包含中性层偏移的应力分布图 

图3为计算流程图 

图中dα×dρ——微元体尺寸， 

ρ_w——反弯曲率半径， 

ρ_la——拉伸区微元体恰好达到屈服极限时的曲率半径， 

ρ_max——压弯后矫直件上表面曲率半径， 

ρ_min——压弯后矫直件下表面曲率半径， 

ρ_ya——压缩区微元体恰好达到屈服极限时的曲率半径， 

2h——矫直件厚度， 

σ_s——矫直件拉伸区屈服应力， 

σ_s'——矫直件压缩区屈服应力， 

y_s'——矫直件压缩区弹性变形厚度， 

y_s——矫直件拉伸区弹性变形厚度， 

e——中性层偏移量， 

A——矫直件几何中心层， 

A＇——矫直件应力中性层。 

具体实施方式

本实施例是用来说明发明的，而不是对本发明做任何限制。 

在实验室11辊矫直机实验平台上进行实验验证。选用Q235钢种做为矫直件材料，矫直件厚度为40mm，宽度为600mm，弹性模量为2.06E+11，泊松比为0.3，入、出口压弯量为-2.08mm/0mm。理论计算矫直力为1921.67KN，实测矫直力为1969.47KN，测量矫直力与计算值误差为2.49%。 

过程如下： 

（1）根据产品规格，利用公式（1）得到中性层偏移量为0.215mm，根据公式（2）求出包含中性层偏移量的矫直力计算公式F₁为1466.4KN； 

（2）根据公式（3，4）和公式（6，7）求出在加载过程中和卸载过程中材料的包辛格值，并根据公式（9,10，11）分别计算出在中性层偏移和包辛格效应共同作用下矫直力F₂为1854.66KN。 

（3）利用公式（12），按照流程图3，利用多极边界元法计算矫直力F₃为2016.95KN。根据公式（13）可以求得矫直力大小。 

F=α₁F₁+α₂F₂+α₃F₃

=0.12*1466.4+0.18*1854.66+0.7*2016.95 

=1921.67KN。