卓英鵬，王 剛，齊朝暉，張 健

（1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023；2.大連理工大學(xué) 海洋科學(xué)與技術(shù)學(xué)院，遼寧 盤錦 124221）

引 言

實(shí)際應(yīng)用中許多結(jié)構(gòu)都可以采用梁的理論分析，例如大型可展開天線、柔性機(jī)械臂、渦輪螺旋槳等.但是隨著材料的輕質(zhì)化以及結(jié)構(gòu)細(xì)長的特點(diǎn)，它們在工作過程中往往會(huì)產(chǎn)生很大的撓度.因此，如何構(gòu)造出一種簡潔高效的空間幾何非線性梁單元是很多領(lǐng)域的共同需要[1-5].

當(dāng)梁的截面發(fā)生大范圍的剛體運(yùn)動(dòng)時(shí)，小變形理論中的線性應(yīng)變不再適用，很多專家學(xué)者將Green 應(yīng)變中的高階項(xiàng)考慮進(jìn)來以便描述幾何非線性效應(yīng)，即早期常用的基于Langrange 坐標(biāo)系的TL 和UL 列式法，但是在推導(dǎo)過程中十分繁瑣，且在分析計(jì)算過程中，經(jīng)常忽略變形后的坐標(biāo)轉(zhuǎn)換矩陣，導(dǎo)致在大轉(zhuǎn)角問題分析中無法收斂甚至出錯(cuò)[6-8].這就引發(fā)了一個(gè)問題：如何可以在利用現(xiàn)有成熟的線性單元基礎(chǔ)上，采用一種更加簡潔形象的策略來描述梁在大變形過程中截面的剛體運(yùn)動(dòng)？因此，Wempner[9]和Belytschko 等[10]最早提出了共旋坐標(biāo)法，其主要思想是將單元的位移場分解為隨單元坐標(biāo)系的剛體運(yùn)動(dòng)以及相對于單元的小位移，即將構(gòu)件的位形認(rèn)為是單元坐標(biāo)系的大范圍運(yùn)動(dòng)與相對于該坐標(biāo)系變形的疊加，這類單元在具體應(yīng)用中展現(xiàn)出了明顯的優(yōu)勢，在此基礎(chǔ)上結(jié)合子結(jié)構(gòu)法的使用和線性系統(tǒng)自由度凝聚理論，解決了實(shí)際工程中的很多問題[11-13]，但是這種策略未考慮單元水平上的幾何非線性問題，用于高速輕質(zhì)系統(tǒng)的分析時(shí)經(jīng)常遇到諸如動(dòng)力剛化項(xiàng)遺失[14-15]、數(shù)值不穩(wěn)定等問題.因此，如何構(gòu)造一種幾何非線性單元愈發(fā)值得重視.

針對幾何非線性單元，Simo 等提出了一種精確幾何模型方法[16-17]，依據(jù)有限轉(zhuǎn)動(dòng)理論導(dǎo)出了具有客觀性的應(yīng)變度量，徹底摒棄了小位移小轉(zhuǎn)動(dòng)的線性梁理論，在處理幾何非線性梁問題時(shí)具有較好的計(jì)算效率和精度[18-20]；Romero 等[21]和Crisfield 等[22]以形心平移和截面轉(zhuǎn)動(dòng)參數(shù)作為節(jié)點(diǎn)坐標(biāo)，對轉(zhuǎn)動(dòng)插值進(jìn)行研究，轉(zhuǎn)動(dòng)和位移插值的獨(dú)立性引起了諸如運(yùn)動(dòng)學(xué)描述冗余和剪切閉鎖等困難；因此Zupan 和Saje 等[23-25]放棄了對轉(zhuǎn)動(dòng)參數(shù)的插值，以應(yīng)變矢量作為節(jié)點(diǎn)坐標(biāo)，對拉伸應(yīng)變和彎曲應(yīng)變進(jìn)行獨(dú)立插值；上述非線性單元均為Timoshenko 梁，解決細(xì)長結(jié)構(gòu)時(shí)沒有體現(xiàn)Euler-Bernoulli 梁的變形耦合關(guān)系.Shabana 提出了一種絕對節(jié)點(diǎn)坐標(biāo)法，將單元節(jié)點(diǎn)的坐標(biāo)定義在全局坐標(biāo)系下，采用斜率矢量代替?zhèn)鹘y(tǒng)有限元中的節(jié)點(diǎn)轉(zhuǎn)角坐標(biāo)，可以避免轉(zhuǎn)角插值，推導(dǎo)的多體系統(tǒng)微分-代數(shù)方程具有常質(zhì)量陣、不存在科氏力和離心力項(xiàng)等特點(diǎn)[26-27]，尤其是近幾年，迅速地得到了廣泛的應(yīng)用，很多基于絕對節(jié)點(diǎn)坐標(biāo)的梁、板、殼單元應(yīng)運(yùn)而生，這種單元雖然比傳統(tǒng)有限元有優(yōu)勢[28-30]，但其節(jié)點(diǎn)參數(shù)較多，無疑對計(jì)算效率會(huì)有影響.

