一种混凝土泵车臂架多姿态固有频率的数值算法

张赵威，吴运新,2，任 武

(1.中南大学 机电工程学院，长沙 410083；2.高性能复杂制造国家重点实验室，长沙 410083)

混凝土泵车作为一种重要的建筑机械，近些年来得到越来越广泛的应用。在工程应用中，泵车臂架的振动性能会直接影响施工质量、人员安全和设备寿命。研究者们在臂架的振动性能上做了大量的工作，代表的有郭立新等[1-2]，其工作主要是基于有限元软件研究臂架特定姿态下的振动相关性能及优化；而后吴运新[3]建立基于MSC. Partran/Nastran 的变姿态臂架振动分析系统实现变姿态的振动分析，而这种基于有限元软件方法计算臂架的振动特性时工作量过大、效率低，在臂架动力学分析以及车载监测[4]中应用存在困难。

传递矩阵法最初应用在计算串联转子的振动特性上，发挥出其建模和计算优势，后经不断发展形成完善的理论体系。由于其计算量小、建模方便、无需系统总体动力学方程、效率高等优点，被多个领域的学者们应用在许多重要机械系统工程问题上[5-9]。

为能够快速、高效的计算臂架任意姿态固有频率，本文提出一种数值方法。首先基于传递矩阵法计算臂架力学模型，以快速计算臂架任意固定姿态的固有频率；而后，通过多元函数拟合[10-11]臂架姿态与固有频率，得到固有频率关于姿态的数值计算公式[12]。并对文中力学模型的准确性，以及数值方法的计算的效果进行相关研究。

1 力学模型

混凝土泵车臂架结构主要包含长臂、连杆以及油缸，忽略其较小的几何特征得到结构简图如图1所示。臂架结构较为细长，而且由多个结构相似六连杆机构部分串联而成，因此将臂架划分成多个子结构能够有效的提高建模效率。

1. 固定平台; 2.液压缸Ⅰ; 3.臂杆Ⅰ; 4.液压缸Ⅱ; 5.连杆 Ⅰ-Ⅱ; 6.臂杆 Ⅱ; 7. 液压缸 Ⅲ; 8.连杆 Ⅱ-Ⅲ; 9.臂杆 Ⅲ; 10. 液压缸 Ⅳ; 11.连杆 Ⅲ-Ⅳ; 12.臂杆 Ⅳ

模型的建立流程：首先将臂架系统分解为两种子结构，而后进一步化分为体元件和铰元件，用矩阵表示各元件的力学特性，对这些元件传递矩阵按系统的联接关系进行组合与拼装而获得子结构传递方程和传递矩阵，而后进一步拼装得到系统总传递方程和总传递矩阵。最后利用边界条件，得到该系统的动力学特性。模型中将长臂部分作为弹性梁元件，连杆、油缸及活塞杆作为刚体元件，液压缸的油液等效为等刚度的弹簧铰元件计算[13]。

文中臂架为平面力学模型，定义各联接点状态矢量均为：

Zi,j=[XYΘzMzQxQy]T

在模态坐标系下，状态矢量中X,Y为线位移；Θz为角位移；Mz为内力矩；Qx,Qy为内力。

1.1 子结构的传递矩阵

根据臂架的结构特征将臂架划分为两种子结构A和子结构B，其简图如图2、3所示。

图2 子结构A简图

子结构A中定义0为子结构输入端，6为子结构输出端，1、2、4为刚性体元件，3为平面弹簧铰，5为等截面欧拉梁元件，则综合元件的传递方程可得式(1), 其中各元件的传递矩阵Ui见文献[5]。

(1)

由交叉点处的位移和受力关系可得：

(2)

由式(3)和式(2)可得：

(3)

式中：

E3=[0 0 0 1 0 0]

(4)

因此由式(3)可得传递方程：

UAZA=0

(5)

式(5)中子结构A的传递矩阵UA以及状态矢量ZA：

图3 子结构B简图

