柔體系統(tǒng)動(dòng)力學(xué)算法中的形函數(shù)選擇問(wèn)題研究

趙　洋　王德石

（海軍工程大學(xué)兵器工程系　武漢　430033）

1　引言

多柔體動(dòng)力學(xué)研究機(jī)械與結(jié)構(gòu)系統(tǒng)中的構(gòu)件彈性、運(yùn)動(dòng)副變形與剛性運(yùn)動(dòng)的耦合動(dòng)力學(xué)特性，其與多剛體動(dòng)力學(xué)有相同的研究目的，即旨在解決動(dòng)力學(xué)系統(tǒng)的仿真與控制正逆兩類(lèi)問(wèn)題［1］。在發(fā)展剛體系統(tǒng)動(dòng)力學(xué)的過(guò)程中，首先關(guān)心的是實(shí)時(shí)控制與仿真計(jì)算所需要的算法效率，所以利用分析力學(xué)的諸多原理［2～5］，建立了各種遞推算法［2～4］；隨著計(jì)算機(jī)芯片技術(shù)的進(jìn)步，算法的效率不再是研究的重點(diǎn)，并且根據(jù)力學(xué)原理的等價(jià)性，無(wú)論采用何種原理都可以獲得系統(tǒng)的一般形式的封閉動(dòng)力學(xué)方程［1，5］。在機(jī)器人［6］、航天器［7］、運(yùn)動(dòng)生物力學(xué)［8］、以及火炮機(jī)構(gòu)［9～13］等諸多領(lǐng)域，利用多剛體動(dòng)力學(xué)理論與算法，已經(jīng)可以建立所需的動(dòng)力學(xué)模型，并且用于分析系統(tǒng)的響應(yīng)或?yàn)橄到y(tǒng)設(shè)計(jì)實(shí)時(shí)動(dòng)力學(xué)控制算法。

在多剛體動(dòng)力學(xué)基礎(chǔ)上進(jìn)一步考慮系統(tǒng)中構(gòu)件或鉸鏈的變形特性，既是工程中的需求，也是提高仿真與控制算法精度的理論需要。周起釗在文獻(xiàn)［13］中表明：“帶有太陽(yáng)能帆板的航天器在將帆板視為剛體的情況下系統(tǒng)是穩(wěn)定的，而考慮帆板的彈性時(shí)系統(tǒng)就可能是不穩(wěn)定的”，并且在柔性體系動(dòng)力學(xué)發(fā)展之早期，就指出了多柔性體系統(tǒng)動(dòng)力學(xué)的重要研究問(wèn)題［14］。其一是浮動(dòng)坐標(biāo)系的選擇問(wèn)題，它決定了系統(tǒng)動(dòng)力學(xué)方程的耦合程度［1，13～15］；其二是彈性體或鉸鏈變形的描述問(wèn)題，它決定著仿真與控制算法的運(yùn)算效率與計(jì)算精度［1，14］。浮動(dòng)坐標(biāo)系和彈性體的形函數(shù)帶來(lái)的這些問(wèn)題，理論上似乎比較完善，但在工程領(lǐng)域遇到的問(wèn)題是，針對(duì)某一領(lǐng)域的仿真與控制問(wèn)題，究竟如何選取它們至今不夠明確。目前在各個(gè)領(lǐng)域如航天器與火炮的仿真與控制，基本上是給定一種方式建立動(dòng)力學(xué)模型，簡(jiǎn)化處理后，設(shè)計(jì)相關(guān)算法。事實(shí)上針對(duì)不同的應(yīng)用情形，浮動(dòng)坐標(biāo)的建立與彈性體的變形及其控制要求是不盡相同的，尤其是形函數(shù)的選擇，不同的應(yīng)用領(lǐng)域的精度要求決定了假定形函數(shù)的復(fù)雜程度以及計(jì)算效率的高低，例如火炮的身管變形與航天器太陽(yáng)能帆板就不必取同一計(jì)算精度的形函數(shù)。這樣的問(wèn)題至今并未得到足夠的研究。

