ANCF/CRBF 平面梁閉鎖問題及閉鎖緩解研究1)

張大羽 羅建軍 王輝 馬小飛

?(中國空間技術研究院西安分院,西安 710000)

?(西北工業大學航天學院,航天飛行動力學技術重點實驗室,西安 710072)

引言

柔性梁被廣泛地應用于航天與機械系統,是柔性多體系統動力學最典型的一類研究對象,建立其動力學模型可為復雜柔性多體系統設計和分析提供重要的理論依據.建立柔性梁模型的關鍵是如何準確地描述梁單元連續的旋轉場,特別是當柔性梁經歷大轉動大變形運動.區分不同類型的非線性柔性梁建模方法的一個重要特征是如何準確地描述系統中的有限轉動[1-10].

共旋坐標法[11-15]在處理有限轉動時通過非線性迭代求解無限小的轉動矢量,并將計算的轉動增量疊加到方向余弦矩陣上.但是這類方法在求解時必須確保迭代步長不能過大,這樣才能將轉動增量近似為小轉動的定軸旋轉.此外,共旋坐標法將全局位移場分解為剛體運動和彈性變形,使得單元的慣性力等首先定義在遷移坐標系內,而后再轉換到全局坐標系,同樣要求了該方法的迭代步長不能過大.但是,過小的迭代步長會造成系統運動方程的線性化,導致柔性體在其純剛體運動下產生非零的應變能,因此共旋坐標法僅適合用于模擬具有小應變的系統.幾何精確梁方法[16-19]采用慣性坐標系下的節點絕對位移來描述結構的整體運動,不再對剛體運動和彈性運動進行區分,避免了增量法中的線性化問題,因此可以處理柔性梁的大變形問題.幾何精確梁方法采用絕對位置向量和有限轉角參數作為節點坐標,位移場和旋轉場采用兩個獨立的插值函數,其中對有限轉動插值方法的研究是該方法研究的一個重要分支.對于幾何精確梁獨立的旋轉場插值策略,Shabana 和Ding[20]指出該策略會導致坐標的冗余及梁內同一物質點在位移場和轉動場上的幾何構型不一致,并通過算例證明了采用兩個獨立的插值函數表述位移和轉動會導致梁在同一時刻存在兩種不同的構型,與剛體動力學假設不相符.相比獨立的旋轉場插值,Shabana[21]提出了絕對節點坐標法(absolute nodal coordinate formulation,ANCF),該方法采用統一的插值函數描述梁內物質點的位置和轉動.通過慣性坐標系下的梯度向量描述柔性體的方向,避免了對轉角參數的相加運算,致使求解過程為非增量形式,這點區別于幾何精確梁方法.此外,該方法還具有常數質量陣及不存在科氏力和離心力項的優點,適用于解決大轉動大變形的柔性體的建模問題.近年來,絕對節點坐標方法已受到越來越多的學者的關注與研究[22-30,53-55].

絕對節點坐標法采用梯度向量描述柔性體轉動和應變,雖然可以描述系統運動過程中的大轉動大變形,但也存在一些缺陷.首先,相比轉角參數,梯度向量不具有明確的轉動物理意義,不利于面向單元旋轉的設計和控制.例如,在變形過程中梯度變化方向不代表實際的系統旋轉方向,對梯度向量的控制最終未必得到想要控制的轉動量.其次,相比轉角,梯度向量的數量更多,導致系統自由度的升高,以至加重系統計算負擔.2015 年Shabana[31]提出基于一致旋轉場列式的絕對節點坐標單元(ANCF consistent rotationbased formulation,ANCF/CRBF).該單元在統一的位移場和旋轉場描述基礎上,通過速度轉換矩陣將梯度向量由轉動參數替代,成功引入了轉角坐標,解決了上述缺陷.此外,ANCF/CRBF 單元采用統一的插值函數解決了位移場和旋轉場分開插值策略而導致坐標和幾何構型冗余的問題.例如,線性插值的位移場得到零曲率,而線性插值的轉動場得到高度非線性曲率,導致單元的曲率定義不統一.隨后,Zheng 等[32]建立了ANCF/CRBF 平面梁單元,并進行了單元動力學分析.