子结构B中0为子结构输入端，10为子结构输出端，2、3、5、7、8、9元件为刚性元件，6元件为平面弹簧铰，1、4元件为等截面欧拉梁元件。同子结构A的推导方式，根据各元件之间的联接关系可以得到传递方程如式(2)，即可以得到子结构B的传递矩阵UB以及状态矢量ZB。

(6)

式中

U1,3=U4U3;U1,5=U4E1U9

U2,1=E4U2U1;U2,4=E4U2E1U8

U3,2=-E5U7U6U5

U5,2=-E2U7U6U5

U6,2=E3U7U6U5;U7,1=E2U1

U7,4=-E2U8;U8,3=E2U3

U8,5=-E2U9;U9,4=E3U8

U10,5=E3U9;U11,1=E3U2U1

U11,4=E3U2E1U8

1.2 总传递矩阵

臂架是由一个子结构A和N-1个子结构B(本文模型N=4)串联拼接组成，因此子结构A的Z5,6与子结构B的Z0, 1,Z0, 5形成联接点，子结构B的Z9, 10与其后同样的子结构B′的Z0, 1,Z0, 5形成新的联接点，根据联接点的位移和受力关系可以得到臂架总传递方程如式(3)，其中臂架系统总传递矩阵Uall和总特征矢量Zall。

UallZall=0

(7)

式中UA_B,UB_A,UB_B′,UB′_B矩阵如下所示：

UA_B=[O6×12I6]

UB_A=[I6E1I6O6×24]

UB_B’=[O6×30I6]

UB’_B=[I6×6E1I6O6×24]

由边界条件知，边界状态矢量如下：

子结构A处边界状态矢量：

非末端子结构B(编号记为i)状态矢量为：

末端子结构B(编号记为N)状态矢量为：

(8)

