2017, Vol. 49引用本文 |
刚度是机械系统设计的重要环节，直接影响系统定位精度和动态特性，国内外学者对其做了大量研究，并取得了显著成果[1-4]。于靖军等[5-6]利用结构力学的方法建立了一种柔性并联机构的刚度矩阵，该解析法可以得到结构规则的任意刚度构件的解析表达式。付铁等[7]提出一种基于静刚度的并联机床误差分析方法，用于提高机床精度。李育文等[8]利用有限元分析了6-UPS并联机床在工作空间内的静刚度性能。
传统刚度分析方法只能得到规则结构的刚度解析表达式，而当构件的结构不规则且非常复杂时，往往由于简化导致刚度矩阵精确度下降。因此，为了便于进行面向刚度的机械系统参数优化设计，必须解决不规则构件的规范化刚度矩阵求解问题。作者提出一种基于多元线性回归原理求解复杂构件刚度矩阵的建模方法。其基本思想是，通过ANSYS软件选取力和位移多样本数对，然后基于多元线性回归理论求解刚度矩阵并进行有效性检验，实现解决复杂构件刚度矩阵的求解难题。
1 刚度矩阵多元线性回归分析一般意义上刚度是抵抗外界变形的能力，空间机构刚度矩阵可由式(1)表述:
| $ \mathit{\boldsymbol{F}} = \mathit{\boldsymbol{KX}} $ | (1) |
式中,
显然，式(1)共包含6个线性方程，令Ki为K矩阵中第i行的刚度向量，则
| $ {\mathit{\boldsymbol{K}}_i} = \left[ {\begin{array}{*{20}{c}} {{k_{i1}}}&{{k_{i2}} \ldots }&{{k_{i6}}} \end{array}} \right],{\rm{ }}i = 1,2, \ldots ,6。$ |
以K矩阵第1行方程为例，其方程为：
| $ {F_x} = {k_{11}}{\Delta _x} + {k_{12}}{\Delta _y} + {k_{13}}{\Delta _z} + {k_{14}}{\Delta _{\theta x}} + {k_{15}}{\Delta _{\theta y}} + {k_{16}}{\Delta _{\theta z}} $ | (2) |
基于ANASYS软件分析，可以获得广义外力及相应特定点、线、面的广义位移数对样本，假定选取n组观测代入式(2)得：
| $ \left\{ \begin{array}{l} {F_{1x}} = {k_{11}}{\Delta _{1x}} + {k_{12}}{\Delta _{1y}} + {k_{13}}{\Delta _{1z}} + {k_{14}}{\Delta _{1\theta x}} + \\ \;\;\;\;\;\;\;\;{k_{15}}{\Delta _{1\theta y}} + {k_{16}}{\Delta _{1\theta z}},\\ {F_{2x}} = {k_{11}}{\Delta _{2x}} + {k_{12}}{\Delta _{2y}} + {k_{13}}{\Delta _{2z}} + {k_{14}}{\Delta _{2\theta x}} + \\ \;\;\;\;\;\;\;\;{k_{15}}{\Delta _{2\theta y}} + {k_{16}}{\Delta _{2\theta z}},\\ \cdots \\ {F_{nx}} = {k_{11}}{\Delta _{nx}} + {k_{12}}{\Delta _{ny}} + {k_{13}}{\Delta _{nz}} + {k_{14}}{\Delta _{n\theta x}} + \\ \;\;\;\;\;\;\;\;{k_{15}}{\Delta _{n\theta y}} + {k_{16}}{\Delta _{n\theta z}} \end{array} \right. $ | (2) |
可以进一步简化为矩阵形式：
| $ \mathit{\boldsymbol{Y}} = \Delta \mathit{\boldsymbol{XK}}_1^{\rm{T}} $ | (3) |
| $ \mathit{\boldsymbol{Y}} = {\left[ {\begin{array}{*{20}{c}} {{F_{1x}}}&{{F_{2x}}}& \cdots &{{F_{nx}}} \end{array}} \right]^{\rm{T}}} $ | (4) |
| $ \Delta \mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{c}} {{\Delta _{1x}}}&{{\Delta _{1y}}}&{{\Delta _{1z}}}&{{\Delta _{1\theta x}}}&{{\Delta _{1\theta y}}}&{{\Delta _{1\theta z}}}\\ {{\Delta _{2x}}}&{{\Delta _{2y}}}&{{\Delta _{2z}}}&{{\Delta _{2\theta x}}}&{{\Delta _{2\theta y}}}&{{\Delta _{2\theta z}}}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {{\Delta _{nx}}}&{{\Delta _{ny}}}&{{\Delta _{nz}}}&{{\Delta _{n\theta x}}}&{{\Delta _{n\theta y}}}&{{\Delta _{n\theta z}}} \end{array}} \right] $ | (5) |
式中，Y为线性回归模型的因变量，末端位移ΔX为自变量，K1T为回归系数。根据力和位移的多组数对样本，利用多元线性回归分析可以求解得到K1T。再通过对回归系数应用t检验进一步化简方程形式，并对简化后的方程进行F检验和拟合优度检验，验证方程的回归效果，最终确定K1。多元线性回归分析流程如图 1所示。
 |