但是,已有ANCF/CRBF 梁在ANCF 全參數平面梁的框架上開發的,同樣存在嚴重的閉鎖問題,特別是泊松閉鎖,會給數值計算帶來困難.ANCF 單元泊松閉鎖問題產生的原因是單元的軸向高階插值和橫向低階插值不一致導致不同方向間的應力相互耦合,致使單元在彎曲過程中由軸向的彎曲變形高階項產生出偽應力,導致單元性能下降,呈現出“彎曲鎖死”的現象.對于傳統ANCF 全參數梁單元的泊松閉鎖問題及其緩解技術,典型的研究有:Gerstmayr等[33]采用縮減積分方法(selectively reduced integration method,SRIM)重新構造了平面梁單元的彈性力矩陣,選擇性地避開了對彎曲變形的高階項的積分,可以有效緩解泊松閉鎖.Mikkola 等[34]在平面ANCF梁單元基礎上采用對橫向變形場獨立插值技術,提出了一種基于混合插值的平面ANCF 梁單元.動力學仿真結果表明該單元具有更高的彎曲變形精度和收斂性.Schwab 和Meijaard[35]采用彈性線方法(elastic center line,ECL)構造了新的彈性力公式,消除了軸向與橫向變形間的高階耦合項,緩解了泊松閉鎖問題.此外,Schwab 和Meijaard[35]分別采用基于Hellinger-Reissner 和Veubeke-Hu-Washizu 變分原理的方法用于消除三維ANCF 全參數梁單元剪切閉鎖問題,改善了單元的彎曲性能.Hussein 等[41]采用Hellinger-Reissner 方法消除了二維ANCF 全參數梁單元的剪切閉鎖問題,并分別測試了Hellinger-Reissner 方法在薄梁和厚梁動力學分析中消除閉鎖的效果.李文雄等[42]提出基于Hellinger-Reissner 變分原理的無閉鎖曲梁單元.與常規曲梁單元相比,該單元具有更高的精度和效率.Patel 和Shabana[36]提出了的一種全新的緩解單元泊松閉鎖的方法,應變分解方法(strain split method,SSM).其核心是將應變分解為單元軸向應變和截面應變,這樣可將由泊松效應引起的軸向彎曲變形高階項與截面變形項區分出來,然后分別施加相應的彈性矩陣,將泊松效應僅計入到軸向上,進而起到泊松閉鎖緩解作用.縮減積分方法和彈性線方法都是從單元積分角度選擇性地避開計算應變?應力關系的高階耦合項,這兩種方法適用于緩解單元的泊松和剪切閉鎖;基于多場變分原理是從單元位移場角度引入單獨的剪切場插值函數,從而消除橫截面剪切變形與軸向拉伸和彎曲變形的耦合項,該方法主要用于解決單元的剪切閉鎖問題;應變分解方法是從單元運動學角度將高階耦合項解耦出來,從而僅將泊松效應施加在應變?應力關系的低階線性項上,該方法用于解決單元的泊松閉鎖問題.對于ANCF/CRBF單元的泊松閉鎖問題鮮有研究,對該單元閉鎖問題的深入探討是十分必要的,有助于提升單元的計算精度和確定單元的應用范圍.

本文系統研究了ANCF/CRBF 單元特有的泊松閉鎖問題及其緩解技術.ANCF/CRBF 梁單元的推導是通過對ANCF 梁單元的節點處橫截面施加運動學約束得到的.由于約束后的橫截面不能變形,消除了節點處的泊松閉鎖,但是非節點處仍然存在泊松閉鎖.這將導致ANCF/CRBF 梁單元與ANCF 梁單元不同的泊松閉鎖特性,同時也使得泊松閉鎖成為ANCF/CRBF 單元主要的閉鎖問題.首先,為了全面地理解ANCF/CRBF 梁單元特有的泊松閉鎖特性以及深入探索施加運動學約束對ANCF/CRBF 梁單元泊松閉鎖特性的影響,在已有ANCF/CRBF 剛性截面梁的基礎上,進一步對單元的軸向梯度進行約束,建立了兩種新的剛性截面的ANCF/CRBF 梁單元.新單元在節點處完全消除了泊松閉鎖和剪切效應,這是與傳統ANCF/CRBF 剛性截面梁單元的不同之處.然后,采用兩種閉鎖緩解技術,彈性線方法和應變分解方法,提升了ANCF/CRBF 單元的收斂性和計算精度.以期為ANCF/CRBF 這類單元的閉鎖問題分析以及單元的應用范圍提供參考數據.

