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

摘要

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

主权项

一种用多极边界元法计算辊式矫直机矫直力的方法，其特征在于实施步骤为：<1>考虑中性层偏移的矫直力计算公式推导对矫直件进行三维弹塑性接触摩擦分析，利用拉伸区和压缩区微元体的力学平衡关系确定中性层偏移计算公式，以几何中心层为界分为拉伸区和压缩区两部分,在拉伸区和压缩区分别取微元体dα×dρ，反弯曲率半径为ρw，公式推导结果如下：拉伸区微元体三向主应力分别为：ρla——拉伸区微元体恰好达到屈服极限时的曲率半径，ρmax——压弯后矫直件上表面曲率半径，σs——矫直件的屈服应力，压缩区微元体三向主应力分别为：ρmin——压弯后矫直件下表面曲率半径，——矫直力，ρya——压缩区微元体恰好达到屈服极限时的曲率半径，中性层向下偏移量为：式中：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>考虑包辛格效应的矫直力计算公式推导在直角坐标系下加载和卸载过程进行分析，分别给出弹性区和弹塑性区应力和矫直力的表达式，具体公式推导如下：①加载过程：弹性加载区：弹塑性加载区：P‾=23A2rc(h-ys′)2-23σs′h-2E13rcys′2---(5)]]>式中：σy——加载时y向应力，E1——弹性模量，rc——原始曲率半径，ys'——压缩区弹塑性界面到中心层的厚度，σs’——加载屈服应力，A2——加工硬化系数；②卸载过程：弹性加载区：弹塑性加载区：P‾=2E13rch(h-2ys′)-23σ*dh---(8)]]>σ*y——卸载时y向应力，σ*d——卸载屈服应力；考虑包辛格效应的矫直力模型为：①初始矫直力为：式中：r0c——矫直件的初始曲率半径，y0s——初始状态弹性厚度；②弹性区矫直力为：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)]]>式中：yis——经过第i辊时钢板压缩区的弹性厚度，当压缩区在钢板中心层下方时取负，压缩区在钢板中心层上方时取正，ric——经过第i辊时矫直件的曲率半径，f1(ε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)]]>式中：f2(εs)——反向屈服应力减少量；<3>利用自主开发的三维弹塑性多物体接触多极边界元法程序包计算矫直力利用fortran语言开发三维弹塑性多物体接触多极边界元法程序包，用于模拟辊式矫直过程；在矫直过程中，将轧件变形视为板弯曲问题的反问题，先建立板弯曲问题的多极边界元法的边界积分方程，利用板弯曲问题的反问题求解矫直力，具体步骤如下：设板的边界为Γ，内部区域为Ω，用w表示板中面的挠度，则板弯曲问题的控制方程可用w做变量表示为：-D▿4w+F=0]]>——板的弯曲刚度，w——挠度，E——材料的弹性模量，ν——泊松比，F——垂直作用于板面的压力，板的边界条件为：在Γ1上：在Γ2上：βn——板边界上任意点处绕法线的转角，βs——板边界上任意点处绕切线的转角，q——剪力，mn——板截面上的弯矩，ms——板截面上的向扭矩，mn=n1m1+n2m2=m11n12+2m12n1n2+m22n22]]>ms=-n2m1+n1m2=-(m11-m22)n1n2+m12(n12-n22)]]>q＝n1q1+n2q2，边界积分方程可写为：F3∫Ωw*dΩ=Ciwi+∫Γ(mn*βn+t*w)dΓ+Σk=1Ltkc*wk-∫Γ(mnβn*+tw*)dΓ-Σk=1Ltkcwk*---(12)]]>w*——单位集中力沿x3轴方向作用于无穷大薄板的ξ点时，板的几何中心层上任意一点产生的挠度，r表示源点与场点之间的距离，q1*=∂m11*∂x1+∂m12*∂x2,q2*=∂m12*∂x1+∂m22*∂x2,]]>Σk=1Ltkc*wk=∫Γ(∂ms*∂sw+ms*∂w∂s)dΓ,Σk=1Ltkcwk*=∫Γ(∂ms∂sw*+ms∂w*∂s)dΓ,]]>计算流程如下所述:①将矫直辊和矫直件划分单元,将单元和节点信息以文本形式读入，并确定边界条件；②输入增量步数据并根据接触条件进行多物体接触检查；③定义摩擦接触约束；④求解单元最大长度，根据单元长度分别采用指数展开和球谐函数展开求多极矩系数；⑤利用GMRES(m)求解方程组(12)，并检验是否满足精度，若满足则执行步骤⑥，否则返回步骤③改变单元长度重新计算；⑥计算表面应力，判断摩擦状态是否改变，若摩擦状态改变则返回步骤④重新计算，否则执行步骤⑦；⑦计算内点应力，判断是否发生穿透，若发生穿透则返回步骤③，细分增量步，否则执行步骤⑧；⑧计算是否为最后增量步，是则结束计算，否则返回步骤②继续迭代计算；<4>综合考虑多极边界元法计算的数值结果，步骤<1>中的中性层偏移量计算公式(1)和步骤<2>中的矫直力计算公式(9)、(10)、(11)，最终得到矫直力计算公式：F＝α1F1+α2F2+α3F3         (13)式中：α1，α2，α3——权重系数，取决于具体工况，F1——由中性层偏移决定的矫直力，F2——由包辛格效应决定的矫直力，F3——由多极边界元法计算得到的矫直力。