節點參數含應變的空間幾何非線性樣條梁單元*

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

（1.大連理工大學 工業裝備結構分析國家重點實驗室，遼寧 大連 116023；2.大連理工大學 海洋科學與技術學院，遼寧 盤錦 124221）

引 言

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

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

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

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

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

1 梁單元位移形函數

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

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

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

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

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

對式(1)求時間導數，可得

對式(2)求時間導數，可得

對式(3)求時間導數，可得

2 梁單元的轉動描述

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

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

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

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

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

式中，分解系數為

對式(14)求弧長導數和時間導數，可得

以及

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

曲率κ的時間導數可用角速度分量的弧長導數表示：

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

其中，弧長伸縮率

對式(19)求時間的導數，可得

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

因此，前兩次轉動的Cardan 角為

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

對式(22)求時間導數，可得

其中

它們與es相互正交，基于此關系，容易得到Cardan 角對時間的導數為

對式(19)、(20)求時間導數，可得

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

因此，式(26)可改寫為

矢量h,b的時間導數可表示為

其中

同理，Cardan 角的弧長導數為

它們的時間導數為

其中

因此，式(33)可改寫為

對式(29)求時間導數，可得

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

它的形函數為

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

對式(38)、(39)求時間導數，可得

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

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

3 梁單元的廣義應變

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

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

由剛體動力學方程可得微元體的運動方程為

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

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

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

其中力和力矩分量來源于

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

相應的本構關系為

它的標準形式通常可改寫為

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

4 梁單元的角速度和曲率

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

其中

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

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

對式(56)、(57)求時間導數，可得

其中，轉換矩陣Tωκ的時間導數為

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

據此，式(59)可改寫為

其中

定義節點參數矢量組成的列向量：

從上述的分析可以看出，θ是節點參數q的函數，因此，它的時間導數可表示為

其中，轉換矩陣Tθq，Tθ′q僅是q的函數，向量aθ，aθ′與加速度項無關.據此，角速度和曲率分量以及它們的時間導數可寫為

其中

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

5 梁單元的質量陣和力陣

根據上述的討論，單元內任意點矢徑和應變的時間導數可以表示為以下的矩陣形式：

單元的慣性力虛功率為

其中，質量陣、力陣為

變形虛功率為

其中，彈性節點力為

經過單元組裝，可得梁系統的虛功率方程為

其中，廣義位移u是所有節點參數組成的列向量，fin,fa分別為廣義內力和廣義外力.為了進一步獲得系統動力學方程，往往需要繼續分析參數u的獨立性.

6 梁單元間的連接光滑化

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

圖4 縮減節點參數后的樣條單元Fig.4 Parameter reduction of spline elements

其中，單元長度

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

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

其中，矩陣

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

這樣，可以采用樣條單元直接建立梁結構的整體方程，內部節點處矢徑對弧長的導數可以不作為節點參數，縮減系統的自由度，并且隨著內部單元的增多，優勢越明顯.

7 樣條梁單元邊界節點參數轉換

通過式(57)，Cardan 角的弧長導數可以表示為

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

圖5 樣條單元的邊界節點參數Fig.5 Boundary nodal parameters in spline elements

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

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

在t+τ時刻，單元應力可以近似表示為

它在極短時間區間(t,t+Δt)內的平均值為

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

式中

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

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

9 數值算例

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

9.1 算例 1：單擺機構自由下落

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

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

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

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

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

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

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

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

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

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

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

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

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

這里將梁作為一個樣條單元，內部取3個節點進行計算，ODE45 求解器參數設置為相對誤差εr=1×10-3，絕對誤差 εa=1×10-6.當端部受彎矩作用，懸臂梁發生純彎曲，存在解析解 κ=M/(EIz)，即發生純彎曲時彎曲曲率κ 與外力矩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

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

圖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

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

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

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

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

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

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

圖14 旋轉柔性梁初始狀態及動態載荷Fig.14 The initial state and dynamic loads of the rotating flexible beam

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

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

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

圖16 時刻15 s 和30 s 時，旋轉梁的廣義應變Fig.16 Generalized strains of the rotating beam at 15 s and 30 s

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

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

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

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

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

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

圖19 “L”形懸臂梁自由端節點的廣義應變Fig.19 Generalized strains at the free end of the L-shaped cantilever beam

9.5 算例5：端部處受組合力作用空間梁結構

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

圖20 空間組合梁結構Fig.20 A spatial structure composed of 3 beams

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

由圖21可知，梁結構端部位移已經超過單個梁的長度尺寸大小，整個結構發生梁截面大位移大轉動的變形，對比結果顯示本文單元針對復雜空間梁結構計算的有效性.

圖21 空間組合梁結構受力端位移Fig.21 Displacements at the end of the spatial structure composed of beams

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

10 結 論

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

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

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

3）大變形梁的幾何非線性效應以及剛性成分會使得計算效率十分低下，特別是在機械領域，往往關注的焦點是系統的宏觀動力學行為，鑒于此，本文將梁結構系統方程進行降噪處理，濾掉其解中的高頻成分，工程設計人員可根據精度要求的需要，適當調整光滑因子，保證要求精度的情況下，實現柔性梁的快速求解，降低求解空間幾何非線性梁單元系統的計算成本.

本文的方法與經典算例的結果吻合度較高.在單元數目上少于文獻中采用的單元數，而且中間節點參數數量較少，有利于縮減較多節點數目的工程梁結構的自由度規模.