1 傳統的ANCF 全參數平面梁單元

傳統的ANCF 全參數平面梁單元是絕對節點坐標方法中的經典單元,其采用統一的插值函數描述單元的平動和轉動.該單元在運動過程中橫截面與梁軸線不保持時刻垂直,截面可變形.

1.1 單元運動學描述

ANCF12 單元采用統一的插值策略,這會導致梁的軸向彎曲變形高階項與橫向低階變形高度耦合,產生嚴重的泊松閉鎖問題,致使計算精度下降.此外,ANCF12 單元采用梯度坐標描述單元的轉動信息,但是它的物理意義不如轉角參數明確.

1.2 單元彈性力公式

圖1 三種變形構型之間的關系Fig.1 Three general configurations

梁單元i的Green-Lagrange 應變張量εi可通過變形梯度張量Ji得到

其中,λ=Eν/[(1+ν)(1 ?2ν)],μ=E/[2(1+ν)];E是彈性模量,ν 是泊松比.

1.3 單元動力學方程

2 ANCF/CRBF 平面梁單元

傳統ANCF 全參數平面梁單元采用梯度坐標描述單元的轉動信息,但是其物理意義不如轉角參數明確,并且采用梯度坐標會使單元的廣義坐標數目增多,導致動力學計算效率降低.基于一致旋轉場列式的ANCF/CRBF 單元的提出可以解決以上問題.下面分別介紹已有的ANCF/CRBF 平面梁單元,以及本文新開發的兩種剛性截面ANCF/CRBF 平面梁單元.

2.1 傳統的ANCF/CRBF 剛性截面梁單元

2.2 ANCF/CRBF 軸向不可伸展的剛性截面梁單元

對于梁單元i的轉換關系如下[43]

2.3 ANCF/CRBF 軸向可伸展的剛性截面梁單元

對于梁單元i的轉換關系如下[43]

2.4 單元動力學方程

值得注意的是,ANCF/CRBF 單元計算轉角是通過積分角度對時間的導數得到的,而非積分角速度.分空間和平面兩種情況進行討論.

(1)在空間情況下,以全參數ANCF 梁單元與ANCF/CRBF 梁單元轉換為例,其速度層面的轉換關系為[31,44-45]

文中還指出ANCF/FFR 單元雖然可準確描述初始彎曲構型,有利于開發一種基于力學的計算機輔助設計和分析(I-CAD-A)融合方法.但是,由于引入了小轉動坐標會導致這類單元不能精確地表示剛體慣性,會影響單元的動力學精度.因此對ANCF/FFR單元的閉鎖緩解技術研究可作為下一步研究的方向.

3 單元閉鎖緩解技術

ANCF 全參數平面梁單元存在著嚴重的泊松閉鎖問題,會嚴重影響單元的計算精度.針對如何緩解ANCF 單元泊松閉鎖問題,已有大量地研究,其中典型的方法有彈性線方法和應變分解方法.本節對兩種方法進行逐一介紹.

3.1 彈性線方法

由彈性線方法(ECL)可知[21,35],應變勢能通過對沿中心軸線xi積分得到,由軸向拉伸勢能剪切勢能和橫向彎曲勢能組成

其中,E,G,A,I分別是彈性模量、剪切模量、橫截面面積和截面慣性矩.剪切模量表示為G=E/[2(1+ν)].As=ks·A,其中,ks是剪切修正因子.對于橫截面是矩形的情況,ks=10(1+ν)/(12+11ν).對于軸向應變、剪切應變、及曲率公式κi分別可表示為

在彈性線方法中,由于僅對梁軸向中心線進行積分,一定程度上可以緩解由高階彎曲變形和低階橫向變形耦合引起的閉鎖問題.彈性線方法適合用于計算較薄的直梁結構.

3.2 應變分解方法

由應變分解方法(SSM)可知[36],對于梁單元i,其位置和梯度向量可分解為

依據Voigt 標記規則,可寫為

對于ANCF 全參數平面梁單元,從式(20)分析得知,其彎曲變形的高階項在單元非軸向分應變εiS中,會產生偽應力,這是導致單元泊松閉鎖的原因.因此,在求解第二類Piola-Kirchhoff 應力時對兩個Green-Lagrange 分應變εiC,εiS分別施加不同的彈性矩陣,可表示為σ=EiCεiC+EiSεiS,其中彈性矩陣分別為