此外，成熟的有限單元建模理論大都仍采用位移元進(jìn)行離散[31]，即直接選擇節(jié)點(diǎn)位移和轉(zhuǎn)角作為描述參數(shù)，這種以虛功原理為基礎(chǔ)的位移元可以保證單元間位移的連續(xù)，一般應(yīng)力對應(yīng)于位移的二階項(xiàng)，理論上，位移元建模單元間的應(yīng)力仍然處于不連續(xù)狀態(tài)，端部施加外力時(shí)，力的邊界條件近似滿足，這些情況只能通過不斷縮小單元長度逐漸逼近連續(xù)和降低近似.大變形梁剛性截面在做大范圍剛體運(yùn)動(dòng)時(shí)，傳統(tǒng)位移元的應(yīng)力不連續(xù)性以及力的邊界條件近似性可能會(huì)在單元數(shù)量選擇不合適的情況下降低分析精度.

因此，一個(gè)值得研究的問題是，是否可以構(gòu)造一種幾何非線性梁單元使其同時(shí)滿足以下幾個(gè)條件：①避免形心位移和轉(zhuǎn)角獨(dú)立插值造成的剪切閉鎖等數(shù)值困難；②在較少自由度情況下，保證單元間應(yīng)力的連續(xù)性以及精確滿足力的邊界條件；③保證一定精度的同時(shí)，加快大變形梁的計(jì)算效率；④采用絕對坐標(biāo)以便于在柔性多體系統(tǒng)動(dòng)力學(xué)各領(lǐng)域的應(yīng)用.鑒于此，本文提出了一種可滿足上述條件的空間幾何非線性樣條梁單元，該單元滿足 Bernoulli-Euler 梁變形耦合關(guān)系，其廣義應(yīng)變與梁截面的剛體運(yùn)動(dòng)無關(guān)，可描述其幾何非線性效應(yīng)；樣條單元的組裝在滿足內(nèi)部節(jié)點(diǎn)應(yīng)力連續(xù)的同時(shí)，縮減了系統(tǒng)自由度；邊界節(jié)點(diǎn)參數(shù)包含應(yīng)變參數(shù)，更加便于施加外力邊界條件；降噪后的運(yùn)動(dòng)方程可采用通用的微分求解器快速求解；最后，通過數(shù)值算例驗(yàn)證了所提單元的有效性.

1 梁單元位移形函數(shù)

如圖1所示，空間梁單元的兩個(gè)節(jié)點(diǎn)分別為左右端面的形心，它們的矢徑及其對弧長的一階二階導(dǎo)數(shù)選作為節(jié)點(diǎn)參數(shù)的一部分.

圖1 空間梁單元Fig.1 A spatial beam element

形心線上的任意點(diǎn)矢徑可以用五次多項(xiàng)式插值擬合得到

為簡化符號(hào)，除非另有說明，本文中變量上標(biāo)“′”表示變量對弧長坐標(biāo)s的偏導(dǎo)數(shù).式中的形函數(shù)

其中L為單元長度，歸一化參數(shù)ξ 定義為

對式(1)求時(shí)間導(dǎo)數(shù)，可得

對式(2)求時(shí)間導(dǎo)數(shù)，可得

對式(3)求時(shí)間導(dǎo)數(shù)，可得

2 梁單元的轉(zhuǎn)動(dòng)描述

廣泛采用的Euler 梁理論通過將截面視為剛性截面，將三維問題轉(zhuǎn)化為一維問題.為了描述截面的運(yùn)動(dòng)，引入一個(gè)固結(jié)于截面形心的坐標(biāo)系，第一根軸es沿截面的法線方向，其余兩根軸et,eb分別指向截面的形心主軸方向，如圖2所示.

圖2 梁單元截面坐標(biāo)系Fig.2 The cross section coordinate system of the beam element

截面坐標(biāo)系的位置和方位是時(shí)間t和形心線弧長坐標(biāo)s的函數(shù)，它們的基矢量保持單位正交性，其對時(shí)間的導(dǎo)數(shù)可以用角速度表示：

類比角速度，引入曲率來描述截面基矢量對弧長坐標(biāo)s的導(dǎo)數(shù)：

根據(jù)這樣的定義，曲率和角速度在截面坐標(biāo)系下分解為

式中，分解系數(shù)為

對式(14)求弧長導(dǎo)數(shù)和時(shí)間導(dǎo)數(shù)，可得

以及

將式(15)、(16)代入到式(13)知，角速度 ω的弧長導(dǎo)數(shù)可用曲率分量的時(shí)間導(dǎo)數(shù)表示：

曲率κ的時(shí)間導(dǎo)數(shù)可用角速度分量的弧長導(dǎo)數(shù)表示：

引入Euler-Bernoulli 梁假設(shè)：梁的截面法線方向始終與形心線的切線方向重合，即

其中，弧長伸縮率

