易燦明， 余海東， 王 皓,b

(上海交通大學(xué) a. 上海市復(fù)雜薄板結(jié)構(gòu)數(shù)字化制造重點(diǎn)實驗室； b. 機(jī)械系統(tǒng)與振動國家重點(diǎn)實驗室， 上海 200240)

軟體機(jī)器人由于其動作靈活、形態(tài)多變和適用性強(qiáng)等特征，已成為目前機(jī)器人發(fā)展的研究熱點(diǎn).軟體機(jī)器人的驅(qū)動主要采用復(fù)合結(jié)構(gòu)，通過將硬度較大的驅(qū)動材料和柔性介質(zhì)材料進(jìn)行耦合，形成軟硬結(jié)合的復(fù)合柔性結(jié)構(gòu).該結(jié)構(gòu)具有的大變形特征使得傳統(tǒng)的有限元建模方法難以對其運(yùn)動和變形進(jìn)行精確的描述[1]；同時考慮驅(qū)動器的存在，多單元耦合也給復(fù)合柔性結(jié)構(gòu)的建模帶來一定的難度.因此，建立精確的大柔性復(fù)合結(jié)構(gòu)動力學(xué)模型對上述結(jié)構(gòu)的運(yùn)動控制及動力學(xué)特性分析具有重要意義.

在多單元耦合的復(fù)合結(jié)構(gòu)建模方面，傳統(tǒng)的基于浮動坐標(biāo)或者轉(zhuǎn)角坐標(biāo)的有限元方法在建立多耦合單元時容易出現(xiàn)耦合條件不充分、在大變形情況下單元間出現(xiàn)穿透現(xiàn)象等問題[2].Carrera[3]總結(jié)了常用于多層板、殼結(jié)構(gòu)的建模方法，其中單一板模型(ESL)假定各層板單元以單層板的形式變形，模型不夠精確，且沒有考慮柔性結(jié)構(gòu)的大變形特征；在Layerwise模型中，通過在各層板之間插入高階邊界條件，提升了模型的精度，但同時也極大地增加了結(jié)構(gòu)的自由度數(shù)量.Zhao等[4]基于不同單元之間的變形協(xié)調(diào)條件，提出了單元間的坐標(biāo)轉(zhuǎn)換矩陣，不僅提高了耦合建模的精度，同時降低了耦合單元的自由度數(shù)量.在處理復(fù)合結(jié)構(gòu)大變形過程中剛體運(yùn)動與柔性變形的耦合問題方面，Shabana等[5]提出絕對節(jié)點(diǎn)坐標(biāo)法為復(fù)合結(jié)構(gòu)建立高階單元模型，并基于該方法建立了三維梁單元模型，進(jìn)而運(yùn)用連續(xù)介質(zhì)力學(xué)方法[6]對單元的剛度矩陣進(jìn)行了推導(dǎo).絕對節(jié)點(diǎn)坐標(biāo)法將單元節(jié)點(diǎn)坐標(biāo)定義在同一個全局坐標(biāo)系下，同時使用節(jié)點(diǎn)的斜率坐標(biāo)代替?zhèn)鹘y(tǒng)的轉(zhuǎn)角坐標(biāo)對節(jié)點(diǎn)的方向進(jìn)行描述，在解決結(jié)構(gòu)大變形、大轉(zhuǎn)動問題時具有更高的精度.隨后，Mikkola等[7]通過該方法為三維板單元和殼單元建立了非增量有限元模型；趙春璋等[8]基于絕對節(jié)點(diǎn)坐標(biāo)法建立了變截面梁單元模型，并針對截面特性對單元動力學(xué)性能的影響進(jìn)行了研究；李彬等[9]通過大變形柔性梁系統(tǒng)驗證了用絕對節(jié)點(diǎn)坐標(biāo)法建模的正確性和精確性.在壓電材料驅(qū)動器建模方面，Gilardi等[10]基于絕對節(jié)點(diǎn)坐標(biāo)法建立了由壓電驅(qū)動器驅(qū)動的大柔性梁模型，并研究其在高速旋轉(zhuǎn)及振動控制下的靜力學(xué)和動力學(xué)特性；Nada等[11]基于絕對節(jié)點(diǎn)坐標(biāo)法對多層含壓電薄片的薄板單元進(jìn)行了有限元建模，并對其運(yùn)動特性進(jìn)行了分析.目前關(guān)于壓電材料驅(qū)動器的研究大多通過傳統(tǒng)的層合理論對壓電驅(qū)動器和主體單元進(jìn)行耦合，單元耦合建模精度有待提高.

