大型重载桁架机器人主要是指桁架式直角坐标机器人，是专门针对各种数控机床与CNC加工中心配置、高效代替人工实现自动上下料的自动化装备，广泛应用于金属加工、零部件装配、仓储搬运等行业^[1]. 大型重载机器人由于负载大于几百公斤，行程较长，如何保证机器人在高速度、高负载、高强度工作条件下的工作性能将成为制约其应用的关键因素，因此对重载桁架机器人进行静动态特性分析十分重要.

目前，国内外学者对结构进行静动态特性分析与优化的研究颇多. Zhu等^[2-5]基于有限元方法，分析桁架式机器人的静态或动态特性. Wu等^[6-12]基于结构的静态或动态特性分析，构建不同结构的优化模型，提高结构的性能. 彭军^[13]通过建立门式起重机主梁有限元模型，对主梁结构进行静动态特性分析. 彭真^[14]采用有限元分析方法，对四自由度机器人进行静动态特性分析. Li等^[15-16]对桁架结构进行拓扑优化. Sorokin等^[17]通过选择结构的设计参数，开展梁单元结构的优化设计，与初始设计形状进行动力学分析对比. 李世杰等^[18]以焊接机器人的薄弱环节为对象，开展模态分析与谐响应分析. Agarwal 等^[19]基于遗传算法和粒子群优化算法，以质量最小化和竖直位移最小化为优化目标，对桁架结构进行优化，对比了2种算法的优化效果.

根据文献[7, 8, 11]的阐述可知，目前国内外对结构的静动态特性进行综合分析以及优化的研究多以机床结构为对象，对桁架类机器人进行综合静动态特性分析及结构优化的研究较少，缺少相应的分析对比，尤其是重载桁架机器人方面，缺少基于静动态特性分析进行综合评估以及结构性能改善方面的研究.

为了探究重载桁架机器人横梁结构在承受极限负载压力下的静动态特性并改善该结构的静力学性能，本文利用有限元方法对机器人横梁结构进行静动态特性分析，将仿真与理论结果进行对比；基于不同实验设计方法对结构进行尺寸优化，综合评估结构的静动态特性与优化设计的可行性，为设计提供参考. 与传统的结构有限元分析与优化方法相比，本文更注重结构的整体性能与优化方法的实用性，综合考虑结构的静动态特性，开展结构的优化改进；通过对比分析，对结构进行更全面的评估分析并指导设计改进. 研究结果为横梁结构的进一步设计与动力学分析提供仿真与理论参考.

重载桁架机器人的主要结构是由立柱、支撑梁、横梁、竖梁及连接部分组成. 如图1(a)所示为重载桁架机器人的三维模型. 横梁结构由于采用两端支撑，中间没有辅助支撑结构且承担着负载与竖梁组合结构的重力，因此横梁结构在机器人运转时更容易发生较大变形与共振现象. 为了研究横梁结构的静动态特性，通过SolidWorks建模软件，建立横梁结构的简化三维模型，如图1(b)所示. 模型的基本参数如表1所示. 表中，S为截面外轮廓尺寸，L为长度，T为厚度，E为弹性模量，ρ为材料密度。总截面惯性矩为4.96×10⁸ mm⁴.

基于ANSYS Workbench有限元分析软件，对横梁结构进行静力学分析. 为了探究横梁在以实体单元或壳单元进行划分时仿真结果的差异，仿真时网格单元分别选择实体单元与壳单元^[20]对横梁结构进行划分. 结构不同实体之间采用绑定接触方式. 横梁原始截面与简化截面以及网格划分结果如图2所示.

网格划分完成后，需要对模型进行边界条件与极限压力的加载. 其中横梁结构两端面采用固定约束. 压力是通过预先在滑轨表面添加的印记面来实现压力分别在结构3个位置上的加载，每个加载位置存在2个加载面，加载面的距离为768 mm. 压力加载的简化模型如图3的位置A、B、C所示. 压力通过对竖梁、负载及连接部件进行测量得到. 其中负载质量为1 000 kg，竖梁与连接部件质量为780 kg.结构所受的总负载压力约为17 800 N. 作用在结构每个加载面的压力为8 900 N.