本文針對(duì)形函數(shù)影響數(shù)值計(jì)算精度和效率問(wèn)題，在分析浮動(dòng)坐標(biāo)系建立于指定鉸鏈處的基礎(chǔ)上，首先建立多柔性體系統(tǒng)一般形式的動(dòng)力學(xué)方程，然后針對(duì)三種形式的形函數(shù)，通過(guò)分析動(dòng)力學(xué)方程中的系數(shù)矩陣，討論不同的形函數(shù)導(dǎo)致的柔性體系統(tǒng)仿真與控制的計(jì)算精度與運(yùn)算效率。最后通過(guò)一個(gè)雙連桿算例對(duì)采用不同類(lèi)型變形函數(shù)的動(dòng)力學(xué)方程進(jìn)行仿真計(jì)算，以期為柔體系統(tǒng)仿真與控制中的形函數(shù)選擇提供參考。

2　多柔體系統(tǒng)運(yùn)動(dòng)學(xué)

考慮系統(tǒng)中任意兩個(gè)鉸鏈連接的柔性體，建立系統(tǒng)的坐標(biāo)系，如圖1所示。慣性坐標(biāo)系位于O點(diǎn)，Oixiyizi為第i個(gè)物體的浮動(dòng)坐標(biāo)系。取浮動(dòng)坐標(biāo)系原點(diǎn)Oi位于第i-1個(gè)柔體與第i個(gè)柔體之間的鉸鏈上，其在慣性坐標(biāo)系中的位置向量為ROi。取第i個(gè)物體為研究對(duì)象，柔性體上任意一點(diǎn)P的位移為剛性運(yùn)動(dòng)ro與彈性運(yùn)動(dòng)uf之疊加。

圖1　兩個(gè)鉸接的柔性體

在動(dòng)坐標(biāo)系中，令uO為未變形時(shí)的位置向量，uf為任意變形引起的位置向量，且令u為P點(diǎn)相對(duì)于動(dòng)坐標(biāo)系的位置向量，則有u=uO+uf。當(dāng)任意的P點(diǎn)遍歷柔體時(shí)，柔體的位移形成一個(gè)位移場(chǎng)Φ(x，y，z)，且該位移場(chǎng)是相容和完備的。為了進(jìn)一步描述物體的變形模式，取形函數(shù)矩陣Φ=[Φ1，Φ2，…，Φn] ，則彈性體上各點(diǎn)的變形量 uf可表示為uf=Φqf，其中：qf為對(duì)應(yīng)于變形的廣義坐標(biāo)。在慣性坐標(biāo)系中，點(diǎn)P的位置矢量可表示為

式中：A為動(dòng)坐標(biāo)系相對(duì)于慣性坐標(biāo)系的旋轉(zhuǎn)變換矩陣。對(duì)式（1）求導(dǎo)，可得P點(diǎn)速度矢量為

由歐拉參數(shù)的定義［15］，旋轉(zhuǎn)變換矩陣可表示為A=E?G?T。設(shè) A?u=-2Au?G?P?=BP?。 P 點(diǎn)速度矢量可表示為矩陣形式：

式（2）確定了任意柔性體上任意一點(diǎn)的速度矢量，可以看出形函數(shù)影響柔性體任意一點(diǎn)的速度矢量和動(dòng)能，形函數(shù)的選擇對(duì)柔性體系統(tǒng)動(dòng)力學(xué)方程的求解至關(guān)重要。

3　動(dòng)力學(xué)基本方程

1）柔性體i的動(dòng)能計(jì)算。由任意P點(diǎn)速度矢量表達(dá)式，可求得柔性體i的動(dòng)能為

2）柔性體i的勢(shì)能計(jì)算。由彈性變形引起的內(nèi)力的虛功為

式中：D*為微分算子矩陣，E為彈性模量矩陣。ε=D*uf=D*Φqf，σ=Eε=ED*Φqf。 Qe=-K為彈性體變形引起的彈性力的廣義力，而Kff=∫v(D*Φ)TED*Φd V 為對(duì)應(yīng)于廣義坐標(biāo) qf的彈性體的剛度矩陣，是對(duì)稱正定矩陣。彈性體的應(yīng)變能為

3）廣義力的計(jì)算。設(shè)作用于柔體上一點(diǎn)的集中力為F=F(q，t)，在其作用點(diǎn)虛位移上的虛功為

