基于絕對節點坐標法的壓電驅動復合結構動力學特性

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

(上海交通大學 a. 上海市復雜薄板結構數字化制造重點實驗室； b. 機械系統與振動國家重點實驗室， 上海 200240)

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

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

鑒于絕對節點坐標法采用的高階單元模型可對具有大變形特征的柔性復合結構運動過程中的剛柔耦合特性進行精確描述，同時考慮到耦合單元精確建模的需求，本文基于絕對節點坐標法，根據變形協調條件推導坐標關聯矩陣，建立了大柔性梁板耦合單元模型，并對其動力學方程進行了求解.此外，引入壓電材料本構方程，將梁單元作為復合結構的驅動器，用Matlab軟件對壓電驅動柔性板結構進行數值仿真分析，研究該結構在外加驅動下的靜力學和動力學特性，并分析各項參數對其性能的影響，以期為軟體機器人的結構及驅動設計提供參考.

1 基于絕對節點坐標法的柔性單元建模

1.1 單元的變形描述

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

圖1 三維四節點板單元模型Fig.1 Three-dimensional four-node plate element

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

(1)

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

通過引入高階插值函數，板單元上任意點P的位置在全局坐標系下可表示為

r=Se

(2)

式中：e為基于絕對節點坐標法建立的單元節點坐標；S為該單元對應的形函數矩陣，

S=[S1IS2I…S16I]

(3)

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

(4)

1.2 板單元質量矩陣

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

(5)

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

(6)

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

(7)

1.3 板單元剛度矩陣

基于連續介質力學，板單元的變形梯度矩陣可表示為

(8)

(9)

式中：

(10)

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

(11)

式中：

(12)

對于彈性體，根據廣義胡克定律，其應變張量與應力張量存在如下對應關系：

σ=cE:ε

(13)

式中：σ為單元應力張量；cE為材料彈性系數張量.根據連續介質力學，對于各向同性的固體材料，彈性系數張量可以通過拉梅常數進行簡化表達.簡化后的彈性系數張量為

(14)

單元剛度矩陣可通過應變能函數求得.應變能函數表達式為

(15)

式中：V0為單元在初始構型下的體積.式(15)對單元節點坐標求偏導，可得單元廣義彈性力為

(16)

由此可推導出剛度矩陣表達式為

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

(17)

式中：

(eTSce-1)Sc1]dV

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

(eTSfe)Sf1]dV

(18)

1.4 板單元的動力學方程

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

(19)

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

(20)

式中：Fa為系統廣義外力陣，可由虛功原理推導得到，

(21)

F為系統所受廣義外力.

2 壓電驅動單元及復合柔性板單元建模

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

2.1 驅動梁單元變形描述

驅動梁單元采用三維兩節點梁單元，其幾何模型如圖2所示.

圖2 三維兩節點梁單元模型Fig.2 Three-dimensional two-node beam element

梁單元包含Ax和Bx節點，根據絕對節點坐標法，其坐標定義在全局坐標系下，每個節點包含12個坐標，因此單元共包含24個絕對節點坐標：

(22)

通過引入高階插值函數，可得到梁單元對應的單元形函數矩陣：

(23)

且有

(24)

2.2 耦合單元變形描述

復合柔性板結構由板單元和梁單元耦合而成，其幾何模型如圖3所示.

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

根據上述假設，由變形協調條件可得

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

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

(25)

G點的位置和斜率坐標可通過形函數矩陣表示為

(26)

同理，對于梁單元上的G′點有

(27)

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

(28)

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

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

同樣，由F和F′點處的變形協調條件可得

(29)

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

(30)

(31)

式中：RBP為關聯板單元與梁單元的坐標轉換矩陣.通過該坐標轉換矩陣，可用板單元的絕對節點坐標描述梁單元節點坐標，在保證耦合單元建模精度的同時，大大地縮減了耦合單元的坐標數量.

2.3 壓電驅動材料本構關系

本文使用壓電材料作為復合柔性單元的驅動器材料.壓電材料本構關系可表示為

(32)

式中：eE為材料壓電耦合系數；Q為施加的電場強度；D為電位移；ζ為相對介電常數.將本構方程展開可得

由材料的對稱性有

c11=c22,c12=c21,c44=c55

c13=c23=c32=c31

e31=e32,ζ11=ζ22

由該本構方程可得驅動單元彈性勢能表達式：

(33)

電場力做功為

(34)

2.4 復合板結構動力學方程

壓電梁驅動器的質量矩陣由梁單元形函數矩陣推導得出:

(35)

基于變形協調條件得出的坐標轉換矩陣，可得復合板單元的質量矩陣為

MBP=M+(RBP)TMBRBP

(36)

壓電梁驅動器剛度矩陣可由式(33)對節點坐標2次求偏導得到，

(37)

電場力做功對節點坐標求偏導可得電場力表達式為

(38)

基于坐標轉換矩陣可得復合板單元的剛度矩陣為

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

(39)

復合板結構動力學方程的建立同樣基于Hamil-ton原理，其表達式為

(40)

對式(40)進行展開并引入質量矩陣式(36)、剛度矩陣式(39)及廣義外力表達式可得復合柔性板結構的動力學方程為

(41)

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

3 壓電驅動板結構數值仿真分析

3.1 靜力學特性分析

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

圖6 復合柔性懸臂板結構Fig.6 Cantilever plane with actuator element

圖7 壓電驅動板結構截面示意圖Fig.7 Cross section of the plane with actuator element

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

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

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

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

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

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

3.2 動力學特性分析

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

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

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

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

圖11 驅動電壓為100 V時B點Z方向位移Fig.11 Displacement of node B at 100 V

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

4 結論

本文基于絕對節點坐標法，對柔性板單元動力學方程進行了推導；同時引入壓電材料本構方程，將梁單元作為驅動器，通過變形協調條件將壓電梁驅動器和柔性板耦合，形成帶驅動器的柔性復合板結構，并推導出該結構的動力學方程.通過對復合板結構進行靜力學和動力學的仿真計算，分析了驅動電壓參數及材料參數對復合板結構動力學特性的影響，得到的主要結論如下：

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

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

(3) 動力學計算中，當驅動梁上施加持續電壓時，柔性板末端產生上下周期性的振動.隨著驅動電壓的增大，板末端振動的幅度相應增大，振動規律基本不變.當驅動電壓一定而彈性模量減小時，板末端的振動幅度增大，同時振動周期也相應增長.