如图4所示为横梁结构的等效变形云图. 图中，D_E、D_S分别为基于实体单元和壳体单元划分的横梁结构变形量. 可以看出，结构在分别以实体单元划分与壳单元进行划分时，云图显示具有很好的一致性，且差值相差很小. 应力分析时，由于模型中存在尖角造成的应力集中，在应力分析对比时选择与压力直接接触的滑轨结构作为应力对比对象，结果对比如表2所示. 表中，D_A、D_B、D_C分别为位置A、B、C的最大变形量，σ_A、σ_B、σ_C分别为位置A、B、C的最大应力，N_d为单元个数，t_a为平均计算时间. 对比结果表明，虽然横梁结构实体单元模型与壳单元模型在应力与变形量方面相差很小，但是在横梁以壳单元进行划分时，单元数量可以缩减约60%，平均计算速度提高约29%. 静力学分析时，采用壳单元代替实体单元对结构进行划分，能够提高建模与仿真效率.

由表2的仿真分析结果可知，结构在未加入重力作用时Y轴的最大变形量为0.470 mm. 为了探究结构在加入重力因素后的综合变形量，在上述仿真的基础上，对横梁结构添加重力因素进行分析. 表3中，S_b为结构简化前变形量，S_a为结构简化后变形量，D_F为差值. 从表3可知，结构未简化前在位置C处最大总变形量达到0.586 mm. 为了进一步提高横梁定位精度，减小结构的总变形量，对结构进行优化. 由于结构优化分析需要多次模型重构，为了方便模型建立与网格划分，对横梁模型进行适当简化处理，简化前后的横梁截面与网格划分如图5所示. 由表3截面简化前、后的仿真结果对比可以看出，简化后截面的仿真结果与简化前的仿真结果差值较小. 通过仿真发现，简化后截面除去齿条结构后在图3位置C处受力时的变形量为0.575 mm，因此齿条结构对横梁结构变形量的影响不大. 结构优化时，为了提高优化效率，可以对横梁进行适当简化处理，以图3 中的C处力加载位置进行优化求解.

由于有限元模型结果的精度受模型简化程度、网格单元质量、边界条件的加载、有限元算法等因素的影响，为了验证仿真模型的准确性，对横梁模型进行理论计算对比. 根据图3中横梁结构的力加载位置，建立两端固定梁理论力学模型，如图6所示.

如图6(c)所示为两端固定、中间对称加载的力学模型. 根据理论力学知识可知，该力学模型可以分解为分别受弯矩以及单个力的梁模型^[21]，如图7所示.

由图7的梁等效力学模型与理论力学公式可知，

(1) $\left[\! {\begin{array}{*{20}{c}} {{\theta _1}}&{{\theta _2}} \\ {{\theta _3}}&{{\theta _4}} \end{array}} \!\right] = \left[ \!{\begin{array}{*{20}{c}} {{M_{\rm{L}}}l/(3EI)}&{{M_{\rm{R}}}l/(6EI)} \\ {Fab(l + b)/(6EIl)}&{Fab(l + a)/(6EIl)} \end{array}} \!\right],$

(2) $\left[\!\!\! {\begin{array}{*{20}{c}} {{\omega _1}}\!\!\!&\!\!\!{{\omega _2}} \\ {{\omega _3}}\!\!\!&\!\!\!{{\omega _4}} \end{array}} \!\!\!\right] = \left[\!\! \!{\begin{array}{*{20}{c}} {{M_{\rm{L}}}{l^2}/(16EI)}\!\!\!&\!\!\!{{M_{\rm{R}}}l^2/(16EI)} \\ {Fb(3{l^2}{\rm{ - 4}}{b^2})/(48EI)}\!\!\!&\!\!\!{Fa(3{l^2}{\rm{ - }}4{a^2})/(48EI)} \end{array}} \!\!\!\right].$

式中：M_L、M_R为左侧和右侧的弯矩，E为材料的弹性模量，I为截面惯性矩，l为梁的长度，a、b分别为力F与梁两端之间的距离，θ₁、θ₂、θ₃、θ₄为转角，ω₁、ω₂、ω₃、ω₄为挠度.