| 图1 刚度矩阵多元线性回归分析流程图 Fig. 1 Flow chart of rigid matrix analysis based on multiple linear regression method |
多元线性回归分析通过最小化误差的平方和寻找数据的最佳函数匹配，有效提高了分析效率和准确度[9]。由式(3)得回归模型中未知参数k11, k12, …, k16的估计，依据最小二乘原理，目标是寻找回归参数
| $ \begin{array}{l} M(\mathop {{k_{11}}}\limits^ \wedge ,\mathop {{k_{12}}}\limits^ \wedge , \cdots ,\mathop {{k_{16}}}\limits^ \wedge ) = \sum\limits_{i = 1}^n {\left( {{{({y_i} - \mathop {{k_{11}}}\limits^ \wedge {x_{i1}} - \mathop {{k_{12}}}\limits^ \wedge {x_{i2}} - \ldots - \mathop {{k_{1n}}}\limits^ \wedge {x_{in}})}^2}} \right)} = \\ \mathop {\min }\limits_{{k_{11}},{k_{12}}, \cdots ,{k_{16}}} 1\sum\limits_{i = 1}^n {{{({y_i} - {k_{11}}{x_{i1}} - {k_{12}}{x_{i2}} - \ldots - {k_{1n}}{x_{in}})}^2}} \end{array} $ | (6) |
通过多元线性回归方法求解完备的6×6空间刚度矩阵需要确定36个回归系数。由于该刚度矩阵主要应用于后续的刚度矩阵叠加，为降低计算难度，提高计算效率，在保证对Fx进行预报和控制质量的前提下，对回归方程中的每一个回归系数作显著性检验[10-11]，剔除不显著变量，重新建立更简单、更精确的线性回归方程就很有必要了。
为了验证刚度矩阵中每个元素的显著性，首先进行t检验，由t分布定义可知:
| $ {t_\alpha } = \frac{{{K_{1j}}}}{{\sqrt {{c_{jj}}} \sqrt {\frac{{{Q_{{\rm{res}}}}}}{{n - m - 1}}} }} $ | (7) |
式中, Qres指剩余离差，cjj为矩阵C的对角线元素，C为ANASYS分析实验中获得的n组观测数据对应的矩阵求逆，m为回归系数的个数，K1j为K的第j个元素(j=1，2，…，m)。在检验过程中，选取显著性水平后，将临界值和由式(7)计算的值进行对比，将不符合条件的元素剔除，简化刚度矩阵。
由于t检验后刚度模型会发生形式改变，需要对其继续进行F检验，以确保刚度模型的显著性要求。由F分布的定义可知：
| $ F = \frac{{{Q_{{\rm{res}}}}/m}}{{{Q_{{\rm{res}}}}/\left( {n - m - 1} \right)}} $ | (8) |
式中, Qreg为回归离差，在进行检验时取定显著性水平，将查表确定的F临界值和式(8)计算值对比，即可确定回归方程的显著性。
最后，为了进一步验证所得模型的精确程度，需进行拟合优度检验，其复相关系数R2可以表示为:
| $ {R^2} = \frac{{{Q_{{\rm{res}}}}}}{{{Q_{\rm{T}}}}} = \frac{{{Q_{\rm{T}}} - {Q_{{\rm{res}}}}}}{{{Q_{\rm{T}}}}} $ | (9) |
式中, QT为总离差，复相关系数的取值范围为[0, 1]，R2越接近1，说明回归直线对观测值的拟合程度越好，即所得的刚度矩阵越精确。
3 数值算例图 2为3-URS并联机构，由3条完全相同的支链组成。每条支链均包括U型轴、上连杆、下连杆以及球铰，最终末端平台可实现空间6自由度运动。
 |
