節(jié)點參數(shù)含應(yīng)變的空間幾何非線性樣條梁單元*

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

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

引 言

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

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

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

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

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

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

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

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

形心線上的任意點矢徑可以用五次多項式插值擬合得到

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

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

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

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

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

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

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

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

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

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

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

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

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

以及

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

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

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

其中，弧長伸縮率

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

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

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

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

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

其中

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

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

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

因此，式(26)可改寫為

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

其中

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

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

其中

因此，式(33)可改寫為

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

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

它的形函數(shù)為

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

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

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

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

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

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

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

由剛體動力學(xué)方程可得微元體的運動方程為

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

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

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

其中力和力矩分量來源于

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

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

它的標準形式通?？筛膶憺?/p>

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

4 梁單元的角速度和曲率

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

其中

因此，它在截面坐標下的分量為

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

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

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

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

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

其中

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

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

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

其中

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

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

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

單元的慣性力虛功率為

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

變形虛功率為

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

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

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

6 梁單元間的連接光滑化

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

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

其中，單元長度

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

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

其中，矩陣

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

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

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

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

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

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

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

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

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

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

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

式中

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

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

9 數(shù)值算例

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

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

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

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

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

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

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

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

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

表1 效率比較（柔性單擺機構(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é)點的應(yīng)變約束為

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

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

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

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

表2 懸臂梁受集中力矩 M=2πL-1EI 時選取不同光滑因子仿真用時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 時，梁的端部位移與不采用降噪方法（Δt=0）時計算的結(jié)果相吻合，最大絕對誤差在0.1 以內(nèi)，但耗時從3 653 s 縮短至982 s，若對精度要求不高，繼續(xù)增大光滑因子，耗時將進一步縮短.圖11、12 給出了在不同彎矩作用下梁最終平衡時的變形狀態(tài)及其端部位移與解析解的對比，文獻[33-34]已經(jīng)指出形心位移和轉(zhuǎn)角獨立插值的傳統(tǒng)梁，解決諸如此類發(fā)生大彎曲變形細長梁時，會出現(xiàn)剪切閉鎖造成誤差增大的現(xiàn)象，對比曲線顯示本文單元計算結(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)時，小變形范圍內(nèi)端部位移存在解析解uy=FL3(3EIz)-1，大變形結(jié)果可采用ANSYS 軟件打開“NLGEOM”命令給出合理的解.表3和圖13給出了采用本文單元建模，不同剪力作用下梁最終平衡狀態(tài)下端部豎向位移，以及和解析解、ANSYS 結(jié)果的對比.由表3和圖13可知，剪力大小在1～3 N 范圍內(nèi)時，梁的變形屬于小變形，與解析解的誤差在1%以內(nèi)；隨著剪力的增大，梁的變形逐漸進入大變形狀態(tài)，在36 N 時變形位移已經(jīng)超出梁長的1/2，此階段，解析解已無法滿足要求，本文計算結(jié)果與ANSYS 大變形計算結(jié)果誤差隨著剪力增大略微增大，但均在0.13%以內(nèi).

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

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

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

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

大范圍運動的旋轉(zhuǎn)柔性梁已被很多學(xué)者用來考察所提方法的有效性.如圖14所示，梁的左端可以繞z軸自由旋轉(zhuǎn)，右端施加y和z方向上兩個隨時間變化的動態(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é)點應(yīng)變約束為：εs=0,κs=0,κt=0,κb=0.

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

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

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

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

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

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

如圖17所示，“L”形懸臂梁的彎頭處受指向面外的動態(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é)點應(yīng)變約束為：εs=0,κs=0,κt=0,κb=0.

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

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

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

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

圖19 “L”形懸臂梁自由端節(jié)點的廣義應(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)的自由端受動態(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é)點應(yīng)變約束為：κs=0,κt=0,κb=0.

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

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

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

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

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

10 結(jié) 論

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

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

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

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

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