對式(19)求時(shí)間的導(dǎo)數(shù)，可得

采用Cardan 角描述剛性截面的轉(zhuǎn)動(dòng)，將截面坐標(biāo)系相對于總體坐標(biāo)系 {g1,g2,g3}的轉(zhuǎn)動(dòng)分為三次相繼定軸轉(zhuǎn)動(dòng)，轉(zhuǎn)角依次為α,β,γ，則截面法向矢量可以表示為

因此，前兩次轉(zhuǎn)動(dòng)的Cardan 角為

可以看出，Cardan 描述的轉(zhuǎn)動(dòng)對梁的截面運(yùn)動(dòng)有明確的物理意義，即前兩次轉(zhuǎn)動(dòng)完全由梁的形心線的空間方位決定，第三次轉(zhuǎn)動(dòng)是繞著形心線的切線方向完成的，可將 α,β理解為繞截面的形心主軸方向的彎曲角度，將γ理解為扭轉(zhuǎn)角度，避免了利用轉(zhuǎn)動(dòng)矢量進(jìn)行插值來描述截面轉(zhuǎn)動(dòng).

對式(22)求時(shí)間導(dǎo)數(shù)，可得

其中

它們與es相互正交，基于此關(guān)系，容易得到Cardan 角對時(shí)間的導(dǎo)數(shù)為

對式(19)、(20)求時(shí)間導(dǎo)數(shù)，可得

由圖5可知:改變條件后PID控制的超調(diào)更大且達(dá)到穩(wěn)態(tài)所需的時(shí)間更長。對于ADRC控制，外加的干擾對其影響很小。

因此，式(26)可改寫為

矢量h,b的時(shí)間導(dǎo)數(shù)可表示為

其中

同理，Cardan 角的弧長導(dǎo)數(shù)為

它們的時(shí)間導(dǎo)數(shù)為

其中

因此，式(33)可改寫為

對式(29)求時(shí)間導(dǎo)數(shù)，可得

為了描述截面繞法線方向es的扭轉(zhuǎn)，第三個(gè)角γ 被引入，采用三次多項(xiàng)式插值得到

它的形函數(shù)為

其中L為單元長度，歸一化參數(shù)ξ 定義為

對式(38)、(39)求時(shí)間導(dǎo)數(shù)，可得

因此，截面坐標(biāo)系形心主軸基矢量可用Cardan 角描述為

由上述可知，截面的轉(zhuǎn)動(dòng)可以采用Cardan 角描述，前兩個(gè)Cardan 角可以表示為形心矢徑的函數(shù)，第三個(gè)Cardan 角及其弧長導(dǎo)數(shù)可以在單元域內(nèi)采用三次Hermite 插值得到.

3 梁單元的廣義應(yīng)變

如圖3所示，從柔性梁中取出一小段原始弧長為ds的微元體，合力 -f和合力矩 -m作用在它的左截面.

圖3 梁的微元體Fig.3 An infinitesimal beam unit

由剛體動(dòng)力學(xué)方程可得微元體的運(yùn)動(dòng)方程為

根據(jù)虛功率原理：外力虛功率為慣性力虛功率與變形虛功率之和，可得微元體的變形虛功率為

將式(46)代入到式(47)，可改寫為

將式(17)和(21)代入到式(48)，可進(jìn)一步簡寫為

其中力和力矩分量來源于

這表明，在截面剛性假設(shè)以及法向矢量與形心切矢共線的Bernoulli 假設(shè)情況下，梁的變形虛功率中的廣義應(yīng)變和應(yīng)力應(yīng)該表示為

相應(yīng)的本構(gòu)關(guān)系為

它的標(biāo)準(zhǔn)形式通常可改寫為

式中，軸向應(yīng)變 εs、曲率在截面坐標(biāo)系中的分量 κs,κt,κb可認(rèn)為是一組廣義應(yīng)變，它們都與梁的剛體運(yùn)動(dòng)無關(guān)，因而適用于幾何非線性問題.

4 梁單元的角速度和曲率

定軸轉(zhuǎn)動(dòng)的角速度矢量的模等于轉(zhuǎn)角的變化率，方向與轉(zhuǎn)軸矢量方向相同.利用角速度疊加原理，采用Cardan 角描述轉(zhuǎn)動(dòng)的截面角速度矢量為

其中

因此，它在截面坐標(biāo)下的分量為

根據(jù)式(11)、(12)中角速度和曲率的定義，它們之間唯一的區(qū)別在于Cardan 角要么隨時(shí)間變化要么隨弧長變化，因此參考式(56)，可得曲率分量為

對式(56)、(57)求時(shí)間導(dǎo)數(shù)，可得

其中，轉(zhuǎn)換矩陣Tωκ的時(shí)間導(dǎo)數(shù)為

為了列式更加簡潔，三個(gè)Cardan 角組成向量

據(jù)此，式(59)可改寫為

其中

定義節(jié)點(diǎn)參數(shù)矢量組成的列向量：

