2. 长春设备工艺研究所, 长春 130012
2. Changchun Equipment Technology Research Institute, Changchun 130012, China
大型对轮强力旋压机床用来制造大直径薄壁空心筒形件,主要是指直径≥2 000 mm、壁厚≤10 mm、长度≥3 000 mm的零件[1],随着我国综合国力的不断提升,我国军工领域、航天航空领域对于该类产品的需求越来越迫切,近些年我国在大型数控强力旋压机床设备设计制造、工艺研究领域发展迅速,并且取得了较大进展。
Luo等[2]采用有限元分析方法研究了旋压机几种不同进给结构的动态性能和静态性能。张健等[3]针对大型超薄筒形件减薄旋压成形时易出现旋压失稳的缺陷,研究了一种全新的同步旋转模环旋压工艺。Wang等[4]通过旋压试验和有限元分析研究了4种不同旋轮进给轮廓对旋压力的影响,在常规金属旋压过程中数值模拟分析了零件壁厚和应力变化。Hartel等[5]通过数值模拟仿真的方法优化了非圆旋压过程。Xu等[6]研究了在管型件旋压过程中不同旋轮分布情况下旋压力和毛坯变形的特点。Wang等[7]分析了3种不同毛坯旋压加工带有内部加强筋的薄壁壳体的成形过程。Xiao等[8]研究了筒形零件旋压过程中等效应变的分布情况,以实现对旋压件沿厚度方向的整体晶粒细化。胡玄通等[9]根据双轮旋压机的工作特点以及相关的技术要求对双轮旋压机的机械部分进行了设计。黄涌等[10]根据旋压成形原理及对轮旋压工艺特点,设计了一种筒形件强力旋压用对轮旋压的成形装置。邹小堤等[11]研究了影响旋压机关键零件主轴的制造加工精度和质量的一些重要因素。
以上研究使我国旋压成形加工技术得到了持续发展与进步,在旋压机结构设计方面、基础理论与旋压工艺研究方面以及工程应用方面都有很大的发展,但是在结构设计时多数采用类比或者选择较大的安全系数方法进行设计,对旋压设备整体或零部件的可靠性研究还没有涉及。
有鉴于此,本文结合最大熵方法、调用结构有限元模型和Taylor展开方法进行结构可靠性分析。以某对轮强力旋压设备旋轮座基体为研究对象,对旋轮座结构工况进行受力分析,建立了旋轮座基体功能函数;采用最大熵方法得到概率密度函数表达式,利用Taylor展开方法与调用有限元程序相结合获得功能函数前四阶中心矩,对于通过调用有限元程序所获得参数进行实验验证;将功能函数标准化得到四阶原点矩,采用MATLAB函数求解最大熵方法约束条件表达式,得到最大熵概率密度函数所需的参数,拟合概率密度曲线;采用积分求解功能函数的失效概率,计算结果与MCS方法对比证明了本文方法的正确性。
1 大型对轮强力旋压设备结构和工作原理某对轮强力旋压设备结构如图 1所示,主要动作包括旋轮横向进给动作、纵向进给动作、主轴旋转动作,横向进给动作采用4对轮在圆周上4等分,相互间隔90°,内、外2组旋轮一组在同一轴线上对顶,其中4个内轮安装在内部旋轮座装置上,4个外轮分别安装在2个导向立柱和2个支撑立柱上,横向进给动作由油缸驱动,采用电液伺服闭环控制。纵向采用双丝杠传动,2台伺服电机同步驱动纵向滑体实现纵向进给,纵向滑体沿4根导轨运动,4根导轨分别安装在2个导向立柱上。纵向滑体上端安装一个转动装置,包括齿轮和轴承,叫做回转支撑轴承,作为机床的主轴,主轴转动采用伺服电机驱动。设备主要参数如表 1所示。
![]() |
图 1 (网络版彩图)对轮旋压设备结构 |
该设备用于大直径筒形件减薄旋压,毛坯为筒形件,工作原理如图 2所示。工作时,毛坯安装在回转支撑轴承上,毛坯随主轴旋转同时做纵向进给运动,内旋轮和外旋轮横向进给,四对轮分别单独位置闭环控制,压下量相同,旋轮与毛坯点接触对毛坯施加压力达到减薄毛坯的目的,内、外旋轮座承受着纵向和横向力共同作用。
![]() |
图 2 对轮旋压工作原理 |
对轮旋压采用内旋轮代替芯模对工件进行加工,内、外旋轮在成形过程中同时从内、外表面对工件进行加工,能够大大减小旋压力,提高工件内表面精度及机床加工能力,同时对于大直径筒形件的加工,节省了芯模,减少了加工成本。
2 旋轮座基体功能函数内、外旋轮座是对轮旋压机的主要受力部件,工作时旋轮与毛坯接触,承受横向力、纵向力和切向力的复合作用,其中横向力最大,并且随着压下量的增大而增大,本文以外旋轮座为研究对象,其结构如图 3所示。
![]() |
图 3 外旋轮座结构 |
如图 3所示,滑体和旋轮座基体之间采用铜板,滑动摩擦实现进给运动,由液压缸驱动,滑体在旋轮座基体内运动。
工作时,旋压力作用在旋轮上,F1为横向力,F2为纵向力,旋压力作用在旋轮上,通过旋轮轴、滑体、活塞杆、后压板和上压板传递到旋轮座基体上。旋轮座基体结构如图 4所示,材料为ZG270-450,弹性模量E=1.75×1011 Pa,长度L=955 mm,宽度B=453 mm,高度H=453 mm,壁厚δ=65 mm。
![]() |
图 4 旋轮座基体结构图 |
旋轮座基体底部通过螺钉固定,在横向力F1和纵向力F2的复合作用下,旋轮座基体会有退让变形,旋轮座基体退让变形量是影响旋压产品质量和精度的主要因素,以旋轮座基体变形建立功能函数如下:
$ Z=d_{\max }-d\left(F_{1}, F_{2}\right). $ | (1) |
其中:dmax表示旋轮座基体在保证加工产品质量可靠的情况下所允许的最大变形,根据工作经验取dmax=0.05 mm;d(F1, F2)表示旋轮座基体实际工作时的变形,F1和F2为随机变量,单位为kN,该参数根据某典型筒形件旋压过程中实际测得的数值统计得到,两者的统计特征如表 2所示。Z>0代表结构可靠,Z < 0代表结构失效。
随机变量 | 均值 | 标准差 | 三阶矩 | 四阶矩 |
F1/kN | 298.99 | 12.65 | -1 390 | 1.76×105 |
F2/kN | 523.38 | 96.07 | 858.07 | 1.37×107 |
3 最大熵法
功能函数(1)若服从概率密度函数为fz(Z)的连续分布,则定义熵为[12]
$ H=-c \int_{-\infty}^{+\infty} f_{z}(Z) \ln f_{z}(Z) \mathrm{d} Z. $ | (2) |
其中参数c>0为常数。考虑将功能函数Z的前m阶原点矩υZi(i=0, 1, …, m)作为约束条件,即在条件:
$ \upsilon_{Z_{i}}=E\left[(Z)^{i}\right]=\int_{-\infty}^{+\infty}(Z)^{i} f_{z}(Z) \mathrm{d} Z. $ | (3) |
得到如下优化:
$ \begin{gathered} \max H=-c \int_{-\infty}^{+\infty} f_{z}(Z) \ln f_{z}(Z) \mathrm{d} Z, \\ \text { s. t. } \upsilon_{Z_{i}}=E\left[(Z)^{i}\right]=\int_{-\infty}^{+\infty}(Z)^{i} f_{z}(Z) \mathrm{d} Z. \end{gathered} $ | (4) |
利用Lagrange乘子法,可得最大熵概率密度函数为
$ f_{z}(Z)=\exp \left(-\sum\limits_{i=0}^{m} a_{i}(Z)^{i}\right). $ | (5) |
其中:a0,a1,…,am为待定参数,为求待定参数,首先求功能函数的前m阶原点矩。
3.1 功能函数统计矩功能函数(1)中,令x =[F1, F2],则功能函数变为Z=g(x)=dmax-d(F1,F2),x的均值为μx,第k个变量的二、三、四阶中心矩分别为μxk2、μxk3、μxk4(k =1, 2),将功能函数在x*处作Taylor展开得到[13-14]:
$ \begin{gathered} Z_{Q} \cong g\left(x^{*}\right)+\left(x-x^{*}\right)^{\mathrm{T}} \nabla g\left(x^{*}\right)+ \\ \frac{1}{2}\left(x-x^{*}\right)^{\mathrm{T}} \nabla^{2} g\left(x^{*}\right)\left(x-x^{*}\right). \end{gathered} $ | (6) |
可近似计算ZQ的前四阶矩如下:
$ \mu_{Z_{Q 1}}=E\left(Z_{Q}\right)=g\left(\mu_{x}\right)+\frac{1}{2} \sum\limits_{k=1}^{2} \frac{\partial^{2} g\left(\mu_{x}\right)}{\partial x_{k}^{2}} \mu_{x k 2}. $ | (7) |
$ \begin{array}{c} \mu_{Z_{Q 2}}=E\left[\left(Z_{Q}-\mu_{Z_{Q 1}}\right)^{2}\right]=\\ \sum\limits_{k=1}^{2}\left[\frac{\partial g\left(\mu_{x}\right)}{\partial x_{k}}\right]^{2} \mu_{x k 2}+\sum\limits_{k=1}^{2} \frac{\partial g\left(\mu_{x}\right)}{\partial x_{k}} \frac{\partial^{2} g\left(\mu_{x}\right)}{\partial x_{k}^{2}} \mu_{x k 3}+ \\ \frac{1}{4} \sum\limits_{k=1}^{2}\left[\frac{\partial^{2} g\left(\mu_{x}\right)}{\partial x_{k}^{2}}\right]^{2}\left(\mu_{x k 4}-3 \mu_{x k 2}^{2}\right)+ \\ \frac{1}{2} \sum\limits_{k=1}^{2} \sum\limits_{o=1}^{2}\left[\frac{\partial^{2} g\left(\mu_{x}\right)}{\partial x_{k} \partial x_{o}}\right]^{2} \mu_{x k 2} \mu_{x o 2}, \end{array} $ | (8) |
$ \begin{array}{c} \mu_{Z_{Q 3}}=E\left[\left(Z_{Q}-\mu_{Z_{Q 1}}\right)^{3}\right]=\sum\limits_{k=1}^{2}\left[\frac{\partial g\left(\mu_{x}\right)}{\partial x_{k}}\right]^{3} \mu_{x k 3}+ \\ \frac{3}{2} \sum\limits_{k=1}^{2}\left[\frac{\partial g\left(\mu_{x}\right)}{\partial x_{k}}\right]^{2} \frac{\partial^{2} g\left(\mu_{x}\right)}{\partial x_{k}^{2}}\left(\mu_{x k 4}-3 \mu_{x k 2}^{2}\right)+ \\ 3 \sum\limits_{k=1}^{2} \sum\limits_{o=1}^{2} \frac{\partial g\left(\mu_{x}\right)}{\partial x_{k}} \frac{\partial g\left(\mu_{x}\right)}{\partial x_{o}} \frac{\partial^{2} g\left(\mu_{x}\right)}{\partial x_{k} \partial x_{o}} \mu_{x k 2} \mu_{x o 2}, \end{array} $ | (9) |
$ \begin{gathered} \mu_{Z_{Q 4}}=E\left[\left(Z_{Q}-\mu_{Z_{Q 1}}\right)^{4}\right]= \\ \sum\limits_{k=1}^{2}\left[\frac{\partial g\left(\mu_{x}\right)}{\partial x_{k}}\right]^{4}\left(\mu_{x k 4}-3 \mu_{x k 2}^{2}\right)+ \\ 3 \sum\limits_{k=1}^{2} \sum\limits_{o=1}^{2}\left[\frac{\partial g\left(\mu_{x}\right)}{\partial x_{k}}\right]^{2}\left[\frac{\partial g\left(\mu_{x}\right)}{\partial x_{o}}\right]^{2} \mu_{x k 2} \mu_{x o 2} . \end{gathered} $ | (10) |
将ZQ转换为标准正态随机变量:
$ Z_{Y}=\frac{Z_{Q}-\mu_{Z_{Q 1}}}{\sigma_{Z_{Q}}}. $ | (11) |
则根据中心矩和原点矩的定义,得到ZQ和ZY的各阶中心矩的关系为
$ \begin{gathered} \mu_{Z_{Q i}}=E\left[\left(Z_{Q}-\mu_{Z_{Q 1}}\right)^{i}\right]=E\left[\left(\sigma_{Z_{Q}}^{i} Z_{Y}\right)^{i}\right]= \\ \sigma_{Z_{Q}}^{i} E\left[\left(Z_{Y}\right)^{i}\right]=\sigma_{Z_{Q}}^{i} \mu_{Z_{Y i}}= \\ \sigma_{Z_{Q}}^{i} \upsilon_{Z_{Y i}}, \quad i=0, 1, \cdots, m. \end{gathered} $ | (12) |
其中
功能函数标准化后,各阶原点矩和中心矩相同,标准正态功能函数四阶矩统计信息如下:
$ \mu_{Z_{Y 0}}=\upsilon_{Z_{Y 0}}=1, $ | (13) |
$ \mu_{Z_{Y 1}}=\upsilon_{Z_{Y 1}}=0, $ | (14) |
$ \mu_{Z_{Y 2}}=\upsilon_{Z_{Y 2}}=1, $ | (15) |
$ \mu_{Z_{Y 3}}=\upsilon_{Z_{Y 3}}=\frac{\mu_{Z_{Q 3}}}{\sigma_{Z_{Q}}^{3}}, $ | (16) |
$ \mu_{Z_{Y 4}}=\upsilon_{Z_{Y 4}}=\frac{\mu_{Z_{Q 4}}}{\sigma_{Z_{Q}}^{4}}. $ | (17) |
功能函数标准化后,最大熵概率密度函数仍然是式(5)的形式,表示如下:
$ f_{Z}\left(Z_{Y}\right)=\exp \left(-\sum\limits_{i=0}^{m} a_{i}\left(Z_{Y}\right)^{i}\right). $ | (18) |
节3.1中得到功能函数的前四阶原点矩,所以式(18)中m=4,将式(13)—(18)代入公式(3)得到如下方程:
$ \int_{-\infty}^{+\infty}\left(Z_{Y}\right)^{i} \exp \left(-\sum\limits_{i=0}^{m} a_{i}\left(Z_{Y}\right)^{i}\right) \mathrm{d} Z_{Y}=\upsilon_{Z_{Y i}} . $ | (19) |
采用MATLAB求解方程(19),得到参数ai,将ai代入式(18)得到概率密度函数,然后积分求解失效概率:
$ \begin{gathered} p_{f}={Pr}\left(Z_{Q} \leqslant 0\right)={Pr}\left(Z_{Y} \leqslant-\frac{\mu_{Z_{Q 1}}}{\sigma_{Z_{Q}}}\right)= \\ \int_{-\infty}^\frac{\mu_{Z_{Q 1}}}{\sigma_{Z_{Q}}} \exp \left(-\sum\limits_{i=0}^{m} a_{i} Z_{Y}\right) \mathrm{d} Z_{Y}. \end{gathered} $ | (20) |
由于旋轮座基体功能函数是隐式函数,无法采用数值计算方法直接求解功能函数的统计矩,因此首先建立旋轮座基体三维模型,然后导入有限元软件,划分网格,建立了旋轮座基体的有限元模型,如图 5所示,对旋轮座基体底面上的16个连接孔施加固定约束,对整个底面施加固定约束,对旋轮座基体与后压板连接面施加载荷F1,对旋轮座基体与上压板连接面施加载荷F2。
![]() |
图 5 (网络版彩图)旋轮座基体有限元模型 |
图 6为某典型筒形件旋压实验,横向力F1通过加装在横向进给油缸两腔的压力传感器反馈电流信号经计算机计算得到,纵向力F2通过加装在丝杠端部的力传感器反馈电流信号经计算机计算得到,由计算机显示并记录力的数值。
![]() |
图 6 (网络版彩图)某典型筒形件旋压实验 |
旋轮座基体工作时的变形可以通过千分表测量得到。如图 7所示,千分表安装在旋轮座基体变形最大位置,千分表数值由人工读取,读取时间间隔为120 s/次。
![]() |
图 7 (网络版彩图)旋轮座基体变形测量 |
实验所用筒形件毛坯长度960 mm,纵向进给速度40 mm/min,主轴转速10 r/min。根据纵向进给速度和毛坯长度可以得到所需的旋压时间为24 min,由于毛坯需要装夹,不能全部旋压,因此旋压时间取22 min,则一共读取千分表数值11个,用δ1表示,而有限元计算值用δ2表示,两者的对比见表 3。
序号 | F1/kN | F2/kN | δ1/mm | δ2/mm | 相对误差/% |
1 | 257.98 | 502.18 | 0.042 | 0.039 5 | 5.95 |
2 | 267.89 | 512.33 | 0.043 | 0.039 8 | 7.44 |
3 | 302.99 | 522.55 | 0.045 | 0.043 3 | 3.78 |
4 | 296.32 | 522.88 | 0.044 | 0.042 9 | 2.50 |
5 | 304.78 | 522.98 | 0.045 | 0.043 7 | 2.89 |
6 | 297.40 | 523.19 | 0.046 | 0.044 3 | 3.70 |
7 | 292.94 | 523.51 | 0.046 | 0.044 1 | 4.13 |
8 | 307.02 | 522.81 | 0.047 | 0.045 4 | 3.40 |
9 | 300.26 | 524.25 | 0.047 | 0.0453 | 3.62 |
10 | 284.49 | 524.46 | 0.048 | 0.0469 | 2.29 |
11 | 306.82 | 524.99 | 0.048 | 0.0474 | 1.25 |
从表 3中可以看出,千分表读取数值与调用ANSYS有限元程序所得数值,两者的最大相对误差为7.44%,平均相对误差为3.72%,通过实验验证了有限元模型的正确性。
4.2 旋轮座基体失效概率借助于ANSYS有限元分析软件,调用有限元程序得到旋轮座基体的位移值,然后结合式(7)—(10)计算功能函数统计矩,计算得到功能函数统计矩结果列于表 4中。
中心矩 | MCS | 本文方法 | 相对误差/% |
一阶 | 5.293×10-4 | 5.290×10-4 | 0.056 |
二阶 | 3.926 4×10-8 | 3.893 1×10-8 | 0.85 |
三阶 | -9.527 3×10-13 | -9.509 4×10-13 | 1.88 |
四阶 | 4.621 3×10-15 | 4.6071×10-15 | 3.07 |
根据统计矩,拟合功能函数概率密度曲线,最后通过积分计算失效概率,在较短运行时间的情况下,本文计算所得功能函数失效概率为0.005 43。
本文将失效概率结果与MCS对比。由于功能函数为隐式,本文采用ANSYS软件中的PDS模块进行有限元分析,采用MCS直接抽样法,抽样次数为105次,MCS方法求出的失效概率为0.005 58,本文方法与MCS方法相对误差仅为1.61%,但是本文方法计算时间仅需0.869 3 s,而MCS方法所需计算时间为46 752.82 s,图 8为本文方法与MCS方法的概率密度曲线图,可以看出吻合度较好。
![]() |
图 8 概率密度曲线 |
通过前面的计算得到旋轮座在工作状态下的可靠度,可以通过改进设计旋轮座基体结构来提高可靠度,在旋轮座基体两边分别加3道加强筋,有限元模型如图 9所示。
![]() |
图 9 (网络版彩图)旋轮座基体改进设计有限元模型 |
改进设计前、后的可靠度指标如表 5所示。
从表 5中可以看出,改进后旋轮座基体结构的失效概率降低11.42%,证明改进后可靠性显著提高,将改进设计的结果量化给设计者提供更直观的参考。
5 结论本文分析了对轮旋压设备旋轮座结构的受力工况,提出旋轮座结构可靠性分析方法。建立了旋轮座基体有限元模型,调用ANSYS有限元程序、Taylor展开方法求解功能函数前四阶中心矩,最后用最大熵方法拟合功能函数概率密度曲线,积分概率密度函数求解失效概率,所得结果与Monte Carlo方法结果对比证明所提出的方法的正确性和计算效率高。本文还对旋轮座结构进行了改进设计,对比改进前、后的失效概率,改进后失效概率降低了11.42%。
本文通过对某典型筒形件的实验研究,采集了旋轮座基体功能函数所需的变形参数,与调用ANSYS有限元程序所得变形量进行对比,验证了有限元分析的正确性。
本文提出的方法对大型旋压设备关键零件的设计分析方面做出了有益的探索,取得了满意的效果,获得了宝贵而丰富的经验,为更大型旋压机的设计和分析提供了借鉴。
[1] |
王大力, 郭亚明, 李亦楠, 等. 大型薄壁筒形件对轮旋压成形数值模拟及成形精度分析[J]. 锻压技术, 2020, 45(3): 47-54. WANG D L, GUO Y M, LI Y N, et al. Numerical simulation and forming precision analysis on counter-roller spinning for large thin-walled cylindrical parts[J]. Forging & Stamping Technology, 2020, 45(3): 47-54. (in Chinese) |
[2] |
LUO W, LU B, ZHENG H W, et al. A research of spinning machine: Effect of different structures of feed system on the static and dynamic characteristics[J]. Advance in Mechanical Engineering, 2019, 11(2): 1-15. |
[3] |
张健, 庞俊忠, 张洪瑞, 等. 同步旋转模环旋压成形的数值模拟分析[J]. 铸造技术, 2016, 37(2): 351-354. ZHANG J, PANG J C, ZHANG H R, et al. Numerical Simulation analysis of synchronous rotary mould ring spinning[J]. Foundry Technology, 2016, 37(2): 351-354. (in Chinese) |
[4] |
WANG L, LONG H. A study of effects of roller path profiles on tool forces and part wall thickness variation in conventional metal spinning[J]. Journal of Materials Processing Technology, 2011, 211(12): 2140-2151. DOI:10.1016/j.jmatprotec.2011.07.013 |
[5] |
HARTEL S, LAUE R. An optimization approach in non-circular spinning[J]. Journal of Materials Processing Technology, 2016, 229: 417-430. DOI:10.1016/j.jmatprotec.2015.09.003 |
[6] |
XU W C, ZHAO X K, MA H, et al. Influence of roller distribution modes on spinning force during tube spinning[J]. International Journal of Mechanical Sciences, 2016, 113: 10-25. DOI:10.1016/j.ijmecsci.2016.04.009 |
[7] |
WANG Z, MA S. Analysis of thin-walled shells with inner ribs formed by inner spinning technology[J]. Materials Research Innovations, 2015, 19(5): 101-105. |
[8] |
XIAO G F, XIA Q X, CHENG X Q, et al. Metal flow model of cylindrical parts by counter-roller spinning[C]//11th International Conference on Technology of Plasticity. Nagoya, Japan: Procedia Engineering, 2014: 2397-2402.
|
[9] |
胡玄通, 夏琴香, 程秀全. HGPX-WSM卧式双旋轮数控旋压机结构设计[J]. 锻压技术, 2014, 39(8): 59-64. HU X T, XIA Q X, CHEN X Q. Structural design of HGPX-WSM horizontal type CNC spinning machine with double rollers[J]. Forging & Stamping Technology, 2014, 39(8): 59-64. (in Chinese) |
[10] |
黄涌, 夏琴香, 程秀全, 等. 筒形件强力旋压用对轮旋压装置的研制[J]. 锻压技术, 2013, 38(6): 62-72. HUANG Y, XIA Q X, CHENG X Q, et al. Developing of counter roller spinning equipment for tube power spinning[J]. Forging & Stamping Technology, 2013, 38(6): 62-72. (in Chinese) |
[11] |
邹小堤, 黄照, 桂林, 等. 大型数控强力旋压机关键零件制造工艺研究[J]. 武汉理工大学学报, 2012, 36(4): 821-824. ZOU X D, HUANG Z, GUI L, et al. Manufacturing process research on the key parts of large NC power spinning machine[J]. Journal of Wuhan University of Technology, 2012, 36(4): 821-824. DOI:10.3963/j.issn.2095-3844.2012.04.039 (in Chinese) |
[12] |
贡金鑫. 工程结构可靠度计算方法[M]. 大连: 大连理工大学出版社, 2003. GONG J X. Engineering structure reliability calculation method[M]. Dalian: Dalian University of Technology Press, 2003. (in Chinese) |
[13] |
张明. 结构可靠度分析[M]. 北京: 科学出版社, 2009. ZHANG M. Structural reliability analysis[M]. Beijing: Science Press, 2009. (in Chinese) |
[14] |
ZHENG H W, MENG G W, LI F, et al. Hybrid structure reliability analysis based on the damped Newton method[J]. Tsinghua Science and Technology, 2020, 25(5): 668-677. DOI:10.26599/TST.2019.9010073 |