張大羽 羅建軍 王輝 馬小飛

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

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

引言

柔性梁被廣泛地應(yīng)用于航天與機(jī)械系統(tǒng),是柔性多體系統(tǒng)動力學(xué)最典型的一類研究對象,建立其動力學(xué)模型可為復(fù)雜柔性多體系統(tǒng)設(shè)計和分析提供重要的理論依據(jù).建立柔性梁模型的關(guān)鍵是如何準(zhǔn)確地描述梁單元連續(xù)的旋轉(zhuǎn)場,特別是當(dāng)柔性梁經(jīng)歷大轉(zhuǎn)動大變形運(yùn)動.區(qū)分不同類型的非線性柔性梁建模方法的一個重要特征是如何準(zhǔn)確地描述系統(tǒng)中的有限轉(zhuǎn)動[1-10].

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

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

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

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

1 傳統(tǒng)的ANCF 全參數(shù)平面梁單元

傳統(tǒng)的ANCF 全參數(shù)平面梁單元是絕對節(jié)點坐標(biāo)方法中的經(jīng)典單元,其采用統(tǒng)一的插值函數(shù)描述單元的平動和轉(zhuǎn)動.該單元在運(yùn)動過程中橫截面與梁軸線不保持時刻垂直,截面可變形.

1.1 單元運(yùn)動學(xué)描述

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

1.2 單元彈性力公式

圖1 三種變形構(gòu)型之間的關(guān)系Fig.1 Three general configurations

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

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

1.3 單元動力學(xué)方程

2 ANCF/CRBF 平面梁單元

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

2.1 傳統(tǒng)的ANCF/CRBF 剛性截面梁單元

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

對于梁單元i的轉(zhuǎn)換關(guān)系如下[43]

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

對于梁單元i的轉(zhuǎn)換關(guān)系如下[43]

2.4 單元動力學(xué)方程

值得注意的是,ANCF/CRBF 單元計算轉(zhuǎn)角是通過積分角度對時間的導(dǎo)數(shù)得到的,而非積分角速度.分空間和平面兩種情況進(jìn)行討論.

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

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

3 單元閉鎖緩解技術(shù)

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

3.1 彈性線方法

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

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

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

3.2 應(yīng)變分解方法

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

依據(jù)Voigt 標(biāo)記規(guī)則,可寫為

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

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

4 算例分析

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

(1)ANCF12-GCM (傳統(tǒng)ANCF 全參數(shù)梁單元+連續(xù)介質(zhì)力學(xué));

(2)ANCF/CRBF10-GCM (傳統(tǒng)剛性截面ANCF/CRBF 梁單元+連續(xù)介質(zhì)力學(xué));

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

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

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

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

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

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

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

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

根據(jù)第2 節(jié)和第3 節(jié)對絕對節(jié)點坐標(biāo)全參數(shù)梁和基于一致旋轉(zhuǎn)場的絕對節(jié)點坐標(biāo)梁單元的描述,詳細(xì)分析各類型單元泊松閉鎖的特點.

(1)ANCF12-GCM

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

②旋轉(zhuǎn)描述中無轉(zhuǎn)角參數(shù),僅有梯度向量;

③節(jié)點處存在剪切變形;

④截面無約束,節(jié)點處和單元內(nèi)均有泊松閉鎖.

(2)ANCF/CRBF10-GCM

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

②旋轉(zhuǎn)描述中既有轉(zhuǎn)角參數(shù)又有梯度向量;

③節(jié)點處存在剪切變形;

④節(jié)點截面全約束,節(jié)點處無泊松閉鎖,單元內(nèi)存在泊松閉鎖.這種不連續(xù)情況導(dǎo)致單元整體的泊松閉鎖問題相較ANCF12 單元更加嚴(yán)重.

(3)ANCF/CRBF8-GCM

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

②旋轉(zhuǎn)描述中僅有轉(zhuǎn)角參數(shù);

③節(jié)點處不存在剪切變形;

④節(jié)點截面全約束,節(jié)點處無泊松閉鎖,單元內(nèi)存在泊松閉鎖.此外,節(jié)點處無剪切變形,單元內(nèi)存在剪切變形.這種不連續(xù)情況導(dǎo)致單元整體的泊松閉鎖相較ANCF10 單元更加嚴(yán)重.

(4)ANCF/CRBF6-GCM

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

②旋轉(zhuǎn)描述中僅有轉(zhuǎn)角參數(shù);

③節(jié)點處不存在剪切變形;

④節(jié)點截面全約束,節(jié)點處無泊松閉鎖,單元內(nèi)存在泊松閉鎖.此外,單元軸線不能伸展.這會導(dǎo)致單元整體的泊松閉鎖相較ANCF8 單元更加嚴(yán)重.

4.1 靜力學(xué)測試:小變形細(xì)長梁算例

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

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

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

表2 小變形情況下Y 向位移對比(使用閉鎖緩解技術(shù))(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 靜力學(xué)測試:大變形深梁算例

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

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

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

表4 大變形情況下Y 向位移對比(使用閉鎖緩解技術(shù))(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 靜力學(xué)測試:懸臂梁受端部集中力矩算例

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

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

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

表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 懸臂梁受端部集中力矩下的構(gòu)型圖(40 個單元)Fig.6 Cantilever beam configurations subjected by a torque(40 elements)

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

圖7 ry 范數(shù)的誤差Fig.7 Error of the norm of ry

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

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

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

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

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

4.5 動力學(xué)測試:柔性單擺算例

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

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

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

圖14 梁中點剪切應(yīng)變變化Fig.14 Shear strain change of beam middle point

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

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

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

5 結(jié)論

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

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

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

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