從上述的分析可以看出，θ是節(jié)點(diǎn)參數(shù)q的函數(shù)，因此，它的時(shí)間導(dǎo)數(shù)可表示為

其中，轉(zhuǎn)換矩陣Tθq，Tθ′q僅是q的函數(shù)，向量aθ，aθ′與加速度項(xiàng)無關(guān).據(jù)此，角速度和曲率分量以及它們的時(shí)間導(dǎo)數(shù)可寫為

其中

從式(58)、(59)和(36)可知，計(jì)算角速度和曲率的時(shí)間導(dǎo)數(shù)涉及矢徑對弧長的二階導(dǎo)數(shù)，為保證單元間應(yīng)力連續(xù)，意味著在單元邊界點(diǎn)處需要滿足矢徑對弧長導(dǎo)數(shù)的二階連續(xù)性，這是選取五次多項(xiàng)式插值作為單元位移形函數(shù)的重要原因.

5 梁單元的質(zhì)量陣和力陣

根據(jù)上述的討論，單元內(nèi)任意點(diǎn)矢徑和應(yīng)變的時(shí)間導(dǎo)數(shù)可以表示為以下的矩陣形式：

單元的慣性力虛功率為

其中，質(zhì)量陣、力陣為

變形虛功率為

其中，彈性節(jié)點(diǎn)力為

經(jīng)過單元組裝，可得梁系統(tǒng)的虛功率方程為

其中，廣義位移u是所有節(jié)點(diǎn)參數(shù)組成的列向量，fin,fa分別為廣義內(nèi)力和廣義外力.為了進(jìn)一步獲得系統(tǒng)動(dòng)力學(xué)方程，往往需要繼續(xù)分析參數(shù)u的獨(dú)立性.

6 梁單元間的連接光滑化

上述構(gòu)造的梁單元每個(gè)節(jié)點(diǎn)上有11個(gè)節(jié)點(diǎn)參數(shù)，眾多的自由度無疑會(huì)影響計(jì)算效率.如圖4所示，為了縮減自由度，n-1個(gè)等長梁單元被組集為一個(gè)樣條單元，內(nèi)部節(jié)點(diǎn)矢徑對弧長的導(dǎo)數(shù)r′k,r′k′（k≠1,k≠n）可以不作為單元的節(jié)點(diǎn)參數(shù).它們是通過節(jié)點(diǎn)處矢徑的三四階弧長導(dǎo)數(shù)連續(xù)性要求確定（五次樣條插值），即根據(jù)式(1)～(7)，存在關(guān)系：

圖4 縮減節(jié)點(diǎn)參數(shù)后的樣條單元Fig.4 Parameter reduction of spline elements

其中，單元長度

方程(75)寫成矩陣形式：

式中，r′=[r′1r′2···r′n]為節(jié)點(diǎn)矢徑對弧長坐標(biāo)一階導(dǎo)數(shù)的矩陣；r′′=[r′1′r′2′···r′n′]為節(jié)點(diǎn)矢徑對弧長坐標(biāo)二階導(dǎo)數(shù)的矩陣；r=[r1r2···rn]為節(jié)點(diǎn)矢徑的矩陣；系數(shù)矩陣都是(n+1)×(n+1)維三對角矩陣，可以表示為

其中，矩陣

所有的系數(shù)矩陣均為常矩陣，對式(81)求導(dǎo)，可得

這樣，可以采用樣條單元直接建立梁結(jié)構(gòu)的整體方程，內(nèi)部節(jié)點(diǎn)處矢徑對弧長的導(dǎo)數(shù)可以不作為節(jié)點(diǎn)參數(shù)，縮減系統(tǒng)的自由度，并且隨著內(nèi)部單元的增多，優(yōu)勢越明顯.

7 樣條梁單元邊界節(jié)點(diǎn)參數(shù)轉(zhuǎn)換

通過式(57)，Cardan 角的弧長導(dǎo)數(shù)可以表示為

因此，邊界節(jié)點(diǎn)處可以使用軸向應(yīng)變的弧長導(dǎo)數(shù) ε′s、Cardan 角α,β,γ 以及廣義應(yīng)變?chǔ)舠,κs,κt,κb代替矢徑對弧長的導(dǎo)數(shù)作為節(jié)點(diǎn)參數(shù)，如圖5所示.這樣做的主要優(yōu)勢在于可以采用式(85)中簡單的約束來精確滿足力的邊界條件，便于提高數(shù)值分析精度.

圖5 樣條單元的邊界節(jié)點(diǎn)參數(shù)Fig.5 Boundary nodal parameters in spline elements

8 樣條梁單元的降噪快速算法