由于在横梁结构加载的2个力F相等，可知M_L=M_R. 又由 ${\theta _1} + {\theta _2} - {\theta _3} - {\theta _4} = 0$，可得

(3) $ \frac{{M}_{{\rm{L}}}l}{3EI}+\frac{{M}_{{\rm{R}}}l}{6EI}-\frac{Fab\left(l+b\right)}{6EIl}-\frac{Fab\left(l+a\right)}{6EIl}=0 .$

计算可得，M_L=M_R=16062.725 N·m.

将M_L=M_R=16062.725 N·m，F=8 900 N，E=2 × 10¹¹ Pa，I=4.96 × 10⁸ mm⁴代入，可得C处挠度：

${\omega _1} + {\omega _2} - {\omega _3} - {\omega _4} = {0.353^{}}\;{\rm{mm}}.$

同理可得，力在位置A、B加载时梁的挠度分别为0.074、0.281 mm.

由上述仿真结果与理论计算结果对比可以发现，理论计算结果与仿真结果有较大的偏差. 进一步分析发现，造成上述偏差较大的主要原因是在仿真时最大变形发生在局部，理论模型为模型的整个截面，因此误差较大. 基于上述原因，提取仿真模型截面在Y方向上变形最大位置处的数据，将数据导出至Excel软件，计算截面变形量的平均值并与理论值进行比较. 如图8所示为横梁结构在3种压力加载位置下的截面变形云图. 图中，D_b为梁挠度最大处截面的变形量。数据处理后得到横梁组合模型在A、B、C 3个位置力加载下梁截面的变形量绝对值的平均值分别为0.084、0.284、0.349 mm. 如表4所示为理论与仿真数据的对比. 表中，D_th为理论值，D_s为仿真值，D_er为误差.

由表4可知，横梁结构的理论计算变形量与仿真最大变形处截面的平均变形量相差很小，最小误差仅为0.003 mm，因此仿真模型具有较高的准确性.

基于静力学划分的有限元模型，对横梁结构进行模态分析. 边界约束条件与静力学分析相同. 如表5所示为结构前10阶固有频率. 表中，f为频率. 结构1~6阶振型云图如图9所示. 图中，D_m为横梁变形量。

由表5可知，横梁结构的前10阶固有频率分布在47~178 Hz，固有频率从第2阶至第3阶发生了突变. 结构前2阶固有频率都小于50 Hz. 由前6阶结构的振型云图可知，与其他振型相比，第3、6阶振型下的变形量较大，分别为1.789 6和1.822 4 mm.

在对该类型机器人进行实验及动力学仿真时发现，由于模型间存在装配间隙、表面不光滑等原因，机器人在运行过程中存在一定的上下振动. 为了考察横梁结构在外界激励载荷下的动态响应，基于结构模态分析结果，针对如图3中C处力加载位置时的横梁组合模型，采集结构的水平面与竖直面. 如图10所示，对该组合结构进行谐响应分析，研究动态响应特性. 由于谐响应分析是研究结构在受简谐载荷作用时的稳态响应^[22]，结构的多自由度运动微分方程为

(4) ${{M}}\ddot {{x}} + {{C}}\dot {{x}} + {{Kx}} = {{F}}.$

由于式(4)中F和x为简谐矩阵，设频率为ω，经过欧拉转化后，联立可得谐响应分析的运动方程：

(5) $( - {\omega ^2}{{M}} + {\rm{i}} \omega {{C}} + {{K}})({{{x}}_1} + {\rm{i}} {{{x}}_2}) = {{{F}}_1} + {\rm{i}} {{{F}}_2}.$

式中：M为质量矩阵，C为阻尼矩阵，K为刚度矩阵，F₁、F₂、 ${x_1}$、 ${x_2}$分别为F与 $x$的实部和虚部.

对式(5)进行求解，可得系统谐响应的位移变化量. 谐响应分析结果如图11所示. 图中，A为振幅.

如图11所示为结构谐响应位移响应分析结果. 可以看出，当频率为52.5 Hz时，结构的水平面与竖直面在Y方向上振幅出现峰值，峰值分别为1.4533、1.446 8 mm. 此时频率接近结构的第2阶固有频率.

