[发明专利]用多极边界元法计算辊式矫直机矫直力的方法

权利要求书

1.一种用多极边界元法计算辊式矫直机矫直力的方法，其特征在于实施步骤为：

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

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

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

σx=23σs(1+lnρlaρmax)]]>

σy=23σslnρlaρmax]]>

σz=σx+σy2=13σs(1+2lnρlaρmax)]]>

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

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

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

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

σ′x=23σs(1+lnρminρya)-P‾]]>

σ′y=23σslnρminρya-P‾]]>

σ′z=σ′x+σ′y2=13σs(1+2lnρminρya)-12P‾]]>

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

——外部受力，

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

中性层向下偏移量为：e=(m-1)ξh2+(1-m2)ξ2h4---(1)]]>

2h——矫直件厚度，

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

ξ——弹区比，

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

F1=4M2b=2σsh2[-mξh-eξh+e(m2ξ2-1)+23(1+m3)ξ2+(1-ξ2)]---(2)]]>

M——矫直弯矩；M=bσsh2[-mξh-eξh+e(m2ξ2-1)+23(1+m3)ξ2+(1-ξ2)];]]>

b——矫直件断面宽度；

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

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

①加载过程：

弹性加载区：σy=2E13rc(y+2ys′)y+23σs′y---(3)]]>

弹塑性加载区：σy=23σs′(y+h)+23A2rc[(y+2ys′)y+h(2ys′-h)]+P‾---(4)]]>

P‾=23A2rc(h-ys′)2-23σs′h-2E13rcys′2---(5)]]>

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

E₁——弹性模量，

r_c——原始曲率半径，

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

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

A₂——加工硬化系数，

——矫直力，

2h——矫直件厚度，

②卸载过程：

弹性加载区：σ*y=2E13rc(y+2ys′)y+23σ*dy---(6)]]>

弹塑性加载区：σ*y=2E13rc[(y+2ys′)y-h(h-2ys′)]+23σ*d(y+h)+P‾---(7)]]>

P‾=2E13rch(h-2ys′)-23σ*dh---(8)]]>σ*_y——卸载时y向应力，

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

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

①初始矫直力：F‾0=P‾0=23A2r0c(h-y0s)2-23σs′h-2E13r0cy0s2---(9)]]>

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

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

②弹性区矫直力：

F‾2i=2A23((h+yis)2ric-(h+y(i-1)s)2r(i-1)c)-2E13(yis2ric+(h+y(i-1)s)2r(i-1)c)-23f1(ϵs)E1yr(i-1)ch---(10)]]>

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

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

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

③弹塑性区矫直力：

F‾2i=2A23((h+yis)2ric-(h+y(i-1)s)2r(i-1)c)-2E13(yis2ric+(h+y(i-1)s)2r(i-1)c)-23[f2(ϵs)σs′+A1+A2yr(i-1)c]h---(11)]]>

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

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

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

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

-D▿4w+F=0]]>

——板的弯曲刚度，

E——材料的弹性模量，

ν——泊松比，

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

板的边界条件为：

在Γ₁上：w-w‾=0βn-β‾n=0βs-β‾s=0]]>在Γ₂上：q-q‾=0mn-m‾n=0ms-m‾s=0]]>

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

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

q——剪力，

mn=n1m1+n2m2=m11n12+2m12n1n2+m22n22]]>

ms=-n2m1+n1m2=-(m11-m22)n1n2+m12(n12-n22)]]>

q=n₁q₁+n₂q₂，βs=-∂w∂s,βn=-∂w∂n,]]>

边界积分方程可写为：

F3∫Ωw*dΩ=Ciwi+∫Γ(mn*βn+t*w)dΓ+Σk=1Ltkc*wk-∫Γ(mnβn*+tw*)dΓ-Σk=1Ltkcwk*---(12)]]>

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

q1*=∂m11*∂x1+∂m12*∂x2,q2*=∂m12*∂x1+∂m22*∂x2,t=q+∂ms∂s,]]>

Σk=1Ltkc*wk=∫Γ(∂ms*∂sw+ms*∂w∂s)dΓ,Σk=1Ltkcwk*=∫Γ(∂ms∂sw*+ms∂w*∂s)dΓ.]]>

计算流程如下:

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

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

③定义摩擦接触约束；

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

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

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

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

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

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

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

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

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

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

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