梁結(jié)構(gòu)在外力作用下，剛性截面發(fā)生大范圍剛體運(yùn)動(dòng)的同時(shí)，還會(huì)有微幅高頻的彈性運(yùn)動(dòng)，從常微分方程分類角度來看，由式(74)得到的梁結(jié)構(gòu)動(dòng)力學(xué)方程屬于其解為慢變與快變分量混合的剛性微分方程.常用的剛性微分方程求解器ODE45 中為了準(zhǔn)確求出高頻振動(dòng)，往往需要將步長調(diào)至很短，這極大地增加了計(jì)算代價(jià)；相關(guān)數(shù)值方法的主要思想是在積分格式中引入數(shù)值阻尼濾掉系統(tǒng)的高頻分量，從而加快求解速度，其實(shí)這個(gè)過程可以在建模過程通過將應(yīng)力替換為極短時(shí)間內(nèi)的平均應(yīng)力來完成[32].

在t+τ時(shí)刻，單元應(yīng)力可以近似表示為

它在極短時(shí)間區(qū)間(t,t+Δt)內(nèi)的平均值為

上述分析中梁單元的變形虛功率改寫為

式中

其中，fe由式(73)定義，加速度無關(guān)項(xiàng)aε來源于

從式(89)可見，采用極短時(shí)間區(qū)間內(nèi)平均應(yīng)力代替瞬時(shí)應(yīng)力后，系統(tǒng)方程中增加的阻尼項(xiàng)和內(nèi)力項(xiàng)會(huì)保持低頻分量不變，而快速衰減模型中的高頻分量.光滑因子的選擇應(yīng)該綜合計(jì)算效率和模型精度后進(jìn)行合理確定，根據(jù)大量的計(jì)算經(jīng)驗(yàn)，光滑因子的推薦值為0.000 1＜Δt＜0.01.

9 數(shù)值算例

為檢驗(yàn)所提單元的合理性，分析了5個(gè)算例，其中算例1、2 檢驗(yàn)了模型的濾除高頻成分，加快計(jì)算效率的可行性，以及柔性部件變形大小對單元的有效性；算例3 驗(yàn)證了單元內(nèi)部廣義應(yīng)變的連續(xù)性；算例4、5 證明了單元適用于任意空間大轉(zhuǎn)動(dòng)問題，且與解析解、相關(guān)文獻(xiàn)進(jìn)行了對比.

9.1 算例 1：單擺機(jī)構(gòu)自由下落

L=5 mA=π ×10-2m2ρ=7 850 kg/m3

如圖6所示，擺長為，截面積均為，材料密度為，彈性模量為E=2.1×1011Pa的機(jī)構(gòu)在y軸方向g=-9.8 m/s2重力加速度下，從圖示靜止位置落下.

圖6 自由下落的柔性單擺機(jī)構(gòu)Fig.6 The free-falling flexible pendulum mechanism

應(yīng)用本文所提單元將單擺作為一個(gè)樣條梁單元，內(nèi)部劃分4個(gè)小單元，采用ODE45 求解器，其中光滑因子分別取 Δt=0,0.001,0.01,0.1進(jìn)行數(shù)值積分，取精度控制參數(shù)相 對誤差 εr=1×10-3和 絕對誤差εa=1×10-6.得到單擺機(jī)構(gòu)上末端點(diǎn)沿y和z方向的速度變化，如圖7所示.

圖7 末端點(diǎn)沿y 和z 方向的速度Fig.7 Velocities along y and z direction at the free end

由圖7可以看出，當(dāng)柔性單擺的材料剛度較大時(shí)（E=2.1×1011Pa），柔性梁的變形很小，當(dāng)不使用快速降噪方法時(shí)（Δt=0），末端點(diǎn)的速度具有高頻振蕩分量，致使計(jì)算效率低下.采用本文所提降噪快速算法可使用非剛性微分方程求解器ODE45 求解該問題，且能夠更好地濾除系統(tǒng)中的高頻，得到更為平滑的變形曲線.

表1對比了上述不同方法的計(jì)算時(shí)間，結(jié)果顯示本文所提快速算法可適用于普通ODE 求解器，且可大幅度提高計(jì)算效率.

表1 效率比較（柔性單擺機(jī)構(gòu)）Table 1 Efficiency comparison (for the flexible pendulum mechanism)

9.2 算例 2：端部受剪力或力矩的懸臂梁

如圖8所示，懸臂梁在端部受到集中剪力F或彎矩M作用，結(jié)構(gòu)參數(shù)采用無量綱記法，梁的長度L=10，線密度 ρA=1，抗拉模量EA=104，抗彎抗扭模量GJ=EIy=EIz=103，截面初始構(gòu)型慣性矩Jp=diag(20,10,10)，梁的端部位移為ux,uy，來源于外力邊界條件的受力端節(jié)點(diǎn)的應(yīng)變約束為

圖8 懸臂梁自由端受力矩作用Fig.8 A cantilever beam subjected to an end moment

① 當(dāng)端部受力矩M作用時(shí)，εs=0,κs=0,κt=M/(EIz),κb=0；

② 當(dāng)端部受剪力F作用時(shí)，εs=0,κs=0,κt=0,κb=0.