由式(22)可以看出,應變分解方法將泊松效應僅施加在低階軸向分應變εiC上,可有效地將梁橫截面低階正應變與由彎曲變形高階項引起的偽應變分離開來,使得泊松效應充分地體現在低階正應變中,有效地緩解了ANCF 梁單元泊松閉鎖問題.在應變分解方法中,梁單元的彈性力向量可計算為

4 算例分析

本文設計5 種算例用于全面測試新單元的性能,分別是小變形細長梁靜力學算例、大變形深梁靜力學算例、懸臂梁末端受彎矩靜力學算例、懸臂梁受重力動力學算例和柔性單擺動力學算例.此外,本文還對兩種閉鎖緩解技術應用在新單元上的效果進行研究.為了方便對比分析,仿真中的模型可歸納為10種類型,分別是:

(1)ANCF12-GCM (傳統ANCF 全參數梁單元+連續介質力學);

(2)ANCF/CRBF10-GCM (傳統剛性截面ANCF/CRBF 梁單元+連續介質力學);

(3)ANCF/CRBF10-ECL (傳統剛性截面ANCF/CRBF 梁單元+彈性線方法);

(4)ANCF/CRBF10-SSM (傳統剛性截面ANCF/CRBF 梁單元+應變分解方法);

(5)ANCF/CRBF8-GCM(軸向可伸展的剛性截面ANCF/CRBF 梁單元+連續介質力學);

(6)ANCF/CRBF8-ECL (軸向可伸展的剛性截面ANCF/CRBF 梁單元+彈性線方法);

(7)ANCF/CRBF8-SSM (軸向可伸展的剛性截面ANCF/CRBF 梁單元+應變分解方法);

(8)ANCF/CRBF6-GCM(軸向不可伸展的剛性截面ANCF/CRBF 梁單元+連續介質力學);

(9)ANCF/CRBF6-ECL (軸向不可伸展的剛性截面ANCF/CRBF 梁單元+彈性線方法);

(10)ANCF/CRBF6-SSM(軸向不可伸展的剛性截面ANCF/CRBF 梁單元+應變分解方法).

根據第2 節和第3 節對絕對節點坐標全參數梁和基于一致旋轉場的絕對節點坐標梁單元的描述,詳細分析各類型單元泊松閉鎖的特點.

(1)ANCF12-GCM

①軸向可伸展、截面可變形;

②旋轉描述中無轉角參數,僅有梯度向量;

③節點處存在剪切變形;

④截面無約束,節點處和單元內均有泊松閉鎖.

(2)ANCF/CRBF10-GCM

①軸向可伸展、截面不可變形;

②旋轉描述中既有轉角參數又有梯度向量;

③節點處存在剪切變形;

④節點截面全約束,節點處無泊松閉鎖,單元內存在泊松閉鎖.這種不連續情況導致單元整體的泊松閉鎖問題相較ANCF12 單元更加嚴重.

(3)ANCF/CRBF8-GCM

①軸向可伸展、截面不可變形;

②旋轉描述中僅有轉角參數;

③節點處不存在剪切變形;

④節點截面全約束,節點處無泊松閉鎖,單元內存在泊松閉鎖.此外,節點處無剪切變形,單元內存在剪切變形.這種不連續情況導致單元整體的泊松閉鎖相較ANCF10 單元更加嚴重.

(4)ANCF/CRBF6-GCM

①軸向不可伸展、截面不可變形;

②旋轉描述中僅有轉角參數;

③節點處不存在剪切變形;

④節點截面全約束,節點處無泊松閉鎖,單元內存在泊松閉鎖.此外,單元軸線不能伸展.這會導致單元整體的泊松閉鎖相較ANCF8 單元更加嚴重.

4.1 靜力學測試:小變形細長梁算例

該算例為典型的懸臂梁受集中載荷問題,考察ANCF/CRBF 單元在小變形靜力學分析中的準確性.圖2 所示為懸臂梁.梁的參數分別為長度l=1 m,高度和厚度h=w=0.01 m,彈性模量E=200 GPa,泊松比ν=0.3,以及末端集中力設為F=?10 N.采用Newton-Raphson 迭代求解其非線性方程組,相對誤差設置為RelT=10?6.Y方向位移具有理論解為δ=?0.02 m.