通过数值计算的方式求解，预置有解区间用二分法逼近计算满足式(8)容差范围(本文容差设为-1e-7～1e-7)的(值，即可得臂架某一固定姿态下任一阶次的固有频率。

为计算其他任意固定姿态的固有频率，需重新修改总传递矩阵。模型中传递矩阵的连接方式不变，仅元件的传递矩阵因元件的角度改变发生变化，其余均不变。各个元件的旋转角度可由臂架的姿态角通过几何解析得到[14]，而后将传递矩阵以该转角进行坐标变换，如式(9)。计算即可得新姿态下的传递矩阵以及传递方程。从而得到臂架任意姿态的任意阶固有频率fk(θ1,θ2,θ3,θ4) (k=1, 2,…)，其中θi为各臂架相对水平面的角度。

(9)

1.3 模型仿真与验证

以某型号37 m四节臂泵车臂架为原型，分别使用Matlab编写的传递矩阵法程序和使用ANSYS软件计算的臂架模型进行计算，并对比分析前5阶固有频率。

其中有限元模型，在表中表示为fFEM。模型中的臂杆等承受弯曲载荷的元件均为Beam单元，其中刚性元件的材料弹性模量设置为所选材料弹性模量的100倍，以达到近似刚性的效果；油缸及各个连杆元件使用Link单元，其中油缸的刚度体现在单元材料的弹性模量上；光滑铰则用相邻节点释放转动自由度连接。

表1 力学模型相对有限元模型计算结果误差分析

对比计算结果可以发现，传递矩阵法的计算结果与ANSYS的计算结果非常接近，除个别较高频的误差超过0.5%，其余均处于较小的误差范围，尤其是较重视的一阶固有频率更为精准，从而证明了传递矩阵法在臂架动力学计算的准确性，由姿态变换仍能保证较高的计算精度，说明该模型的可靠性。并且该模型效率很高，使用普通微机计算一个固有频率值耗时约为0.03～0.05 s(精度为1e-7)，远优于有限元软件的计算效率，使得大规模计算成为可能。

2 拟合回归方程

臂架在运动过程中姿态变化幅度大，而其固有频率跟姿态有着很密切的对应关系，因而可以假想存在一个函数如式(10)，使得对于任意姿态都可快速计算其固有频率值。

f=Fk(θ1,θ2,θ3,θ4)

(10)

因而可使用数值拟合的方法获取一个工程上适用的固有频率计算公式，使得输入任意姿态角度即可得到相应的固有频率。

介于在臂架姿态变换中，是相邻臂架之间的夹角对固有频率起到直接作用，所以公式拟合时采用α1,α2,α3,α4作为直接变量，其中αi,θi之间的关系如式(11)，将f=Fk(θi)关于αi的函数方程写为式(12)：

α1=θ1;α2=θ2-θ1

α3=θ3-θ3;α4=θ4-θ3

(11)

f=Gk(α1,α2,α3,α4)

(12)

由于直接进行数据拟合的计算量太大，而且会有一些频率值随着αi的变化产生突变，直接拟合效果较差，因而将函数分成三部分进行推算如式(13)：

f=Gk(α1,α2,α3,α4)=fk,0+

(13)

式中：fk,0和αi均为0时的基姿态固有频率；Gk,i(αi)是单一角度变量下的增量函数；δk(α1,α2,α3,α4)是补偿函数。

2.1 计算Gk,i (αi)各项系数

以α1影响的一阶固有频率增量函数G1,1(α1)为例，保持α2,α3,α4等于零，α1在0～π/2范围内取n(n>50)个点，可以得到n个离散点，对这些点进行多项式拟合[15]并去除基项，即可得相对于基量的变化量：

(14)

图4 G1,1(α1)离散点与曲线

(15)

图5 G1,4 (α4)离散点与曲线(β1,4=1.134)

突变的情况一般为分为两段函数，式中βk,i为分段函数的分界点。经过数值计算发现此分界点的角度为臂架之间的死角，其中α3,α4影响较为明显，图5为G1,4(α4)的拟合曲线，在有突变的情况下得到的拟合曲线。将四种角度的增量函数相加，即可得到臂架任意姿态第k阶固有频率的增量项∑Fk,i(αi)。

2.2 计算δk(α1, α2, α3, α4)各项系数

增量函数的叠加所得公式会产生误差较大，所以在公式前两项推算之后添加一个补偿函数，意义为：

(16)

使αi在0～π/2中每项变量均取m个点，共计m4种姿态进行二次多项式回归，多项式的一般形式为：

(l=9,…,15)

(17)

通过这m4个点使用最小二乘拟合，即可得到臂架任意姿态第k阶固有频率的补偿项δk(α1,α2,α3,α4)。

由上述各项系数及式(11)，带入式(16)即可得到f=Fk(θ1,θ2,θ3,θ4)，臂架任意姿态第k阶固有频率的计算值。为评估拟合产生的偏差，对回归多项式与力学模型的计算结果任取10 000种姿态进行误差统计，结果如表2，其中η=|fk-Fk|/fk。

表2 经验公式相对力学模型计算结果误差统计

从表2中可以看出回归公式与力学模型的误差绝大多数都小于3%，一阶频率的误差更小，所以使用这种拟合方式可以保留力学模型的计算精度，同时再一次大幅提高求解固有频率的效率。

3 结 论

(1) 使用该方法得到的公式计算臂架任意姿态的固有频率，可以在保证计算结果准确性的基础上，极大地降低求解臂架固有频率的计算量，提高计算效率，使得车载固有频率计算成为可能；

(2) 该多元函数的混合拟合方法，可以解决函数拟合中的突变问题，同时保证了回归公式的计算精度。

[1]郭立新, 陈长征, 张国忠，等.混凝土泵车布料机构水平工况的动态分析[J]. 振动与冲击, 1999, 18(3): 79-82.

GUO Li-xin, CHEN Chang-zheng, ZHANG Guo-zhong,et al. The dynamic analysis of level operating mode of placing boom of truck concrete pump [J]. Journal of Vibration and Shock, 1999, 18(3): 79-82.

[2]宋建安, 董忠红, 吕彭民.水泥混凝土输送泵车整机模态[J]. 长安大学学报(自然科学版), 2004, 24(5): 104-106.

SONG Jian-an, DONG Zhong-hong, LÜ Peng-min. Modalling for concrete-auto-pump vehicle [J]. Journal of Chang’an University (Natural Science Edition), 2004, 24(5): 104-106.

[3]吴运新, 钟志宏, 滑广军. 基于 MSC.Patran/Nastran的泵车臂架分析系统的研究[J]. 郑州大学学报(工学版), 2010, 31(6): 83-86.

WU Yun-xin, ZHONG Zhi-hong, HUA Guang-jun. Study on anaysis system for boom of truck mounted concrete pump based on MSC.Patran/Nastran [J]. Journal of Zhengzhou University (Engineering scinece), 2010, 31(6): 83-86.

[4]黄毅, 王佳茜, 邝昊. 混凝土泵车便携式状态监测与故障诊断系统研究[J]. 中国工程机械学报, 2010, 8(4): 441-443.

HUANG Yi, WANG Jia-qian, KUANG Hao. Reaserach and development on condition monitoring and fault diagnosis system for boom pumps [J]. Chinese Journal of Construction Machinery, 2010, 8(4): 441-443.

[5]芮筱亭, 何斌, 王国平，等.多体系统传递矩阵法及其应用[M]. 北京: 科学出版社, 2008.

[6]芮筱亭, 戎保. 多体系统传递矩阵法研究进展[J]. 力学进展, 2012, 42(1):1-17.

RUI Xiao-ting, RONG Bao. Advances in transfer matrix method for multibody system dynamics [J]. Advances in Mechanics, 2012, 42(1): 1-17.

[7]刘杨, 刘挺, 冯霏，等.基于传递矩阵法的振动压实系统动力学分析[J]. 东北大学学报(自然科学版), 2010, 31(1): 96-98.

LIU Yang, LIU Ting, FENG Fei, et al. Dynamic analysis based on transfer matrix for vibration compaction system [J]. Journal of Northeastern University (Natural Science), 2010, 31(1): 96-98.

[8]贺军义, 芮筱亭, 王国平，等.多管火箭发射过程中定向器振动特性研究[J]. 振动与冲击, 2012, 31(1): 35-38.

HE Jun-yi, RUI Xiao-ting, WANG Guo-ping, et al. Vibration characteristic of launch guider of a MLRS in firing process [J]. Journal of Vibration and Shock, 2012, 31(1): 35-38.

[9]孙建鹏, 李青宁. 求解结构自振频率的精细传递矩阵法[J]. 世界地震工程, 2000, 25(2): 140-145.

SUN Jian-peng, LI Qing-ning. Precise transfer matrix method for resolving natural frequencies of structures [J]. World Earthquake Engineering, 2000, 25(2): 140-145.

[10]Bezemer E，Rutan S C. Multivariate curve resolution with non-linear fitting of kinetic profiles [J]. Chemometrics and Intelligent Laboratory Systems, 2001,59(1-2): 19-31.

[11]邱竹贤, 邱天爽, 王兆文，等. 铝电解实验数据的回归分析和经验公式拟合[J]. 东北大学学报, 2003, 24(4): 352-357.

QIU Zhu-xian, QIU Tian-shuang, WANG Zhao-wen, et al. Regression analysis of aluminum electroyte data and their empirical expressions [J]. Journal of Northeastern University (Natural Science), 2003, 24(4): 352-357.

[12]Chattha H T,Huang Y，Zhu X，et al. An empirical equation for predicting the resonant frequency of planar inverted-F antennas [J]. Antennas and Wireless Propagation Letters, IEEE, 2009,8: 856-860.

[13]戴云飞. 液压缸液压刚度的计算[J]. 有色金属设计, 1999, 26(1): 61-63.

DAI Yun-fei. The calculation of hydraulic cylinder stiffness[J]. Nonferrous Metals Design, 1999, 26(1): 61-63.

[14]胡凡. 泵车臂架低周应力软测量模型的建立及其在疲劳损伤计算中的应用[D]. 长沙: 中南大学, 2012.

[15]Chapra S C, Canale R P, Numerical methods for engineers [M]. New York: McGraw-Hill, 2010.