這里將梁作為一個(gè)樣條單元，內(nèi)部取3個(gè)節(jié)點(diǎn)進(jìn)行計(jì)算，ODE45 求解器參數(shù)設(shè)置為相對誤差εr=1×10-3，絕對誤差 εa=1×10-6.當(dāng)端部受彎矩作用，懸臂梁發(fā)生純彎曲，存在解析解 κ=M/(EIz)，即發(fā)生純彎曲時(shí)彎曲曲率κ 與外力矩M成正比，外力矩M=2πL-1EIz時(shí)，梁彎成一個(gè)圓.表2和圖9、10 分別給出了懸臂梁受彎矩M=2πL-1EI作用時(shí)選取不同光滑因子 Δt仿真所用時(shí)間及得到的端部位移隨時(shí)間變化的曲線.

表2 懸臂梁受集中力矩 M=2πL-1EI 時(shí)選取不同光滑因子仿真用時(shí)Table 2 The simulation time of different smoothing factors for the cantilever beam subjected to an end bending load M=2πL-1EI

圖9 懸臂梁受集中力矩作用的位移ux 曲線Fig.9 End displacement ux of the cantilever beam with an end bending moment

圖10 懸臂梁受集中力矩作用的位移uy 曲線Fig.10 End displacement uy of the cantilever beam with an end bending moment

從圖中可知，當(dāng)光滑因子 Δt為0.001 時(shí)，梁的端部位移與不采用降噪方法（Δt=0）時(shí)計(jì)算的結(jié)果相吻合，最大絕對誤差在0.1 以內(nèi)，但耗時(shí)從3 653 s 縮短至982 s，若對精度要求不高，繼續(xù)增大光滑因子，耗時(shí)將進(jìn)一步縮短.圖11、12 給出了在不同彎矩作用下梁最終平衡時(shí)的變形狀態(tài)及其端部位移與解析解的對比，文獻(xiàn)[33-34]已經(jīng)指出形心位移和轉(zhuǎn)角獨(dú)立插值的傳統(tǒng)梁，解決諸如此類發(fā)生大彎曲變形細(xì)長梁時(shí)，會(huì)出現(xiàn)剪切閉鎖造成誤差增大的現(xiàn)象，對比曲線顯示本文單元計(jì)算結(jié)果與解析解吻合較好，無剪切閉鎖的發(fā)生.

圖11 懸臂梁在不同彎矩作用下的變形曲線Fig.11 Large deformation curves of the cantilever beam under different bending moments

圖12 懸臂梁受集中力矩的位移曲線Fig.12 The displacements of the cantilever beam with an end bending moment

當(dāng)端部受剪力作用處于平衡狀態(tài)時(shí)，小變形范圍內(nèi)端部位移存在解析解uy=FL3(3EIz)-1，大變形結(jié)果可采用ANSYS 軟件打開“NLGEOM”命令給出合理的解.表3和圖13給出了采用本文單元建模，不同剪力作用下梁最終平衡狀態(tài)下端部豎向位移，以及和解析解、ANSYS 結(jié)果的對比.由表3和圖13可知，剪力大小在1～3 N 范圍內(nèi)時(shí)，梁的變形屬于小變形，與解析解的誤差在1%以內(nèi)；隨著剪力的增大，梁的變形逐漸進(jìn)入大變形狀態(tài)，在36 N 時(shí)變形位移已經(jīng)超出梁長的1/2，此階段，解析解已無法滿足要求，本文計(jì)算結(jié)果與ANSYS 大變形計(jì)算結(jié)果誤差隨著剪力增大略微增大，但均在0.13%以內(nèi).

表3 懸臂梁不同剪力作用下平衡狀態(tài)時(shí)豎向位移Table 3 Equilibrium deflections at the cantilever beam end under different shear forces

圖13 懸臂梁平衡狀態(tài)時(shí)豎向位移Fig.13 Deflections at the end in equilibrium

從上述分析結(jié)果可以看出，柔性梁無論處于小變形范圍內(nèi)還是發(fā)生大變形，乃至將梁彎成一個(gè)圓，本文所提單元計(jì)算的結(jié)果都可以得到準(zhǔn)確的解，由此說明本文單元處理柔性梁部件變形的有效性.

9.3 算例 3：端部受剪力作用的旋轉(zhuǎn)空間柔性梁

大范圍運(yùn)動(dòng)的旋轉(zhuǎn)柔性梁已被很多學(xué)者用來考察所提方法的有效性.如圖14所示，梁的左端可以繞z軸自由旋轉(zhuǎn)，右端施加y和z方向上兩個(gè)隨時(shí)間變化的動(dòng)態(tài)載荷Fy(t)和Fz(t).梁的參數(shù)采用無量綱記法，長度L=10，線密度 ρA=1，抗拉模量EA=104，抗彎抗扭模量GJ=EIy=EIz=5×102，截面初始構(gòu)型慣性矩Jp=diag(20,10,10).光滑因子取0.001，ODE45 求解器參數(shù)設(shè)置為相對誤差 εr=1×10-3，絕對誤差 εa=1×10-6，來源于外力邊界條件的受力端節(jié)點(diǎn)應(yīng)變約束為：εs=0,κs=0,κt=0,κb=0.