| 1. 动平台；2. 上连杆；3. 滑动槽；4. 下连杆；5. U型轴。 图2 3-URS并联机构 Fig. 2 3-URS parallel mechanism |
3-URS并联机构中球铰、动平台、底座支架等均为规则组件，可以采用传统方法建立刚度模型。而上下连杆、滑动槽、U型轴等构件结构复杂，均为不规则形状。
图 3为U型轴的3维模型，显然很难用普通方法求解其刚度矩阵，进一步导致无法获得3-URS并联机构整体刚度的解析表达式。
 |
| 图3 U型轴不规则构件 Fig. 3 U-shape axis irregular component |
根据图 2所示3-URS并联机构的结构特点和受力特征，确定U型轴的输入力在A、B两圆周面上，C面为输出端，在C面上建立如图 3所示的坐标系。然后在C面施加多组广义外力并通过ANASYS软件测量该面上的广义变形值，选择24次样本实验，如表 1所示。
| 表1 变形量和外力样本数据 Tab. 1 Datum of the force and deformation |
 |
将表 1的数据分别代入式(6)，对其求极值，可分别求解得到6组线性方程的回归系数，整理后得到刚度矩阵K。
| $ \begin{array}{l} \mathit{\boldsymbol{K}} = \\ \left[ {\begin{array}{*{20}{c}} {5.338{\rm{ }}78 \times {{10}^8}}&{7.237{\rm{ }}12 \times {{10}^4}}&{ - 1.222{\rm{ }}14 \times {{10}^5}}&{7.927{\rm{ }}35 \times {{10}^2}}&{2.187{\rm{ }}72 \times {{10}^2}}&{2.182{\rm{ }}18 \times {{10}^5}}\\ { - 1.288{\rm{ }}96 \times {{10}^5}}&{4.382{\rm{ }}29 \times {{10}^7}}&{ - 1.099{\rm{ }}33 \times {{10}^5}}&{1.274{\rm{ }}74 \times {{10}^3}}&{4.704{\rm{ }}96 \times {{10}^2}}&{8.741{\rm{ }}35 \times {{10}^4}}\\ {1.635{\rm{ }}09 \times {{10}^5}}&{ - 2.815{\rm{ }}65 \times {{10}^4}}&{1.739{\rm{ }}04 \times {{10}^8}}&{ - 1.474{\rm{ }}37 \times {{10}^6}}&{ - 4.801{\rm{ }}19 \times {{10}^5}}&{5.169{\rm{ }}35 \times {{10}^2}}\\ { - 2.963{\rm{ }}54 \times {{10}^5}}&{9.118{\rm{ }}17 \times {{10}^2}}&{ - 2.681{\rm{ }}78 \times {{10}^6}}&{3.004{\rm{ }}08 \times {{10}^4}}&{8.196{\rm{ }}44 \times {{10}^3}}&{ - 7.069{\rm{ }}42}\\ { - 7.773{\rm{ }}49 \times {{10}^2}}&{1.329{\rm{ }}65 \times {{10}^3}}&{ - 6.414{\rm{ }}28 \times {{10}^5}}&{5.228{\rm{ }}89 \times {{10}^3}}&{1.292{\rm{ }}13 \times {{10}^4}}&{ - 2.152{\rm{ }}21 \times {{10}^{ - 1}}}\\ {3.627{\rm{ }}36 \times {{10}^5}}&{1.706{\rm{ }}88 \times {{10}^5}}&{ - 2.008{\rm{ }}89 \times {{10}^3}}&{1.801{\rm{ }}09 \times 10}&{5.523{\rm{ }}76}&{3.954{\rm{ }}19 \times {{10}^3}} \end{array}} \right]。\end{array} $ |
同时，需要检验每一个回归系数的显著性，应用式(7)对其进行t检验计算，刚度矩阵对应的t检验矩阵T为:
| $ \begin{array}{l} \mathit{\boldsymbol{T}} = \\ \left[ {\begin{array}{*{20}{c}} {551.777{\rm{ }}6}&{0.891{\rm{ }}4}&{ - 0.421{\rm{ }}9}&{0.3231}&{0.2730}&{663.173}\\ { - 3.888{\rm{ }}1}&{1{\rm{ }}575.420{\rm{ }}4}&{ - 1.107{\rm{ }}8}&{1.516{\rm{ }}5}&{1.713{\rm{ }}8}&{775.382{\rm{ }}1}\\ {1.045{\rm{ }}3}&{ - 0.214{\rm{ }}5}&{371.457{\rm{ }}1}&{ - 371.759{\rm{ }}9}&{ - 370.682{\rm{ }}4}&{0.971{\rm{ }}8}\\ { - 0.110{\rm{ }}2}&{0.406{\rm{ }}9}&{ - 335.113}&{443.663{\rm{ }}2}&{370.649{\rm{ }}9}&{ - 0.778{\rm{ }}4}\\ { - 0.597{\rm{ }}5}&{1.208{\rm{ }}0}&{ - 164.726{\rm{ }}0}&{158.519{\rm{ }}3}&{1{\rm{ }}199.428{\rm{ }}2}&{ - 0.048{\rm{ }}5}\\ {240.816{\rm{ }}9}&{135.048{\rm{ }}1}&{ - 0.445{\rm{ }}5}&{0.471{\rm{ }}5}&{0.442{\rm{ }}8}&{771.913{\rm{ }}0} \end{array}} \right]。\end{array} $ |
取显著性置信度α=0.05和n=24，根据t的分布临界值表可知其对应的临界值t0=1.710 9。寻找T矩阵中tij值小于t0的对应位置并记录全部的(i，j)数对，然后将刚度矩阵K中相应的Kij值置零，并进行二次线性回归，得到新刚度矩阵为:
| $ \begin{array}{l} \mathit{\boldsymbol{K}} = \\ \left[ {\begin{array}{*{20}{c}} {5.309{\rm{ }}59 \times {{10}^8}}&0&0&0&0&{2.168{\rm{ }}17 \times {{10}^5}}\\ { - 5.801{\rm{ }}015 \times {{10}^5}}&{4.442{\rm{ }}67 \times {{10}^7}}&0&0&{4.100{\rm{ }}05 \times {{10}^2}}&{8.680{\rm{ }}89 \times {{10}^4}}\\ 0&0&{1.732{\rm{ }}34 \times {{10}^8}}&{ - 1.468{\rm{ }}79 \times {{10}^6}}&{ - 4.785{\rm{ }}68 \times {{10}^5}}&0\\ 0&0&{ - 2.652{\rm{ }}76 \times {{10}^6}}&{2.979{\rm{ }}80 \times {{10}^4}}&{8.126{\rm{ }}20 \times {{10}^3}}&0\\ 0&0&{ - 6.359{\rm{ }}86 \times {{10}^5}}&{8.183{\rm{ }}10 \times {{10}^3}}&{1.290{\rm{ }}71 \times {{10}^4}}&0\\ {3.622{\rm{ }}74 \times {{10}^5}}&{1.716{\rm{ }}40 \times {{10}^5}}&0&0&0&{3.954{\rm{ }}39 \times {{10}^3}} \end{array}} \right]。\end{array} $ |
简化后的线性方程需要进一步进行F检验，其中显著性置信度α取0.05，则可以得到各行的F临界值，如表 2所示。
| 表2 F检验值及其临界值 Tab. 2 F value and its critical value |
 |