通過(guò)推導(dǎo)柔體系統(tǒng)動(dòng)力學(xué)方程的質(zhì)量矩陣M和剛度矩陣K表達(dá)式可以看出，在指定浮動(dòng)坐標(biāo)系下，柔體系統(tǒng)動(dòng)力學(xué)方程形式以及各系數(shù)矩陣元素相同，所以方程的解耦程度相同，說(shuō)明形函數(shù)的選擇與方程的解耦程度無(wú)關(guān)。不同形式的形函數(shù)導(dǎo)致不同復(fù)雜程度的系數(shù)矩陣以及不同數(shù)量的廣義坐標(biāo)，因而形函數(shù)影響動(dòng)力學(xué)方程計(jì)算精度和效率。按照目前采用的假設(shè)方法，雖然可以得到問(wèn)題的解，但是其精度估算是非常困難的，誤差也是時(shí)變函數(shù)，這正是柔體動(dòng)力學(xué)的困難所在。

4　柔性體的形函數(shù)

柔體在運(yùn)動(dòng)過(guò)程中本身的構(gòu)型是不斷變化的，采用離散化手段，把柔性體相對(duì)變形表示成N個(gè)廣義坐標(biāo)qi(i=1，…，N)的線性組合，獲得變形一般位移表達(dá)式為

其中：Φ為一般形式的形函數(shù)，它只與位置有關(guān)。假定合適的形函數(shù)是離散化的核心問(wèn)題，期望選取的形函數(shù)盡可能描述物體的實(shí)際變形，靜力變形、有限元模型和振動(dòng)模態(tài)是主要的三種形函數(shù)。

考慮圖2所示的雙連桿機(jī)構(gòu)，m1和m2為端部集中質(zhì)量，兩連桿均為均勻歐拉梁，其質(zhì)量分別為M1、M2，長(zhǎng)度分別為 l1、l2，剛度分別為 EI1、EI2。下面分別以三種形函數(shù)描述連桿變形。

圖2　雙連桿平面運(yùn)動(dòng)機(jī)構(gòu)

情形1：以靜力變形為形函數(shù)。靜力變形模型是取物體在外力包括重力作用下的靜變形作為其假設(shè)變形模式。連桿可等效為懸臂梁模型，受自身重力和端部集中力作用。由材料力學(xué)知識(shí)可得，連桿在受均布載荷（自身重力）和端部集中力作用下的柔度分別為

采用靜力變形模型描述柔性體運(yùn)動(dòng)中的彈性變形，比使用剛體模型計(jì)算時(shí)的精度有所提高。靜力變形模型簡(jiǎn)單，所需廣義坐標(biāo)數(shù)量少，計(jì)算效率高。但以靜力變形代替動(dòng)力變形描述柔性體彈性變形，計(jì)算精度必然存在較大誤差。

情形2：以有限元模型為形函數(shù)。有限元法是將具有復(fù)雜形狀、邊界條件和載荷的物體化整為零，分割為有限數(shù)量、有限大小且有一定規(guī)則形狀的單元。屬于i物體第 j單元上任意一點(diǎn)P的位移向量可表示為

式中：Nij為 j單元的變形模式或假設(shè)位移場(chǎng)，稱為 j單元的形函數(shù)，為該單元的節(jié)點(diǎn)位移向量。在將所有單元拼裝后，物體上所有節(jié)點(diǎn)的位移向量，就構(gòu)成了該物體的彈性廣義坐標(biāo)。在圖2所示模型中，將連桿等效為梁模型并劃分為i個(gè)單元，單元長(zhǎng)l，單元節(jié)點(diǎn)位移是節(jié)點(diǎn)的撓度和轉(zhuǎn)角：

通過(guò)多項(xiàng)式擬合，可得到以單元節(jié)點(diǎn)位移表示的形函數(shù)：

在得到梁?jiǎn)卧獎(jiǎng)偠染仃?、質(zhì)量矩陣后，根據(jù)初始條件可組裝成系統(tǒng)的剛度矩陣和質(zhì)量矩陣，然后進(jìn)行振動(dòng)分析。在有限元模型計(jì)算中，單元變形模式Nij的選取有多種形式，可根據(jù)柔性體的材料、形狀以及計(jì)算精度和效率限制等因素酌情選取。劃分單元數(shù)量越多，廣義坐標(biāo)數(shù)量越多，計(jì)算精度越高，隨之帶來(lái)的問(wèn)題是計(jì)算效率的下降。