圖14 旋轉(zhuǎn)柔性梁初始狀態(tài)及動(dòng)態(tài)載荷Fig.14 The initial state and dynamic loads of the rotating flexible beam

圖15給出了利用本文算法計(jì)算的結(jié)果與文獻(xiàn)[35]對比的結(jié)果曲線.

從圖15可以看出，旋轉(zhuǎn)柔性梁端部的水平位移、豎直位移以及垂直紙面方向的位移呈周期性變化，且由于繞z軸轉(zhuǎn)角沒有限制的原因，水平位移幅值較大.本文內(nèi)部含三個(gè)節(jié)點(diǎn)的一個(gè)樣條單元（即內(nèi)部有4個(gè)小單元）的仿真結(jié)果與文獻(xiàn)[35]中20個(gè)單元仿真結(jié)果吻合得比較好.圖16給出了15 s 和30 s 時(shí)梁截面的廣義應(yīng)變隨弧長的變化曲線，數(shù)據(jù)顯示梁的廣義應(yīng)變沿形心線是光滑連續(xù)的，依據(jù)本構(gòu)關(guān)系可知廣義應(yīng)力同樣沿形心線是連續(xù)的.

圖15 旋轉(zhuǎn)梁自由端的位移-時(shí)間曲線Fig.15 Displacement-time curves at the free end of the rotating beam

圖16 時(shí)刻15 s 和30 s 時(shí)，旋轉(zhuǎn)梁的廣義應(yīng)變Fig.16 Generalized strains of the rotating beam at 15 s and 30 s

9.4 算例4：彎頭處受剪力作用的“L”形空間懸臂梁

如圖17所示，“L”形懸臂梁的彎頭處受指向面外的動(dòng)態(tài)載荷Fz(t)的作用.載荷大小遵循圖中所示的函數(shù)關(guān)系.結(jié)構(gòu)參數(shù)采用無量綱記法，梁的長度L=10，線密度 ρA=1，抗拉模量EA=104，抗彎抗扭模量GJ=EIy=EIz=103，截面初始構(gòu)型慣性矩Jp=diag(20,10,10).光滑因子取0.001，ODE45 求解器參數(shù)設(shè)置為相對誤差εr=1×10-3，絕對誤差εa=1×10-6.來源于外力邊界條件的自由端節(jié)點(diǎn)應(yīng)變約束為：εs=0,κs=0,κt=0,κb=0.

圖17 “L”形懸臂梁初始狀態(tài)及動(dòng)態(tài)載荷Fig.17 The initial state and dynamic loads of the L-shaped cantilever beam

圖18給出了利用本文算法計(jì)算的結(jié)果與文獻(xiàn)[36]對比的結(jié)果曲線.

算例4 來源于文獻(xiàn)[36]，本文將“L”形懸臂梁的每個(gè)分支都用相同的樣條梁單元表示，樣條內(nèi)部為一個(gè)節(jié)點(diǎn)(即內(nèi)部含2個(gè)小單元，“L”形結(jié)構(gòu)共4個(gè)單元)進(jìn)行仿真.從圖18可以看出，仿真結(jié)果中自由端與彎頭處的z向位移與文獻(xiàn)中采用20個(gè)單元的仿真結(jié)果相符合.“L”形懸臂梁自由端不受任何載荷作用，理論上此處廣義應(yīng)力均為零，圖19中給出了自由端節(jié)點(diǎn)處軸線應(yīng)變和截面曲率隨時(shí)間的變化曲線，數(shù)據(jù)顯示自由端節(jié)點(diǎn)處的軸向應(yīng)變和截面曲率很小，在計(jì)算誤差范圍內(nèi)，說明外力邊界條件約束的準(zhǔn)確性.

圖18 “L”形懸臂梁自由端與彎頭處的位移曲線Fig.18 The displacement curves at the free end and the elbow of the L-shaped cantilever beam

圖19 “L”形懸臂梁自由端節(jié)點(diǎn)的廣義應(yīng)變Fig.19 Generalized strains at the free end of the L-shaped cantilever beam

9.5 算例5：端部處受組合力作用空間梁結(jié)構(gòu)

如圖20所示，三根單位長度梁垂直連接組成的空間結(jié)構(gòu)的自由端受動(dòng)態(tài)載荷F的作用.載荷大小遵循圖中所示的函數(shù)關(guān)系，結(jié)構(gòu)參數(shù)采用無量綱記法，梁的長度均為L=2，截面形式為0.1 × 0.1的正方形，彈性模量E=1×109Pa，Poisson 比 ν=0.3，密度 ρ=1 000 kg·m-3，光滑因子取0.001，ODE45 求解器參數(shù)設(shè)置為相對誤差εr=1×10-3，絕對誤差εa=1×10-6.來源于外力邊界條件的自由端節(jié)點(diǎn)應(yīng)變約束為：κs=0,κt=0,κb=0.

圖20 空間組合梁結(jié)構(gòu)Fig.20 A spatial structure composed of 3 beams