鑒于絕對節(jié)點(diǎn)坐標(biāo)法采用的高階單元模型可對具有大變形特征的柔性復(fù)合結(jié)構(gòu)運(yùn)動過程中的剛?cè)狁詈咸匦赃M(jìn)行精確描述，同時考慮到耦合單元精確建模的需求，本文基于絕對節(jié)點(diǎn)坐標(biāo)法，根據(jù)變形協(xié)調(diào)條件推導(dǎo)坐標(biāo)關(guān)聯(lián)矩陣，建立了大柔性梁板耦合單元模型，并對其動力學(xué)方程進(jìn)行了求解.此外，引入壓電材料本構(gòu)方程，將梁單元作為復(fù)合結(jié)構(gòu)的驅(qū)動器，用Matlab軟件對壓電驅(qū)動柔性板結(jié)構(gòu)進(jìn)行數(shù)值仿真分析，研究該結(jié)構(gòu)在外加驅(qū)動下的靜力學(xué)和動力學(xué)特性，并分析各項參數(shù)對其性能的影響，以期為軟體機(jī)器人的結(jié)構(gòu)及驅(qū)動設(shè)計提供參考.

1 基于絕對節(jié)點(diǎn)坐標(biāo)法的柔性單元建模

1.1 單元的變形描述

根據(jù)絕對節(jié)點(diǎn)坐標(biāo)法，柔性板單元的所有坐標(biāo)均定義在全局坐標(biāo)系下，三維四節(jié)點(diǎn)板單元模型如圖1所示.圖中：O-XYZ為全局坐標(biāo)系；o-xyz為與板單元固連的浮動坐標(biāo)系，描述板單元上任意一點(diǎn)的相對位置.

圖1 三維四節(jié)點(diǎn)板單元模型Fig.1 Three-dimensional four-node plate element

該板單元模型包含A、B、C和D4個節(jié)點(diǎn)，每個節(jié)點(diǎn)由3個位置坐標(biāo)和9個斜率坐標(biāo)進(jìn)行描述，板單元上節(jié)點(diǎn)i(i=A,B,C,D)的坐標(biāo)可定義為

(1)

式中：ri為節(jié)點(diǎn)i的全局位置坐標(biāo)；?ri/?k(k=x,y,z)為節(jié)點(diǎn)i處的浮動坐標(biāo)系相對于全局坐標(biāo)系在k軸上的斜率坐標(biāo).由于每個節(jié)點(diǎn)包含12個絕對節(jié)點(diǎn)坐標(biāo)，所以整個板單元包含48個絕對節(jié)點(diǎn)坐標(biāo).

通過引入高階插值函數(shù)，板單元上任意點(diǎn)P的位置在全局坐標(biāo)系下可表示為

r=Se

(2)

式中：e為基于絕對節(jié)點(diǎn)坐標(biāo)法建立的單元節(jié)點(diǎn)坐標(biāo)；S為該單元對應(yīng)的形函數(shù)矩陣，

S=[S1IS2I…S16I]

(3)

I為3×3的單位矩陣，且有

(4)

1.2 板單元質(zhì)量矩陣

將式(2)對時間求導(dǎo)，可得到板單元上任意點(diǎn)的速度表達(dá)式：

(5)

因此板單元的動能可表示為

(6)

式中：ρ為材料密度；V為梁單元體積；M為單元常數(shù)質(zhì)量矩陣，

(7)

1.3 板單元剛度矩陣

基于連續(xù)介質(zhì)力學(xué)，板單元的變形梯度矩陣可表示為

(8)

(9)

式中：

(10)

由于該應(yīng)變張量為對稱的二階張量，因此可將其改寫成向量形式：

(11)

式中：

(12)

對于彈性體，根據(jù)廣義胡克定律，其應(yīng)變張量與應(yīng)力張量存在如下對應(yīng)關(guān)系：

σ=cE:ε

(13)

式中：σ為單元應(yīng)力張量；cE為材料彈性系數(shù)張量.根據(jù)連續(xù)介質(zhì)力學(xué)，對于各向同性的固體材料，彈性系數(shù)張量可以通過拉梅常數(shù)進(jìn)行簡化表達(dá).簡化后的彈性系數(shù)張量為

(14)

單元剛度矩陣可通過應(yīng)變能函數(shù)求得.應(yīng)變能函數(shù)表達(dá)式為

(15)

式中：V0為單元在初始構(gòu)型下的體積.式(15)對單元節(jié)點(diǎn)坐標(biāo)求偏導(dǎo)，可得單元廣義彈性力為

(16)

由此可推導(dǎo)出剛度矩陣表達(dá)式為

K=(λ+2μ)K1+λK2+2μK3