表1 為不同單元類型采用連續介質力學方法時計算的Y向位移.首先,發現ANCF12 單元、ANCF/CRBF10 單元、ANCF/CRBF6 和ANCF/CRBF8 單元均不能收斂到理論解,這是因為采用連續介質力學方法計算的彈性力會造成軸向和截面的正應變高度耦合,使得泊松效應更加地明顯.其次,ANCF 全參數平面梁單元和ANCF/CRBF 單元在小變形靜力學問題上表現相當.表2 為不同單元類型采用閉鎖緩解方法時計算的Y向位移.以ANCF/CRBF10 單元作為參考解,可以發現ANCF/CRBF6 和ANCF/CRBF8 單元的計算精度低于參考解.這是因為這兩個單元都對節點處梯度向量進行了正交約束,可捕捉到的截面變形類型減少,所以計算精度不如傳統ANCF 全參數單元.此外,約束截面雖然可以消除節點處的泊松閉鎖,但單元內的截面可以變形,仍然存在閉鎖問題.圖3 所示為ANCF/CRBF8 單元采用三種彈性力計算方式下的Y向位移誤差收斂率,可以發現在小變形情況下,采用彈性線方法或應變分解法時的收斂率相當.由于本算例的特色為細長梁,截面變形小,泊松閉鎖問題不明顯.下一個算例將采用深梁模型進一步研究單元的閉鎖問題.

表1 小變形情況下Y 向位移對比(m)Table 1 Tip vertical displacement of beam in case of small deformation(m)

表2 小變形情況下Y 向位移對比(使用閉鎖緩解技術)(m)Table 2 Tip vertical displacement of beam in case of small deformation using locking alleviation methods(m)

圖3 ANCF/CRBF8 單元Y 向位移相對誤差Fig.3 Relative error of tip vertical displacement of ANCF/CRBF8 beam in case of small deformation

4.2 靜力學測試:大變形深梁算例

本算例采用深梁模型,如圖2 所示.梁的參數分別為長度l=2 m,高度h=0.5 m 和厚度w=0.1 m,彈性模量E=207 GPa,以及泊松比ν=0.3.由于深梁的橫截面有著較大的變形,可完全捕捉泊松效應,因此單元的泊松閉鎖問題更加嚴重,該算例常被用于測試單元受較大變形時的靜力學計算精度.為了驗證大變形算例,設置末端集中力為F=?5.0×108×h3N,與文獻[33]中的參數相同.采用Newton-Raphson 迭代求解其非線性方程組,相對誤差設置為RelT=10?6.本算例中的參考解為采用商業有限元軟件ANSYS 計算的數值收斂解δ=?0.713 420 m.

表3 為不同單元類型采用連續介質力學方法時計算的Y向位移.首先,以傳統ANCF12 單元作為參考,發現ANCF/CRBF 單元由于對梯度向量施加了約束,消除了部分的變形模態,因此其計算精度不如ANCF12 單元.其次,發現兩種剛性截面的ANCF/CRBF 單元的計算精度不如ANCF/ANCF10 單元.對于深梁結構,由于長寬比較小,梁內部的泊松閉鎖更加嚴重,導致ANCF/CRBF 單元的靜力學計算精度更差.另外,在ANCF/CRBF 單元中,對梯度向量約束的越多,靜力學計算精度越差.因此,相比ANCF 全參數梁單元,ANCF/CRBF 單元不具有靜力學計算優勢,特別是對于大變形深梁結構.表4 為不同單元類型采用閉鎖緩解方法時計算的Y向位移.對比發現僅有應變分解法計算的結果可收斂到數值解.圖4 所示為ANCF/CRBF8 單元的Y向位移誤差收斂率,發現在大變形情況下,應變分解法的收斂率最好.對于大變形深梁模型,截面泊松效應明顯,應變分解法可將梁軸線方向正應變和截面方向正應變分解,可完全解除軸向高階彎曲變形和截面低階變形間的高度耦合關系,泊松效應僅計入在軸向正應變上.而彈性線方法僅可去除掉一些不重要的高階項,但不能將軸向和截面變形解耦,單元仍然存在一定程度上的泊松閉鎖,因此計算精度不如應變分解法的計算精度高.