情形3：以振動(dòng)模態(tài)為形函數(shù)。振動(dòng)模態(tài)是結(jié)構(gòu)動(dòng)力學(xué)中普遍采用的一種描述物體空間變形的方法。它通過(guò)引入模態(tài)向量和模態(tài)坐標(biāo)來(lái)描述柔性體上任意點(diǎn)的變形，即

其中：φi(x)為桿件的第i階振型，qi(t)為模態(tài)坐標(biāo)。根據(jù)圖2所示連桿的邊界條件，連桿1、2的振型函數(shù)可取簡(jiǎn)支梁的振型函數(shù)，即(x)=sin。

可以根據(jù)先驗(yàn)的響應(yīng)特征和精度要求，來(lái)考慮模態(tài)截?cái)喾秶?。一般情況下，低階模態(tài)響應(yīng)貢獻(xiàn)大，高階模態(tài)響應(yīng)貢獻(xiàn)較小在保證求解精度的前提下，可把貢獻(xiàn)小的模態(tài)截去，以最大限度地縮減求解規(guī)模。

5　案例分析

設(shè)r1、r2分別為連桿1上點(diǎn)P1和連桿2上點(diǎn)P2在慣性坐標(biāo)系中的位置向量，則系統(tǒng)的動(dòng)能為

1）計(jì)算連桿1的動(dòng)能。

坐標(biāo)系的選取如圖3所示，可得連桿1上任意一點(diǎn)P1的位置向量r1：

圖3　連桿1示意圖

連桿1的動(dòng)能為

式中：v1為連桿1的撓度，連桿1的密度ρ(x)=11+mδ(x-l)。111

2）計(jì)算桿2的動(dòng)能。

坐標(biāo)系的選取如圖4所示，可得連桿2上任意一點(diǎn)P2的位置向量r2：

圖4　連桿2示意圖

連桿2的動(dòng)能為

式中：v2為連桿2的撓度，連桿2的質(zhì)量密度ρ(x)=22+mδ(x-l) 。222

情形1：以靜力更正模型為例進(jìn)行推導(dǎo)，將式（11）、（23）、（26）代入拉格朗日方程，可得質(zhì)量矩陣M中非零元素的表達(dá)式為

由于以靜變形代替桿件動(dòng)變形，系統(tǒng)并未增加廣義坐標(biāo)，剛度矩陣K=0。雙連桿系統(tǒng)重力勢(shì)能為

雙連桿系統(tǒng)的驅(qū)動(dòng)力矩τ1(t)、τ2(t)分別為連桿1和連桿2的主動(dòng)力，利用虛功原理可求得對(duì)應(yīng)的廣義力為

情形2：式（17）、（18）已推導(dǎo)出有限元模型中單元的質(zhì)量矩陣M和剛度矩陣K，將單元質(zhì)量矩陣和剛度矩陣進(jìn)行組裝，即可得到系統(tǒng)的質(zhì)量矩陣和剛度矩陣，此處不再贅述。

情形3：以振動(dòng)模態(tài)為例進(jìn)行推導(dǎo)，取模態(tài)截?cái)鄶?shù)n=3，桿件的前三階振型函數(shù)為

將式（20）、（28）、（23）、（26）代入拉格朗日方程，可得質(zhì)量矩陣M中非零元素的表達(dá)式為

雙連桿系統(tǒng)的驅(qū)動(dòng)力矩τ1(t)、τ2(t)分別為連桿1和連桿2的主動(dòng)力，利用虛功原理可求得對(duì)應(yīng)的廣義力為

分別以靜力變形模型、有限元模型和振動(dòng)模態(tài)描進(jìn)行數(shù)值仿真，仿真參數(shù)如下：

假設(shè)初始時(shí)刻兩桿均處于水平位置，且變形為零，初始條件為q(0)=0 ，q?(0)=0。

取前三階振動(dòng)模態(tài)為參照，分別以靜力變形模型、一階振動(dòng)模態(tài)、有限元模型和剛體模型，計(jì)算連桿角速度和末端位移與其之間的誤差，仿真結(jié)果如圖5～6所示。