(17)

式中：

(eTSce-1)Sc1]dV

(eTSbe-1)(Sa1+Sc1)]dV+

(eTSfe)Sf1]dV

(18)

1.4 板單元的動力學(xué)方程

基于Hamilton原理建立板單元動力學(xué)方程：

(19)

式中：Wa為廣義外力做功.對方程中各項進(jìn)行展開可得到簡化后的柔性板單元動力學(xué)模型：

(20)

式中：Fa為系統(tǒng)廣義外力陣，可由虛功原理推導(dǎo)得到，

(21)

F為系統(tǒng)所受廣義外力.

2 壓電驅(qū)動單元及復(fù)合柔性板單元建模

為了引入軟體機(jī)器人的驅(qū)動方式，考慮在柔性板單元基礎(chǔ)上增加壓電驅(qū)動梁單元.

2.1 驅(qū)動梁單元變形描述

驅(qū)動梁單元采用三維兩節(jié)點(diǎn)梁單元，其幾何模型如圖2所示.

圖2 三維兩節(jié)點(diǎn)梁單元模型Fig.2 Three-dimensional two-node beam element

梁單元包含Ax和Bx節(jié)點(diǎn)，根據(jù)絕對節(jié)點(diǎn)坐標(biāo)法，其坐標(biāo)定義在全局坐標(biāo)系下，每個節(jié)點(diǎn)包含12個坐標(biāo)，因此單元共包含24個絕對節(jié)點(diǎn)坐標(biāo)：

(22)

通過引入高階插值函數(shù)，可得到梁單元對應(yīng)的單元形函數(shù)矩陣：

(23)

且有

(24)

2.2 耦合單元變形描述

復(fù)合柔性板結(jié)構(gòu)由板單元和梁單元耦合而成，其幾何模型如圖3所示.

板單元與梁單元的耦合通過變形協(xié)調(diào)條件實現(xiàn)，滿足2點(diǎn)基本假設(shè)：① 耦合單元變形前后，板單元與梁單元端面的切向平面始終共面；② 在運(yùn)動過程中，板單元與梁單元的接觸面不發(fā)生相對滑移.基于這2點(diǎn)基本假設(shè)，以圖4所示的AD端面為例，A′、G、G′和Ax4點(diǎn)始終共線，并且G與G′點(diǎn)始終重合.

根據(jù)上述假設(shè)，由變形協(xié)調(diào)條件可得

圖3 柔性梁板耦合單元模型Fig.3 Flexible composite beam-plate element

圖4 耦合單元AD端截面示意圖Fig.4 Cross section at side AD of the coupled element

(25)

G點(diǎn)的位置和斜率坐標(biāo)可通過形函數(shù)矩陣表示為

(26)

同理，對于梁單元上的G′點(diǎn)有

(27)

綜合式(25)～(27)可得

(28)

對于耦合單元BC端面，其示意圖如圖5所示.

圖5 耦合單元BC端截面示意圖Fig.5 Cross section at side BC of the coupled element

同樣，由F和F′點(diǎn)處的變形協(xié)調(diào)條件可得

(29)

綜合式(28)和(29)，有

(30)

(31)

式中：RBP為關(guān)聯(lián)板單元與梁單元的坐標(biāo)轉(zhuǎn)換矩陣.通過該坐標(biāo)轉(zhuǎn)換矩陣，可用板單元的絕對節(jié)點(diǎn)坐標(biāo)描述梁單元節(jié)點(diǎn)坐標(biāo)，在保證耦合單元建模精度的同時，大大地縮減了耦合單元的坐標(biāo)數(shù)量.

2.3 壓電驅(qū)動材料本構(gòu)關(guān)系

本文使用壓電材料作為復(fù)合柔性單元的驅(qū)動器材料.壓電材料本構(gòu)關(guān)系可表示為

(32)

式中：eE為材料壓電耦合系數(shù)；Q為施加的電場強(qiáng)度；D為電位移；ζ為相對介電常數(shù).將本構(gòu)方程展開可得

由材料的對稱性有

c11=c22,c12=c21,c44=c55

c13=c23=c32=c31

e31=e32,ζ11=ζ22

由該本構(gòu)方程可得驅(qū)動單元彈性勢能表達(dá)式：

(33)

電場力做功為

(34)

2.4 復(fù)合板結(jié)構(gòu)動力學(xué)方程

壓電梁驅(qū)動器的質(zhì)量矩陣由梁單元形函數(shù)矩陣推導(dǎo)得出:

(35)

基于變形協(xié)調(diào)條件得出的坐標(biāo)轉(zhuǎn)換矩陣，可得復(fù)合板單元的質(zhì)量矩陣為