从上述静动态分析结果可知，横梁结构在极限负载与重力的双重作用下在位置C处的最大总变形量为0.586 mm. 谐响应分析显示，结构第2阶固有频率附近对结构竖直方向上的变形量影响较大.为了减小横梁结构在极限负载压力下的总变形量，同时实现结构的轻量化并维持结构较高的一阶固有频率. 以横梁以及上下滑轨组成的横梁结构作为优化对象. 利用SolidWorks三维建模软件提取横梁结构10个设计参数进行优化求解. 横梁结构简化后的截面与优化尺寸参数及网格划分截面如图12所示. 其中横梁结构的尺寸信息包括：横梁截面的尺寸信息(宽度/高度/相邻面的厚度，共4个尺寸参数)、滑轨截面尺寸信息(宽度(上、下各1个，共2个)/高度(上、下各1个，共2个)/定位尺寸(上、下各1个，共2个),共6个尺寸，非标准化)、齿条结构(已标准化，不对尺寸参数进行更改). 由于横梁结构设计时要求上、下滑轨尺寸信息保持一致，需要对上、下滑轨的尺寸参数进行相关性设置，使优化时尺寸参数同步变化. 截面尺寸参数具体设置为：P₈=P₅，P₉=P₆，P₁₀=P₇，因此只需对7个尺寸参数进行变化范围的设置. 如表6所示为各个尺寸参数的变化范围.

从文中结构的静力学分析可知，结构的等效应力约为22 MPa，远小于材料的许用应力. 对结构进行优化时，不将结构的应力作为优化目标. 根据选用电机的参数可知，伺服电机的额定转速为1 500 r/min，频率为0~25 Hz. 为了避开电机的激振频率，结合结构的静动态特性分析. 综合考虑结构的质量、最大总变形量及一阶固有频率等性能指标，建立结构的优化模型如下：

(6) $\left. \begin{array}{l} \min\; \left\{ {F(x)} \right\};\\ {\rm{s}}{\rm{.t}}.\;\;\;\;f(x) \geqslant 35\;{\rm{Hz,}}\\ \;\;\;\;\;\;\;\;m(x) \leqslant {m_0},\\ \;\;\;\;\;\;\;\;x \in ({x_{\rm{L}}},{x_{\rm{H}}}). \end{array} \right\}$

式中：F(x)为结构的最大总变形量，f(x)为结构的一阶固有频率，m(x)为结构的总质量，m₀为结构的原始质量， $x$为结构的设计变量， ${x_{\rm{L}}}$为设计变量的下限， ${x_{\rm{H}}}$为设计变量的上限.

为了创建准确的响应面模型，在设计变量可行域内取得一定数量的样本点进行拟合. 分别采用中心复合设计(central composite design，CCD)与最佳填充空间设计(optimal space-filling design，OSF)^[23]2种方式，生成所需的样本点. CCD实验设计是Workbench默认的DOE类型. 它提供了一个筛选集来确定元模型的总体趋势，更专注设计区域边界附近的参数、具有序贯、高效、灵活的特点，但最多支持20个输入参数. OSF实验设计是根据某些指定标准，创建最佳空间填充实验设计计划，能够在整个设计空间中平均分配设计参数，从而以最少的数量获得对设计的最大了解. OSF设计缺点是未必涵盖极端情况（即设计空间的角落），且选择的设计点太少时会导致响应预测的质量降低. 为了对比2种设计生成样本点的最终优化效果，样本点数量设置为79组。P₁设计变量在2种实验设计方式下生成样本点的分布情况如图13(a)所示，由2种实验设计方式生成的质量点分布如图13(b)所示. 图中，m为质量，N为生成点数量.

响应面法(response surface method，RSM)是利用函数关系拟合仿真模型，采用筛选试验来确定优化方向. 以RSM为指向，按照满足约束条件和质量最轻方向寻找最优解，设定好收敛条件，可以搜索到最佳设计^[24]. 非参数回归模型虽然在计算速度上较慢，但具有较高的精度，且本文的样本点数量较少，因此选择非参数回归模型^[25]模拟复杂的响应关系. 如表7所示为响应面模型的拟合精度. 表中，d为变形量，σ为应力. 评判标准中测定系数最佳值为1，均方根误差最佳值为0，最大相对残差最佳值、相对均方根误差最佳值、相对最大绝对误差、相对平均绝对误差最佳值都为0. 由表7可见，响应面模型满足精度要求.