表3 大變形情況下Y 向位移對比(m)Table 3 Tip vertical displacement of beam in case of large deformation(m)

表4 大變形情況下Y 向位移對比(使用閉鎖緩解技術)(m)Table 4 Tip vertical displacement of beam in case of large deformation using locking alleviation methods(m)

圖4 ANCF/CRBF8 單元Y 向位移相對誤差Fig.4 Relative error of tip vertical displacement of ANCF/CRBF8 beam in case of large deformation

4.3 靜力學測試:懸臂梁受端部集中力矩算例

本算例驗證ANCF/CRBF 梁單元在強幾何非線性靜力學分析中的有效性,如圖5 所示.梁的參數分別為長度l=10 m,高度h=0.1 m 和厚度w=0.1 m,彈性模量E=210 GPa,以及考慮泊松效應,泊松比設為ν=0.3,與參考文獻[18]相同.在梁的末端所施加的理論力矩為M=πEI/l.該算例具有理論解,當M=2πEI/l時,直梁最終被彎曲成一個圓.采用牛頓迭代法求解,相對誤差設置為10?6,載荷迭代步數設為20.本算例彈性力的計算采用應變分解方法.分別采用ANCF12,ANCF/CRBF10,ANCF/CRBF8 和ANCF/CRBF6 模擬懸臂梁.為了全面測試ANCF/CRBF 單元承受彎曲的情況,所加載的彎矩由小到大,并分為多組進行測試.

圖5 懸臂梁末端受集中力矩作用Fig.5 Cantilever beam subjected by a torque

對于該強幾何非線性靜力學算例,對比表5 和表6 發現:首先,ANCF 單元由于存在閉鎖問題,測試發現采用40 個單元劃分梁模型才使得懸臂梁被彎成整圓(所受力矩為ML/(πEI)=2 時).其次,對三種ANCF/CRBF 單元進行測試,采用與ANCF 單元彎成整圓相同的單元數,發現ANCF/CRBF 單元在受載小彎矩時的結果與理論解相近(ML/(πEI)=0.2 ～1.0).當彎矩增大,ANCF/CRBF 單元的結果開始與理論解呈現較大的偏差(ML/(πEI)=1.2 ～2.0).圖6 是ANCF 單元和ANCF/CRBF 單元在所受力矩為ML/(πEI)=2 時的構型圖.

表5 梁末端X 向位移對比(m)Table 5 Comparison of tip axial displacement using strain split locking alleviation methods(m)

表6 梁末端Y 向位移對比(m)Table 6 Comparison of tip vertical displacement using strain split locking alleviation methods(m)

圖6 懸臂梁受端部集中力矩下的構型圖(40 個單元)Fig.6 Cantilever beam configurations subjected by a torque(40 elements)

根據從文獻[45]可知,只有單元上兩個節點之間的相對轉動不超過30?時,得到的可認為是單位向量,從而ANCF/CRBF 單元才有效.可以通過以下的計算證明:例如讓單元的一個節點處的轉角為0?,另一個節點處的轉角從0?變化到360?,計算ξ=0.5 處的范數中的誤差,如圖7 所示.可以得到,在兩個節點相對轉角為30?時,誤差為0.034,與文獻[45]中的結論相符.此時ANCF/CRBF 單元的可近似認為單位向量,并且單元有效.但是當相對轉角繼續增大時,范數的誤差也隨之增大,這也解釋了ANCF/CRBF 這類單元在承受較大力矩時計算精度差的原因.

圖7 ry 范數的誤差Fig.7 Error of the norm of ry

4.4 動力學測試:懸臂梁受重力作用算例

本小節采用文獻[39]中的懸臂深梁受到重力作用的算例,該算例專門用于分析泊松閉鎖問題對單元動力學特性的影響,如圖8 所示.本算例彈性力的計算分別采用連續介質力學方法、彈性線方法和應變分解方法.懸臂梁的參數分別為長度l=1 m,高度和厚度h=0.1 m,w=0.07 m,彈性模量E=60 MPa,密度ρ=7200 kg/m3,以及泊松比ν=0.3.在重力作用下,懸臂梁在XY平面內往復運動,仿真時間為1 s.柔性梁被劃分為10 個單元.采用顯式積分方法,軸向和橫向的高斯積分點分別為5 個和3 個.