圖5～6表示剛體模型以及三種形函數(shù)下柔體模型中的連桿角速度誤差?？梢钥闯觯河捎谠陟o力變形模型中，是以靜力變形代替動(dòng)力變形作為柔體的形函數(shù)，所以靜力變形和剛體模型對(duì)柔體角速度的影響程度相同。剛體模型和靜力變形誤差較大，達(dá)到±0.14rad/s，柔性變形影響桿件的剛性轉(zhuǎn)角，在精度要求較高的運(yùn)動(dòng)機(jī)械中，結(jié)構(gòu)的柔性變形不可忽略。一階振動(dòng)模態(tài)的計(jì)算誤差小于±0.05rad/s，計(jì)算精度明顯提高，對(duì)于一般機(jī)械結(jié)構(gòu)，在保證求解精度的前提下，可截去高階模態(tài)，以最大限度地縮減求解規(guī)模。有限元模型誤差小于±0.02rad/s，計(jì)算精度最高，但由于每個(gè)單元都對(duì)應(yīng)著數(shù)個(gè)廣義坐標(biāo)，系統(tǒng)動(dòng)力學(xué)方程系數(shù)矩陣維數(shù)高，計(jì)算效率最低。

圖5　連桿1角速度誤差

圖6　連桿2角速度誤差

圖7　端點(diǎn)X方向位移誤差

圖8　端點(diǎn)Y方向位移誤差

通過(guò)對(duì)比圖7、8中剛體模型以及三種形函數(shù)下柔體模型中的端點(diǎn)位移誤差，我們可以看出：由于桿件柔性的存在，桿件末端速度與按剛性所推算出來(lái)的末端速度存在偏差，X方向和Y方向的偏差可達(dá)到±0.07m/s。靜力變形模型中，由于考慮了桿件的彈性變形，誤差較剛體模型有所減小，但以靜力變形代替動(dòng)力變形作為柔體的形函數(shù)，計(jì)算誤差相對(duì)于一階振動(dòng)模態(tài)和有限元模型仍然較大。一階振動(dòng)模態(tài)和有限元模型誤差較小，并且隨著模態(tài)截?cái)鄶?shù)增大和劃分單元數(shù)增多，計(jì)算精度提高，但廣義坐標(biāo)也隨之增多，計(jì)算效率降低。此外我們還可以看到，除連桿的變形效應(yīng)是必須加以考慮之外，桿的運(yùn)動(dòng)中包含著高頻部分的振動(dòng)，從而導(dǎo)致桿件的抖動(dòng)，增大了對(duì)桿件的控制難度。

6　結(jié)語(yǔ)

本文針對(duì)形函數(shù)影響數(shù)值計(jì)算精度和效率問(wèn)題，在分析浮動(dòng)坐標(biāo)系建立于指定鉸鏈處的基礎(chǔ)上，通過(guò)建立一般形式的動(dòng)力學(xué)方程，分別討論了形函數(shù)的三種形式并分析了三種形函數(shù)的特點(diǎn)。分析和案例表明：

在指定浮動(dòng)坐標(biāo)系下，柔體系統(tǒng)動(dòng)力學(xué)方程解耦程度相同，即質(zhì)量矩陣、剛度矩陣和阻尼矩陣有相同元素，說(shuō)明形函數(shù)的選擇與方程解耦程度無(wú)關(guān)。柔體系統(tǒng)動(dòng)力學(xué)方程計(jì)算精度與浮動(dòng)坐標(biāo)系和方程耦合程度的關(guān)系不大，耦合程度主要影響方程復(fù)雜度，而形函數(shù)對(duì)方程計(jì)算效率和精度影響較大。

靜力變形模型方程形式簡(jiǎn)單，計(jì)算效率最高，在計(jì)算桿件末段位移時(shí)，誤差較剛體模型減小，計(jì)算角速度時(shí)與剛體模型誤差相同。相比之下，有限元模型應(yīng)用靈活，適合復(fù)雜結(jié)構(gòu)的仿真計(jì)算。振動(dòng)模態(tài)模型中模態(tài)截?cái)鄶?shù)不同，導(dǎo)致計(jì)算結(jié)果不同。有限元模型和振動(dòng)模態(tài)模型中的廣義坐標(biāo)數(shù)量和計(jì)算精度隨著劃分單元個(gè)數(shù)和模態(tài)截?cái)鄶?shù)的增加而增加，計(jì)算效率也隨之下降。