优化设计中所涉及到的3个优化目标不可能同时达到最优解，需要对获得的Pareto解集进行求解得到妥协解，从妥协解中选取合理的优化方案. 本文基于多目标遗传算法(multi-objective genetic adventure，MOGA)，求解Pareto解集^[26]，对结果进行优选，得到优化方案. 如图14所示为横梁结构的优化流程. 优化求解设置相关的参数如下: 初始样本数为100，每次迭代的样本数为100，最大允许Pareto百分比为70%，收敛稳定性百分比为2%，最大迭代次数为20，突变概率为0.01，交叉概率为0.98.

如图15所示为CCD和OSF实验设计的Pareto解集分布图。可以看出，在以变形量为主要优化目标的条件下，底部的点集具有更好优化效果.由解集点分布情况可知，OSF实验设计与CCD实验设计相比所得到的点集更分散，点集数值分布范围更广，因此OSF实验设计能够得到更好的优化目标值，优化效果更佳.

为了进一步对比2种实验设计的优化效果，分别提取2种实验设计方案下各5组设计推荐方案，进行量化比较，如表8~10所示。如图16所示为2种实验设计方式的优化效果图. 图中，n为优化方案编号. 如图17所示为CCD设计与OSF设计优化结果的对比曲线图. 其中图17(a)曲线上数据为最大总变形量，图17(b)曲线上数据为质量. 如表8、9所示为2种设计优化结果的具体尺寸. 如表10所示为这2种设计优化结果的频率、变形量、质量。表中，v_f、v_d、v_m分别为它们对应的变化量百分比，即为(优化后数值−原始数值)/原始数值.

如图16所示，2种实验设计方式下的优化结果与原模型相比，质量、最大总变形量、频率都降低. 由图17可知，CCD设计与OSF设计优化结果相比，OSF设计优化结果在质量与最大总变形量方面优化效果更好，在频率方面与CCD结果相比较小. 从表8~10所示的结构优化后尺寸参数及变化百分比可知，CCD实验设计的5个备选优化方案尺寸变化范围较小，优化结果相差较小；OSF实验设计优化结果变化明显，优化效果更好. 从OSF优化结果可以看出，结构经过优化后，从最大变形量方面考虑，优化后结构的最大总变形量降低了19.13%，质量降低了0.14%，频率降低了12.74%. 从轻量化方面考虑，优化后结构质量最大可以降低1.47%，最大总变形量降低了18.41%，频率降低了13.48%. 优化后结构的最小一阶固有频率为41.13 Hz，远离电机激振频率. 综合考虑质量、频率、最大总变形量等因素，OSF实验设计优化结果中的方案3优化效果最好.

（1）静力学分析时，横梁结构采用壳单元划分与结构采用实体单元进行划分时的仿真结果相差较小. 壳单元能够降低建模与单元划分时间，单元数量少，计算效率高. 壳单元更适合横梁结构静力学分析.

（2）根据横梁模型建立的理论力学模型计算可知，理论计算结果与仿真最大变形截面处的平均变形量有很好的一致性，最小误差为0.003 mm.

（3）横梁结构的前6阶固有频率分布在47~134 Hz，对结构进行谐响应分析发现，引起结构出现共振响应的频率为结构前3阶固有频率. 其中位移响应中第2阶固有频率的影响最大.

（4）OSF实验设计与CCD实验设计相比，得到的Pareto解集优化目标的分布范围更大，更容易得到良好的优化效果. 在OSF实验设计下，通过结构尺寸优化分析可知，结构质量降低了1.47%，总变形量降低了18.41%，频率降低了13.48%.

（5）从优化结果尺寸变化可以看出，结构尺寸参数P₃、P₄所代表的尺寸对结构的性能影响更大，可以将P₃、P₄代表的尺寸参数作为设计变量，进一步扩大尺寸搜寻范围，对结构进一步优化.