圖8 懸臂梁動力學模型Fig.8 Cantilever beam dynamic problem

圖9 為梁末端端點處X方向位移變化.可以看到當單元采用連續介質力學方法計算的彈性力時,單元泊松閉鎖問題嚴重.當采用兩種閉鎖緩解技術均可以緩解閉鎖問題.從靜力學分析和本小節的動力學分析可以得到,應變分解方法既適用于細長梁結構也適用于深梁結構,是最有效的方法,而彈性線方法更適用于模擬細長梁結構.

圖9 梁端點處X 方向位移變化Fig.9 Tip horizontal displacement of the cantilever beam

4.5 動力學測試:柔性單擺算例

本小節采用經典的柔性單擺算例分析ANCF/CRBF10、ANCF/CRBF8 和ANCF/CRBF6 三種單元的動力學特性,如圖10 所示.根據前面靜力學分析的結論,本算例彈性力的計算采用應變分解方法.單擺的參數分別為長度l=1 m,高度和厚度h=w=0.02 m,彈性模量E=20 MPa,密度ρ=7200 kg/m3,以及泊松比ν=0.3.在重力作用下,單擺在XY平面內運動,仿真時間為2 s.采用20 個ANCF/CRBF 單元模擬.采用商業有限元軟件LS-DYNA 中的Belytschko-Schwer beam(BS beam)單元作為對比[40],采用30 個BS beam單元模擬.另外,為最大程度上保證自編程序和有限元軟件的一致性,均采用顯式積分方法,這樣可以避免采用隱式積分算法時的數值阻尼對計算結果的影響.此外,自編程序和有限元軟件在截面方向的高斯積分點的選取均為3 個.

圖10 柔性單擺模型Fig.10 Flexible pendulum

圖11 梁端點處Y 方向位移變化Fig.11 Tip vertical displacement of the flexible pendulum

圖12 ANCF/CRBF8 單元能量變化Fig.12 Energy change of ANCF/CRBF8 element

圖13 梁中點軸向應變變化Fig.13 Axial strain change of beam middle point

圖11 為梁末端端點處Y方向位移變化,可以看出ANCF/CRBF10 和ANCF/CRBF8 單元計算的結果與LS-DYNA 軟件計算的結果有較好的吻合度,驗證了單元在大范圍剛體運動時可準確描述彈性變形的能力.而ANCF/CRBF6 單元由于軸向不能伸展,計算結果錯誤,不適用于模擬大變形梁問題.圖12 是ANCF/CRBF8 單元的能量變化曲線.可以看出單元能量守恒,從能量的角度驗證了模型的正確性.圖13 是梁中心點處軸向應變變化情況.ANCF/CRBF10 單元和ANCF/CRBF8 單元軸向均可以伸展,數值變化較為吻合,而ANCF/CRBF6 單元完全約束了軸向變形,因此無軸向應變變化.另外可以發現LSDYNA BS beam 單元軸向應變與ANCF/CRBF梁單元軸向應變的變化趨勢相同,但是幅值稍微偏大.這可以歸因于四點:第一,LSDYNA BS beam 單元與ANCF/CRBF 梁單元采用了根本不同的位移場插值策略.ANCF/CRBF 梁單元選擇對軸向和截面方向采用統一的插值策略,使得軸向、彎曲和剪切變形模態高度耦合,導致了系統剛度的上升,造成了軸向應變的減小.LSDYNA BS beam 單元基于共旋坐標法,對位移和轉動采用了獨立插值的策略,各個變形模態不會產生高度耦合現象,單元的系統剛度小于ANCF/CRBF 梁單元,所以軸向應變變化較大.第二,LSDYNA BS beam 單元與ANCF/CRBF 單元采用不同的坐標系統.由于ANCF/CRBF10 單元含有梯度向量以及ANCF/CRBF8 單元含有伸展系數,因此在進行邊界轉動副約束時也與LSDYNA BS beam 單元的邊界約束條件不同.第三,LSDYNA BS beam 單元和ANCF/CRBF 單元的求解過程不同.LSDYNA 軟件采用增量求解過程,而ANCF/CRBF 采用非增量求解過程.第四,兩種模型采用的顯示積分方法分別為MATLAB 和LSDYNA 軟件自設的算法,在誤差判斷和步長選擇上有所差異,隨著模擬時間的延長,兩種解之間的差異是可以預見的.