由表 2可知，K1至K6行的F检验值都远大于其相应的临界值，故其6行均可通过检验，所有方程通过了F检验，满足线性回归的显著性要求。
根据式(9)可进一步得到6个线性方程的复相关系数，如表 3所示。
| 表3 复相关系数R2值 Tab. 3 Multiple correlation coefficient R2 |
 |
由表 3可知每一行的复相关系数R2均接近1，表明其拟合优度很高，同时验证了线性回归方程具有很好的回归效果。
为了进一步验证上述方法求解得到刚度矩阵的准确性，任意选取2组外力，将同一外力下的刚度计算值和ANSYS分析值变形量进行对比，结果如表 4所示。
| 表4 ANSYS值和计算值 Tab. 4 Results of ANSYS analysis and theoretical calculation |
 |
根据上述方法求得的刚度矩阵，建立了载荷与变形的解析表达式，从而能够方便快速的求解这类不规则的复杂构件在任意外部载荷下的6维变形量。对比分析表明，计算结果与ANSYS分析结果相差很小，最大相对误差仅为0.11%，完全满足工程需求。
4 结论1) 基于多元线性回归理论提出求解不规则构件刚度矩阵的统一解法。首先采用ANSYS分析得到力和位移值的多组数对，利用多元线性回归分析求解回归系数，然后应用t检验对方程形式进行简化，并对简化后的方程进行F检验和拟合优度检验，验证方程的回归效果，从而解决了传统方法无法求解不规则构件刚度矩阵的难题。
2) 应用该方法求解了U型轴复杂构件的刚度矩阵，在相同外力作用下，对构件变形量的理论计算值和ANSYS分析值进行了对比，最大误差仅为0.11%，验证了方法的准确性。
3) 该方法为后续对整机做以刚度为目标的优化设计提供了重要依据。可将复杂构件的刚度矩阵作为常量与规则构件的刚度矩阵统一建模得到整机的刚度矩阵解析模型，解决了多数刚度优化设计中需要对不规则构件进行简化处理的问题，能够有效提高分析精度。
| [1] |
Ai Qinglin, Huang Weifeng, Zhang Hongtao, et al. Review of stiffness and statics analysis of parallel robot[J]. Advances in Mechanics, 2012, 42(5): 583-592. [艾青林, 黄伟锋, 张洪涛, 等. 并联机器人刚度与静力学研究现状与进展[J]. 力学进展, 2012, 42(5): 583-592. DOI:10.6052/1000-0992-11-073] |
| [2] |
Tannous M, Caro S, Goldsztejn A. Sensitivity analysis of parallel manipulators using an interval linearization method[J]. Mechanism and Machine Theory, 2014, 71: 93-114. DOI:10.1016/j.mechmachtheory.2013.09.004 |
| [3] |
Gosselin C. Stiffness mapping for parallel manipulators[J]. IEEE Transactions on Robotics and Automation, 1990, 6(3): 377-382. DOI:10.1109/70.56657 |
| [4] |
Chen Xiaogang, Sun Yu, Peng Binbin, et al. Distribution characteristics of static stiffness for 6-UPS parallel machine tool[J]. Journal of Nanjing University of Science and Technology, 2013, 37(6): 926-933. [陈小岗, 孙宇, 彭斌彬, 等. 6-UPS并联机床静刚度分布特性[J]. 南京理工大学学报, 2013, 37(6): 926-933.] |
| [5] |
Yu Jingjun, Bi Shusheng, Zong Guanghua. Stiffness matrix method for displacement analysis of fully spatial compliant mechanisms[J]. Journal of Beijing University of Aeronautics and Astronautics, 2002, 28(3): 323-326. [靖军, 毕树生, 宗光华. 空间全柔性机构位置分析的刚度矩阵法[J]. 北京航空航天大学学报, 2002, 28(3): 323-326.] |
| [6] |
Yu Jingjun, Bi Shusheng, Zong Guanghua, et al. Analysis for the static stiffness of a 3DOF parallel compliant micromanipulator[J]. Chinese Journal of Mechanical Engineering, 2002, 38(4): 7-10. [于靖军, 毕树生, 宗光华, 等. 3自由度柔性微机器人的静刚度分析[J]. 机械工程学报, 2002, 38(4): 7-10.] |
| [7] |
Fu Tie, Ding Hongsheng, Pang Siqin, et al. Simulation of the machining error of variable-axis NC machine tool based on static stiffness[J]. Journal of Beijing Institute of Technology, 2002, 22(6): 672-674. [付铁, 丁洪生, 庞思勤, 等. 基于静刚度的变轴数控机床加工误差仿真研究[J]. 北京理工大学学报, 2002, 22(6): 672-674.] |
| [8] |
Li Yuwen, Zhang Hua, Yang Jianxin, et al. Finite element analysis and experimental study for the stiffness of a 6-UPS parallel kinematic machine[J]. China Mechanical Engineering, 2004, 15(2): 112-115. [李育文, 张华, 杨建新, 等. 6-UPS并联机床静刚度的有限元分析和实验研究[J]. 中国机械工程, 2004, 15(2): 112-115.] |
| [9] |
Jia Xiaoyong, Xu Chuansheng, Bai Xin. The invention and way of thinking on least squares[J]. Journal of Northwest University (Natural Science Edition), 2006, 36(3): 507-511. [贾小勇, 徐传胜, 白欣. 最小二乘法的创立及其思想[J]. 西北大学学报(自然科学版), 2006, 36(3): 507-511.] |
| [10] |
Fu Huimin, Wu Qiong. Multivariate regression analysis method for linear process[J]. Journal of Mechanical Strength, 2011, 33(3): 343-347. [傅惠民, 吴琼. 多元线性过程回归分析[J]. 机械强度, 2011, 33(3): 343-347.] |
| [11] |
Li Mintong, Zhao Jizheng, Yang Qing, et al. Study on tractor vibration based on wavelet and multiple linear regression[J]. Journal of Agricultural Mechanization Research, 2013(3): 40-45. [李敏通, 赵继政, 杨青, 等. 基于小波和多元线性回归的拖拉机振动分析[J]. 农机化研究, 2013(3): 40-45.] |