MBP=M+(RBP)TMBRBP

(36)

壓電梁驅(qū)動器剛度矩陣可由式(33)對節(jié)點(diǎn)坐標(biāo)2次求偏導(dǎo)得到，

(37)

電場力做功對節(jié)點(diǎn)坐標(biāo)求偏導(dǎo)可得電場力表達(dá)式為

(38)

基于坐標(biāo)轉(zhuǎn)換矩陣可得復(fù)合板單元的剛度矩陣為

KBP=K+(RBP)TKBRBP+(RBP)TKERBP

(39)

復(fù)合板結(jié)構(gòu)動力學(xué)方程的建立同樣基于Hamil-ton原理，其表達(dá)式為

(40)

對式(40)進(jìn)行展開并引入質(zhì)量矩陣式(36)、剛度矩陣式(39)及廣義外力表達(dá)式可得復(fù)合柔性板結(jié)構(gòu)的動力學(xué)方程為

(41)

式(41)所示的代數(shù)微分方程可采用Runge-Kutta方法進(jìn)行數(shù)值求解.本文采用商用軟件Matlab中的ode45函數(shù)進(jìn)行求解，該函數(shù)是一種自適應(yīng)步長的Runge-Kutta數(shù)值解法，適用于求解具有非剛性特性的常微分方程.

3 壓電驅(qū)動板結(jié)構(gòu)數(shù)值仿真分析

3.1 靜力學(xué)特性分析

數(shù)值仿真分析所選用的復(fù)合柔性板結(jié)構(gòu)如圖6所示.板的AD端完全固定，BC端自由，為懸臂結(jié)構(gòu).柔性板長度、寬度和厚度分別為 0.6、0.6 和 0.01 m.壓電梁驅(qū)動器下表面與柔性板上表面貼合，Ax和Bx點(diǎn)分別位于柔性板AD端和BC端中點(diǎn)上方.壓電梁驅(qū)動器長度、寬度和高度分別為 0.6、0.01 和 0.01 m.該結(jié)構(gòu)縱向截面如圖7所示.壓電梁驅(qū)動器上下表面接通電源，從而產(chǎn)生壓電驅(qū)動力.

圖6 復(fù)合柔性懸臂板結(jié)構(gòu)Fig.6 Cantilever plane with actuator element

圖7 壓電驅(qū)動板結(jié)構(gòu)截面示意圖Fig.7 Cross section of the plane with actuator element

本文的靜力學(xué)計算案例中，柔性板的材料參數(shù)：E=100 MPa；ν=0.3；ρ=1 300 kg/m3.壓電梁驅(qū)動器選用壓電陶瓷材料(PZT)，其材料屬性可由靜態(tài)法試驗測得，具體數(shù)值如表1所示.

表1 壓電梁驅(qū)動器材料屬性Tab.1 Material properties of piezoelectric beam

圖8 不同驅(qū)動電壓下柔性板的變形Fig.8 Deformations of the plane at different voltage

靜力學(xué)計算中，在柔性板和驅(qū)動梁尺寸及材料參數(shù)不變的情況下，對驅(qū)動梁分別施加U=10，20，50，100，200 V的電壓.柔性板在不同驅(qū)動電壓值下的變形如圖8所示.由圖可知，驅(qū)動梁上施加不同的電壓值時，柔性板產(chǎn)生不同程度的變形.這是由于壓電材料具有特殊的逆壓電效應(yīng)，施加在壓電驅(qū)動梁上下表面的電壓產(chǎn)生的電場使得梁結(jié)構(gòu)在厚度方向發(fā)生相應(yīng)的結(jié)構(gòu)變形，從而使得與其相連的柔性板隨著壓電驅(qū)動梁的結(jié)構(gòu)變形而發(fā)生彎曲變形，且該變形量隨驅(qū)動電壓的增大而逐漸增大.

基于上述靜力學(xué)計算結(jié)果，在E=100，50，10 MPa的3種不同條件下計算柔性板末端B點(diǎn)Z方向位移(dBZ)隨驅(qū)動電壓變化的情況，結(jié)果如圖9所示.由圖可見：柔性板在壓電驅(qū)動作用下產(chǎn)生彎曲變形，柔性板末端位移量隨驅(qū)動電壓的增大呈增大趨勢；柔性板末端位移增大幅度隨彈性模量的減小而變大.

圖9 不同驅(qū)動電壓下柔性板末端B點(diǎn)Z方向位移Fig.9 Z-displacement of node B at different voltage

3.2 動力學(xué)特性分析