圖14 是梁中心點處剪切應變變化情況.由于ANCF/CRBF8 單元和ANCF/CRBF6 單元節點處同時施加了軸向和截面方向的運動學約束,因此節點處不存在剪切變形.ANCF/CRBF10 單元僅約束了截面的變形,未約束截面與軸線的轉動,因此節點處存在剪切效應.

圖14 梁中點剪切應變變化Fig.14 Shear strain change of beam middle point

相比傳統ANCF/CRBF10 單元,本文構造的ANCF/CRBF8 和ANCF/CRBF6 單元的自由度更少,使得模型維度更低.例如,在柔性單擺算例中,ANCF/CRBF10,ANCF/CRBF8 和ANCF/CRBF6 單元建立的模型自由度分別為105,84 和63,仿真時間分別為139.2 min,121.3 min 和48.7 min(仿真環境為16 GB 內存,Intel Core i7 計算平臺).此外,相比傳統ANCF/CRBF10 單元,本文構造的ANCF/CRBF8 和ANCF/CRBF6 單元坐標沒有梯度向量,可直接對其轉角參數施加運動約束或力約束,便于面向旋轉的設計操作.

湯惠穎等[46]基于李群SE(3)局部標架,研究了幾類梁單元的閉鎖緩解方法.文中分別采用Hu-Washizu 三場變分原理和應變分解法,緩解了SE(3)局部標架幾何精確梁單元的剪切閉鎖和絕對節點坐標梁單元的泊松閉鎖,提升了單元的計算精度.對比本文與文獻[46]可以發現:(1)文獻[46]的特點是引入李群SE(3)局部標架,可消除剛體運動帶來的幾何非線性[47],極大地減少迭代矩陣的更新次數,有效地提升了計算效率;本文的特點是引入單元截面的運動約束,可消除截面處的泊松閉鎖,縮減了單元的自由度,在低維動力學問題上提升了計算效率.(2)文獻[46]與本文均采用應變分解法處理絕對節點坐標單元的泊松閉鎖問題,在靜力學和動力學問題上的計算精度均得到了較好地提升.

另外,湯惠穎等[46]發現高階幾何精確梁單元比低階幾何精確梁單元的計算精度高.這可啟發本研究進一步拓展到絕對節點坐標高階梁研究方面,例如三維三節點梁單元[48]和三維高階梁單元[49-50].由于絕對節點坐標高階梁沒有閉鎖問題,因此可結合本文的一致旋轉場列式方法開發出精度較高的三維ANCF/CRBF 高階梁單元.

5 結論

本文系統地研究了基于一致旋轉場列式的絕對節點坐標平面梁單元的泊松閉鎖問題.文中首先開發了兩種剛性截面的絕對節點坐標平面梁單元ANCF/CRBF8 和ANCF/CRBF6,并對比了泊松閉鎖問題在新單元與已有ANCF/CRBF10 梁單元上的特點.隨后,結合兩種閉鎖緩解技術,彈性線方法和應變分解方法,研究了單元的收斂性.通過典型的靜力學和動力學算例測試了泊松閉鎖問題對ANCF/CRBF 單元性能的影響,可得到以下結論:

(1)對比ANCF 全參數梁單元,ANCF/CRBF 單元對節點梯度向量進行了約束,導致節點處不能完全捕捉泊松效應,但單元內部可以完全捕捉泊松效應,這種不連續情況會加重單元整體的泊松閉鎖問題,極大地降低靜力學計算精度.ANCF/CRBF 單元不具有靜力學的計算優勢,特別是在單元受集中力矩的強幾何非線性情況.

(2)對比三種ANCF/CRBF 單元,對節點梯度向量約束的越多,導致計算深梁結構的靜力學精度越差.對于長寬比較大的細長梁結構,梁的彎曲變形不會導致軸向伸展,剛性截面致使整體單元近似滿足歐拉梁假設.

(3)在ANCF/CRBF 單元的基礎上結合應變分解方法,可以更加有效地緩解單元泊松閉鎖問題,動力學計算結果具有良好的精度,ANCF/CRBF 單元更適用于模擬動力學問題.