本文將空間組合梁結(jié)構(gòu)的每個(gè)分支都用相同的樣條梁單元表示，樣條內(nèi)部為3個(gè)節(jié)點(diǎn)(即內(nèi)部含4個(gè)小單元，結(jié)構(gòu)共12個(gè)單元)進(jìn)行仿真.方便對比，采用ANSYS Workbench 對同樣的結(jié)構(gòu)進(jìn)行動(dòng)力學(xué)分析，提取受力端節(jié)點(diǎn)位移；得到本文計(jì)算結(jié)果與ANSYS 結(jié)果的對比曲線，如圖21所示.

由圖21可知，梁結(jié)構(gòu)端部位移已經(jīng)超過單個(gè)梁的長度尺寸大小，整個(gè)結(jié)構(gòu)發(fā)生梁截面大位移大轉(zhuǎn)動(dòng)的變形，對比結(jié)果顯示本文單元針對復(fù)雜空間梁結(jié)構(gòu)計(jì)算的有效性.

圖21 空間組合梁結(jié)構(gòu)受力端位移Fig.21 Displacements at the end of the spatial structure composed of beams

上述幾個(gè)算例表明，在較大的載荷變化范圍內(nèi)，本文所得的結(jié)果與解析解、商業(yè)軟件結(jié)果吻合得很好；相較傳統(tǒng)梁單元離散方式，本文所提單元避免了位移和轉(zhuǎn)角獨(dú)立插值帶來的剪切閉鎖問題，采用樣條插值，節(jié)點(diǎn)處導(dǎo)數(shù)連續(xù)性高、光滑性好，可以有效地求解大變形梁結(jié)構(gòu)，保證了各分支梁內(nèi)部節(jié)點(diǎn)處廣義應(yīng)力連續(xù)，濾掉剛性較強(qiáng)的梁結(jié)構(gòu)中高頻成分以在滿足精度的要求下提高計(jì)算效率.此外，樣條單元內(nèi)部每一個(gè)節(jié)點(diǎn)參數(shù)為5個(gè)，要比傳統(tǒng)的節(jié)點(diǎn)參數(shù)數(shù)量少，節(jié)點(diǎn)參數(shù)數(shù)目的優(yōu)勢會(huì)隨著單元的增多而愈發(fā)明顯.

10 結(jié) 論

實(shí)際工程梁結(jié)構(gòu)中，諸如大型起重機(jī)伸縮臂結(jié)構(gòu)、桁架式臂架結(jié)構(gòu)等，它們的梁分支較多，大部分屬于細(xì)長梁，往往伴隨著大變形現(xiàn)象的發(fā)生，單元節(jié)點(diǎn)數(shù)可能成百上千.本文在大變形梁虛功率原理的基礎(chǔ)上提出了一種節(jié)點(diǎn)參數(shù)含應(yīng)變的幾何非線性空間樣條梁單元，它可以滿足工程要求精度的同時(shí)，降低計(jì)算成本，具有以下幾個(gè)特點(diǎn)：

1）傳統(tǒng)幾何非線性單元分析細(xì)長梁時(shí)，剪切閉鎖造成的誤差無法滿足工程要求，本文考慮Euler-Bernoulli 梁變形耦合關(guān)系，建立梁單元轉(zhuǎn)動(dòng)參數(shù)和位移參數(shù)間的函數(shù)方程，可準(zhǔn)確地仿真細(xì)長梁.

2）具有與截面剛體運(yùn)動(dòng)無關(guān)的廣義應(yīng)變（弧長伸縮率、截面曲率），可用于描述單元水平上的幾何非線性問題；保證單元間應(yīng)力連續(xù)情況下，采用樣條插值進(jìn)行自由度縮減式組裝單元，將廣義應(yīng)變納入到邊界節(jié)點(diǎn)參數(shù)當(dāng)中，便于精確地施加力的邊界條件.

3）大變形梁的幾何非線性效應(yīng)以及剛性成分會(huì)使得計(jì)算效率十分低下，特別是在機(jī)械領(lǐng)域，往往關(guān)注的焦點(diǎn)是系統(tǒng)的宏觀動(dòng)力學(xué)行為，鑒于此，本文將梁結(jié)構(gòu)系統(tǒng)方程進(jìn)行降噪處理，濾掉其解中的高頻成分，工程設(shè)計(jì)人員可根據(jù)精度要求的需要，適當(dāng)調(diào)整光滑因子，保證要求精度的情況下，實(shí)現(xiàn)柔性梁的快速求解，降低求解空間幾何非線性梁單元系統(tǒng)的計(jì)算成本.

本文的方法與經(jīng)典算例的結(jié)果吻合度較高.在單元數(shù)目上少于文獻(xiàn)中采用的單元數(shù)，而且中間節(jié)點(diǎn)參數(shù)數(shù)量較少，有利于縮減較多節(jié)點(diǎn)數(shù)目的工程梁結(jié)構(gòu)的自由度規(guī)模.