以圖6所示的帶驅(qū)動器的柔性板為例，分析不同電壓參數(shù)及材料參數(shù)對板結(jié)構(gòu)動力學(xué)特性的影響.柔性板長度、寬度和厚度分別取 0.6、0.6 和 0.01 m，驅(qū)動器長度、寬度和高度分別取 0.6、0.01 和 0.01 m.

壓電材料參數(shù)見表1，驅(qū)動電壓分別取50和100 V，柔性板的材料參數(shù)同前(E=100 MPa，ν=0.3，ρ=1 300 kg/m3).仿真時間為 0.1 s，輸出柔性板末端B點(diǎn)X、Y、Z方向位移.由于仿真結(jié)果數(shù)據(jù)量龐大，所以只在仿真時間區(qū)間內(nèi)均勻選取其中的30個點(diǎn)，且保證相鄰點(diǎn)時間間隔相同，得到B點(diǎn)3個方向位移隨時間(τ)變化的曲線如圖10所示.由圖可見：當(dāng)對柔性板施加驅(qū)動電壓時，板末端產(chǎn)生周期性的上下擺動；B點(diǎn)在X和Y方向上位移量較小，而Z方向的位移以近似余弦規(guī)律隨時間變化；驅(qū)動電壓由50 V增加至100 V時，末端位移明顯增大，且末端運(yùn)動規(guī)律相同，振動周期不變.

圖10 驅(qū)動電壓為50和100 V時B點(diǎn)位移隨時間變化曲線Fig.10 Displacement of node B at 50 and 100 V

基于上述動力學(xué)計算結(jié)果，在E=100，50，10 MPa，U=100 V的條件下，對柔性板進(jìn)行計算，以研究不同材料參數(shù)對板結(jié)構(gòu)動力學(xué)特性的影響，計算得到的板末端B點(diǎn)Z方向位移隨時間變化的曲線如圖11所示.由圖可見，隨著柔性板彈性模量的減小，板末端在驅(qū)動電壓作用下產(chǎn)生的位移量增大，且當(dāng)彈性模量發(fā)生變化時，柔性板末端上下振動的幅度以及振動周期均發(fā)生改變.B點(diǎn)在不同彈性模量下的振動幅度(dBmax)及振動周期(T)如表2所示.由表可知，驅(qū)動電壓不變時，柔性板末端B點(diǎn)Z方向位移隨彈性模量的減小呈增大趨勢，且振動周期相應(yīng)增大.

圖11 驅(qū)動電壓為100 V時B點(diǎn)Z方向位移Fig.11 Displacement of node B at 100 V

表2 柔性板B點(diǎn)在不同彈性模量下的振動特性Tab.2 Kinetic statistics of node B under actuation

4 結(jié)論

本文基于絕對節(jié)點(diǎn)坐標(biāo)法，對柔性板單元動力學(xué)方程進(jìn)行了推導(dǎo)；同時引入壓電材料本構(gòu)方程，將梁單元作為驅(qū)動器，通過變形協(xié)調(diào)條件將壓電梁驅(qū)動器和柔性板耦合，形成帶驅(qū)動器的柔性復(fù)合板結(jié)構(gòu)，并推導(dǎo)出該結(jié)構(gòu)的動力學(xué)方程.通過對復(fù)合板結(jié)構(gòu)進(jìn)行靜力學(xué)和動力學(xué)的仿真計算，分析了驅(qū)動電壓參數(shù)及材料參數(shù)對復(fù)合板結(jié)構(gòu)動力學(xué)特性的影響，得到的主要結(jié)論如下：

(1) 靜力學(xué)計算中，當(dāng)驅(qū)動梁上施加電壓時，柔性板在驅(qū)動電壓作用下產(chǎn)生彎曲變形.隨著驅(qū)動電壓的增大，柔性板彎曲變形程度增大，板末端位移量呈近似線性增大的趨勢.

(2) 彈性模量對復(fù)合板結(jié)構(gòu)的靜力學(xué)特性具有一定影響.材料彈性模量越小，結(jié)構(gòu)末端在驅(qū)動電壓作用下產(chǎn)生的位移增幅隨驅(qū)動電壓的增大而增大.相同驅(qū)動電壓條件下，柔性板末端位移隨彈性模量的減小而非線性增大.

(3) 動力學(xué)計算中，當(dāng)驅(qū)動梁上施加持續(xù)電壓時，柔性板末端產(chǎn)生上下周期性的振動.隨著驅(qū)動電壓的增大，板末端振動的幅度相應(yīng)增大，振動規(guī)律基本不變.當(dāng)驅(qū)動電壓一定而彈性模量減小時，板末端的振動幅度增大，同時振動周期也相應(yīng)增長.