一種混凝土泵車臂架多姿態固有頻率的數值算法

張趙威，吳運新,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.