上述工作從理論公式推導(dǎo)與分析，并通過(guò)案例分析的不同側(cè)面，詳細(xì)給出了柔體系統(tǒng)彈性變形選取與方程計(jì)算精度和效率之間的變化規(guī)律，根據(jù)不同需求，可供柔體系統(tǒng)動(dòng)力學(xué)的仿真與控制參考選用。

［1］王德石.火炮振動(dòng)理論［M］.北京：兵器工業(yè)出版社，2015：111-114.WANG Deshi.Theory of Gun vibration［M］.Beijing：Weapons industry press，2015：111-114.

［2］遇立基，陳循介.操作機(jī)器人動(dòng)力學(xué)與算法［M］.北京：機(jī)械工業(yè)出版社，1983：1-5.Popov.，E.L..Dynamics and Its Algorism for Manipulation Robot［M］.Beijing：Mechanical industry press，1983：1-5.

［3］WALKERMW.Efficient dynamic computer simulation of robotic mechanisms［J］.Measurement and Control，1982，104：205-211.

［4］王德石.慣量耦合補(bǔ)償控制的遞推算法［M］.北京：科學(xué)出版社，2002：122-130.WANG Deshi.Recursive algorism for the inertia compensation control［M］.Beijing：Science Press，2002：122-130.

［5］KANE TR，LEVINSON D A.The use of Kane's dynamical equations in robots［M］.New York：Mc Graw-Hill Book Company，1983：97-119.

［6］WITTERBURG J.Dynamics of systems of rigid bodies［M］.Stuttgart：B.G.Teubner，1977：107-111.

［7］ROBE TR，KANE TR，etal.Dynamics ofan elastic satellite［J］.Solid and Structure，1996（3）：333-352，691-703，1031-1051.

［8］HUSTON R L.Human body dynamics：impact，occupational，and athletic aspects［M］.Oxford：Claremont press，1982：59-62.

［9］WANG Deshi.Vibrations analysis for gun systems bymultibody dynamics［M］.New York：Academic Press，1999：111-116.

［10］王德石，史躍東.火炮振動(dòng)分析與多體系統(tǒng)模型研究［J］.動(dòng)力學(xué)與控制學(xué)報(bào)，2012，10（4）：304-323.WANG Deshi，SHIYuedong.Research on Artillery Vibration Analysis and Multi-body System Model［J］.Journalof Dynamicsand Control，2012，10（4）：304-323.

［11］WANG Deshi.Dynamics ofmechanical system solved by quadratic programming［M］.New York：Academic Press，1993：63-67.

［12］芮筱亭.有任意個(gè)集中質(zhì)量的轉(zhuǎn)管炮的固有振動(dòng)［J］.兵工學(xué)報(bào)，1994（2）：1-5.RUIXiaoting.The natural vibration ofa rotating gun with arbitrary lumped mass［J］.Acta Weapon Industry，1994（2）：1-5.

［13］周起釗.柔性系統(tǒng)力學(xué)中的主要課題［J］.力學(xué)進(jìn)展，1989，19（4）：464-476.ZHOUQizhao.Themain topics in themechanics of flexible systems［J］.Advances in Mechanics，1989，19（4）：464-476.

［14］周起釗.多體系統(tǒng)力學(xué)的一個(gè)重要研究方向［R］.北京：北京大學(xué)，1988：1-4.ZHOU Qizhao.An important research direction ofmultibody system dynamics［R］.Beijing：Peking University，1988：1-4.

［15］陸佑方.多柔性體系動(dòng)力學(xué)［M］.北京：科學(xué)出版社，1993.124-126.LU Youfang.Flexible Multi-Body System Dynamics［M］.Beijing：Higher Education Press，1993.124-126.

［16］曾攀.有限元分析及應(yīng)用［M］.北京：清華大學(xué)出版社，2003.25-31.ZENG Pan.Finite element analysis and application［M］.Beijing：Tsinghua university press，2003.25-31.