兩類(lèi)基于局部標(biāo)架梁?jiǎn)卧拈]鎖緩解方法1)

湯惠穎 張志娟 劉 鋮,2) 劉紹奎

?(北京理工大學(xué)宇航學(xué)院,北京 100081)

?(北京空間飛行器總體設(shè)計(jì)部,北京 100094)

引言

幾何精確方法[1](geometrically exact formulation,GEF)與絕對(duì)節(jié)點(diǎn)坐標(biāo)方法[2](absolute nodal coordinate formulation,ANCF)是描述大轉(zhuǎn)動(dòng)、大變形柔性多體系統(tǒng)的兩類(lèi)常用的建模方法.20 世紀(jì)80 年中期,李群/李代數(shù)基本理論逐漸與非線(xiàn)性有限元方法融合,以解決梁、板/殼結(jié)構(gòu)的大轉(zhuǎn)動(dòng)、大變形動(dòng)力學(xué)問(wèn)題.幾何精確建模方法在有限元理論框架下利用節(jié)點(diǎn)位置矢量與轉(zhuǎn)動(dòng)偽矢量作為節(jié)點(diǎn)參數(shù)描述三維梁大轉(zhuǎn)動(dòng)、大變形運(yùn)動(dòng),對(duì)位置矢量與轉(zhuǎn)動(dòng)偽矢量獨(dú)立插值[3].其中涉及的兩個(gè)核心算法在于:第一、剛體旋轉(zhuǎn)矩陣在李群SO(3)內(nèi)采用乘法更新,對(duì)應(yīng)的轉(zhuǎn)動(dòng)偽矢量由對(duì)數(shù)映射確定,滿(mǎn)足剛體轉(zhuǎn)動(dòng)合成的幾何意義[4].第二、通過(guò)轉(zhuǎn)動(dòng)參數(shù)定義梁截面局部標(biāo)架,并在局部標(biāo)架下研究梁微元的動(dòng)能與彈性勢(shì)能,即拉伸、剪切、扭轉(zhuǎn)、彎曲應(yīng)變以及截面角速度均定義在局部標(biāo)架內(nèi)[5-6].在多體系統(tǒng)動(dòng)力學(xué)領(lǐng)域,Shabana[7]于1996 年提出了絕對(duì)節(jié)點(diǎn)坐標(biāo)方法,其本質(zhì)上也屬于一種非線(xiàn)性有限元方法[8].該方法在處理物體轉(zhuǎn)動(dòng)問(wèn)題上,摒棄復(fù)雜的轉(zhuǎn)動(dòng)參數(shù),采用節(jié)點(diǎn)位置矢量以及沿物質(zhì)坐標(biāo)的斜率矢量來(lái)描述物體轉(zhuǎn)動(dòng),可方便地建立大轉(zhuǎn)動(dòng)、大變形柔性多體系統(tǒng)動(dòng)力學(xué)模型[9-10].

已有研究表明,在局部標(biāo)架思想下,可構(gòu)造出一類(lèi)SE(3)幾何精確梁?jiǎn)卧猍11]以及絕對(duì)節(jié)點(diǎn)坐標(biāo)梁?jiǎn)卧?局部標(biāo)架梁?jiǎn)卧軌蛞?guī)避剛體運(yùn)動(dòng)帶來(lái)的幾何非線(xiàn)性問(wèn)題,離散數(shù)值模型中廣義質(zhì)量矩陣與切線(xiàn)剛度矩陣滿(mǎn)足剛體變換的不變性,可明顯地提高柔性多體系統(tǒng)動(dòng)力學(xué)問(wèn)題的計(jì)算效率.在有限元方法中,梁、板/殼單元普遍存在剪切、薄膜以及泊松閉鎖問(wèn)題.多年來(lái),有限元領(lǐng)域?qū)W者提出了許多單元技術(shù)來(lái)有效地緩解閉鎖.本文主要關(guān)注如何提高局部標(biāo)架梁?jiǎn)卧氖諗啃?關(guān)于李群局部標(biāo)架建模與計(jì)算方法涉及的其它細(xì)節(jié)本文不再贅述,可參見(jiàn)作者的相關(guān)研究工作.

縮減積分是解決閉鎖問(wèn)題最廣泛使用的技術(shù)之一.Zienkiewicz 等[12]最早將縮減積分應(yīng)用于殼單元、等參四邊形單元和梁?jiǎn)卧猍13],提高了單元精度.Hughes 等[14]將選擇縮減積分應(yīng)用于板單元,緩解了薄板的剪切閉鎖.Simo 等[15]采用縮減積分技術(shù)緩解了幾何精確梁中的剪切閉鎖.Malkus 和Hughes[16]證明了某些混合列式與使用縮減積分的位移公式的等價(jià)性,并且表明縮減積分單元可以達(dá)到多變量有限元的性能.Noor 和Peters[17]在曲梁中應(yīng)用選擇縮減積分,并通過(guò)比較剛度矩陣討論了混合列式和位移列式的等價(jià)性和“近等價(jià)性”.

混合法不僅假設(shè)位移場(chǎng),同時(shí)假設(shè)應(yīng)力場(chǎng)或應(yīng)變場(chǎng),是一種高性能的多變量有限元方法,已被廣泛用于解決體積閉鎖、剪切閉鎖和薄膜閉鎖等問(wèn)題.1985 年,Pian[18]提出了一種基于Hellinger-Reissner兩場(chǎng)變分原理的混合單元,采用假設(shè)應(yīng)力和位移場(chǎng)改善了梁和板的彎曲性能.1988 年,Liu 等[19]提出了基于Hu-Washizu 三場(chǎng)變分原理的彎曲超收斂單元,該單元在粗網(wǎng)格中顯示出良好的精度.1994 年,Dorfi和Busby[20]提出了一種基于Hellinger-Reissner 兩場(chǎng)變分原理的混合曲梁?jiǎn)卧?該單元的位移和應(yīng)力具有良好的收斂性.1990 年,Simo 和Rifai[21]提出了增強(qiáng)假設(shè)應(yīng)變法,該方法是一種具有變分基礎(chǔ)的閉鎖緩解技術(shù).隨后,Simo 和Armero[22]將增強(qiáng)應(yīng)變法推廣到了幾何非線(xiàn)性問(wèn)題中.1993 年,Andelfinger 和Ramm[23]利用增強(qiáng)應(yīng)變法開(kāi)發(fā)了二維和三維板和殼單元,并證明了該方法與基于Hellinger-Reissner 兩場(chǎng)變分原理的假設(shè)應(yīng)力法[24]是等價(jià)的.

ANCF 全參數(shù)(fully parameterized)梁?jiǎn)卧谌S連續(xù)介質(zhì)力學(xué)方法計(jì)算廣義彈性力,當(dāng)泊松比不為零時(shí),其計(jì)算結(jié)果不能收斂到正確解.該問(wèn)題即為泊松閉鎖[25].起初,學(xué)者們令泊松比為零[26]或通過(guò)簡(jiǎn)化應(yīng)力張量[27]除去泊松效應(yīng),緩解了泊松閉鎖.隨后,Gerstmayr 等[28]采用選擇縮減積分技術(shù)緩解了泊松閉鎖,同時(shí)在單元的變形模式中保留了泊松效應(yīng).2010 年,Matikainen 等[29]提出了一種高階三維ANCF 梁?jiǎn)卧?該單元對(duì)截面進(jìn)行二次插值,有效地緩解了泊松閉鎖,但是單元自由度過(guò)多.2018 年,Patel 和Shabana[30]提出了一種新型閉鎖緩解技術(shù)——應(yīng)變分解法(strain split method,SSM).該方法通過(guò)分解Green-Lagrange 應(yīng)變張量并修改本構(gòu)模型,假設(shè)僅與梁中線(xiàn)變形相關(guān)的低階應(yīng)變具有泊松效應(yīng),忽略正應(yīng)變的高階項(xiàng)的泊松效應(yīng),顯著改善了梁的彎曲性能.

本文研究局部標(biāo)架下幾類(lèi)梁?jiǎn)卧拈]鎖緩解方法.采用Hu-Washizu 三場(chǎng)變分原理緩解局部標(biāo)架下李群SE(3)幾何精確梁的剪切閉鎖,采用應(yīng)變分解法緩解局部標(biāo)架下ANCF 全參數(shù)梁的泊松閉鎖.通過(guò)對(duì)比幾何精確與絕對(duì)節(jié)點(diǎn)坐標(biāo)兩類(lèi)建模方法,分析基于局部標(biāo)架的李群SE(3)幾何精確梁在計(jì)算精度與效率方面的優(yōu)勢(shì).本文中的李群SE(3)幾何精確梁與ANCF 梁?jiǎn)卧诰植繕?biāo)架,為了方便起見(jiàn),下文省去“基于局部標(biāo)架”幾個(gè)字.

1 基于李群SE(3)局部標(biāo)架的幾何精確梁理論

1.1 梁的運(yùn)動(dòng)學(xué)描述及本構(gòu)模型

基于經(jīng)典的Timoshenko 梁假設(shè),李群幾何精確梁理論融合非線(xiàn)性有限元方法與幾何力學(xué)思想,在李群SE(3)內(nèi)建立系統(tǒng)的平衡方程,可精確地描述柔性梁的大轉(zhuǎn)動(dòng)、大變形剛?cè)狁詈蟿?dòng)力學(xué)特性,同時(shí)可規(guī)避剛體運(yùn)動(dòng)帶來(lái)的幾何非線(xiàn)性[31].

梁的初始構(gòu)型和當(dāng)前構(gòu)型如圖1 所示.引入?yún)⒖甲鴺?biāo)系{O;e1,e2,e3}.在當(dāng)前構(gòu)型中,在梁中線(xiàn)的每一點(diǎn)定義隨體基矢量(t,n,b),引入梁中線(xiàn)位置矢量x和截面旋轉(zhuǎn)矩陣R,初始構(gòu)型對(duì)應(yīng)的量為(·)0.

圖1 幾何精確梁初始構(gòu)型與當(dāng)前構(gòu)型Fig.1 Initial and current configurations of the geometrically exact beam

在梁的當(dāng)前構(gòu)型中,截面上任意一點(diǎn)p的位置矢量為

式中,彈性系數(shù)矩陣D=diag(EA,ks2GA,ks3GA,GJs,EI2,EI3),其中E和G分別為楊氏模量和剪切模量,A,I2和I3分別表示梁初始構(gòu)型的截面面積和截面慣性矩,Js表示圣維南扭轉(zhuǎn)常數(shù),ks2和ks3為截面的剪切修正系數(shù).

1.2 SE3 幾何精確梁控制方程及線(xiàn)性化

系統(tǒng)的動(dòng)力學(xué)平衡方程為

式中,δWint為內(nèi)力所做的虛功,δWext為外載荷所做的虛功,δWine為慣性力所做的虛功.

為獲得一個(gè)對(duì)稱(chēng)的切線(xiàn)剛度矩陣,本文在前一個(gè)收斂步所在切平面進(jìn)行物理量的一次與二次變分.因此,應(yīng)變的變分可以表示為

式中,T 為SE(3)群的切空間算符,δh=[δdTδηT]T,δd和δη分別為局部標(biāo)架下的虛位移和虛轉(zhuǎn)角,B 為應(yīng)變算符.內(nèi)部虛功可以表示為

外部虛功可以表示為

式中,pext為作用于p點(diǎn)的外載荷.慣性虛功可以表示為

式中,v=[UTωT]T,U和ω分別為局部標(biāo)架下的線(xiàn)速度與角速度,頂標(biāo)·表示變量對(duì)時(shí)間t的全導(dǎo)數(shù),ρA為單位長(zhǎng)度的材料密度,I3表示3×3 的單位矩陣,J為截面形心慣性矩陣.

對(duì)內(nèi)部虛功做線(xiàn)性化可得

應(yīng)變算符B 的線(xiàn)性化的轉(zhuǎn)置為

為了進(jìn)一步簡(jiǎn)化?(δWint)的幾何部分,引入兩個(gè)矩陣ΞT和ΞT′,它們的定義如下

式中,A為一個(gè)任意的6×1 的矢量,記號(hào)和分別是矩陣ΞT′的前半部分和后半部分,其具體形式可參見(jiàn)文獻(xiàn)[11].因此,幾何部分為

最后,內(nèi)部虛功的線(xiàn)性化可以寫(xiě)為

需要指出的是,由于δh為局部標(biāo)架下的平動(dòng)與轉(zhuǎn)動(dòng)虛位移,?h為局部標(biāo)架下的平動(dòng)與轉(zhuǎn)動(dòng)無(wú)限小位移,它們均與該點(diǎn)的任意剛體位移無(wú)關(guān).因此,單元的彈性力及其切線(xiàn)剛度矩陣均滿(mǎn)足剛體變換的不變性.對(duì)于廣義慣性力及其廣義質(zhì)量矩陣也同時(shí)具有類(lèi)似的形式.

1.3 SE3 幾何精確梁有限元離散格式

對(duì)梁中線(xiàn)位置矢量x和截面相對(duì)旋轉(zhuǎn)偽矢量Θr[11]進(jìn)行空間離散

式中,NI是節(jié)點(diǎn)I的形函數(shù),xI是梁中線(xiàn)上節(jié)點(diǎn)I的位置矢量,Rr是參考點(diǎn)所在截面的旋轉(zhuǎn)矩陣,RI是節(jié)點(diǎn)I所在截面的旋轉(zhuǎn)矩陣,旋轉(zhuǎn)偽矢量為轉(zhuǎn)軸單位矢量與轉(zhuǎn)角的乘積.

離散應(yīng)變算符B 得到應(yīng)變?位移矩陣BI,它的表達(dá)式為

彈性力可以表示為

材料剛度和幾何剛度矩陣可以表示為

對(duì)于非常細(xì)長(zhǎng)的梁結(jié)構(gòu),采用一階插值描述幾何精確梁彎曲時(shí)將會(huì)產(chǎn)生剪切閉鎖.由于梁非常細(xì)長(zhǎng),截面內(nèi)橫向剪切應(yīng)變處處接近于0.而由于平動(dòng)參數(shù)x和轉(zhuǎn)動(dòng)參數(shù)θ采用了相同的插值函數(shù),使式(5)中橫向剪切應(yīng)變?chǔ)胹中的R(θ)和x′兩項(xiàng)具有不同階次,γs=0 不可能處處滿(mǎn)足,導(dǎo)致能量泛函中剪切變形能項(xiàng)的量級(jí)不正確,由此帶來(lái)了剪切閉鎖.此時(shí),單元表現(xiàn)得十分剛硬,即計(jì)算出來(lái)的位移值遠(yuǎn)小于參考解,總體剛度矩陣很大.

2 基于局部標(biāo)架的絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的三維全參數(shù)梁理論

本節(jié)簡(jiǎn)要介紹基于局部標(biāo)架的絕對(duì)節(jié)點(diǎn)坐標(biāo)全參數(shù)梁?jiǎn)卧７椒?某質(zhì)點(diǎn)從初始位置X運(yùn)動(dòng)到當(dāng)前位置r,變形梯度可以表示為

Green-Lagrange 應(yīng)變張量可以表示為

應(yīng)變能可以表示為

式中,D為彈性矩陣,V為梁的體積.應(yīng)變能對(duì)廣義坐標(biāo)e求偏導(dǎo)得彈性力為

彈性力對(duì)廣義坐標(biāo)求偏導(dǎo)得剛度陣為

相比于幾何精確梁?jiǎn)卧?ANCF 單元沒(méi)有轉(zhuǎn)動(dòng)參數(shù),不能直接定義局部標(biāo)架,需要通過(guò)平動(dòng)位移場(chǎng)定義轉(zhuǎn)動(dòng)場(chǎng).此處采用ANCF 斜率矢量定義該點(diǎn)的局部標(biāo)架基矢量

式中,rX和rY為當(dāng)前位置矢量r對(duì)物質(zhì)坐標(biāo)的偏導(dǎo)數(shù).由于平動(dòng)位移場(chǎng)在局部標(biāo)架中描述,單元的剛體運(yùn)動(dòng)可自然被消除.因此,系統(tǒng)質(zhì)量矩陣以及切線(xiàn)剛度矩陣滿(mǎn)足剛體運(yùn)動(dòng)的不變性,從而可極大地減少系統(tǒng)質(zhì)量矩陣與剛度矩陣在仿真過(guò)程中的更新次數(shù).

ANCF 全參數(shù)梁?jiǎn)卧獣?huì)遇到泊松閉鎖問(wèn)題,閉鎖產(chǎn)生的原因如下.由于ANCF 全參數(shù)梁是連續(xù)介質(zhì)單元,梁的中線(xiàn)為三階Hermite 插值,在泊松效應(yīng)的作用下,截面應(yīng)該變形為曲面.然而橫截面為一階插值,不能描述截面的曲面變形模式,因此導(dǎo)致泊松閉鎖.

3 兩類(lèi)局部標(biāo)架梁?jiǎn)卧拈]鎖緩解方法

3.1 基于Hu-Washizu 三場(chǎng)變分原理的剪切閉鎖緩解方法

Hu-Washizu 三場(chǎng)變分原理不僅假設(shè)位移場(chǎng),同時(shí)假設(shè)應(yīng)力和應(yīng)變場(chǎng),且三類(lèi)場(chǎng)變量的插值相互獨(dú)立.假設(shè)剪切應(yīng)變和假設(shè)剪切應(yīng)力的分布為

式中,Ss=[0,,0,0,0]T.剪切部分的彈性力可以表示為

剪切部分的剛度矩陣可以表示為

式中,

離散方程組可以表示為

式中,Km和Kb分別為薄膜和彎曲部分的剛度矩陣,?hN為節(jié)點(diǎn)運(yùn)動(dòng)矢量的增量,Res為殘差.將上式靜力縮聚化為單場(chǎng)形式

一階插值的幾何精確梁會(huì)產(chǎn)生嚴(yán)重的剪切閉鎖,二階及二階以上插值的幾何精確梁剪切閉鎖較輕.因此,本文關(guān)于幾何精確梁剪切閉鎖處理方法著重于一階梁?jiǎn)卧?Hu-Washizu 三場(chǎng)變分原理中,位移場(chǎng)、應(yīng)變場(chǎng)和應(yīng)力場(chǎng)均獨(dú)立插值.式(28)假設(shè)橫向剪切應(yīng)變?yōu)槌?shù),可避免產(chǎn)生過(guò)大的剪切應(yīng)變能,有效緩解了剪切閉鎖.

此外,對(duì)于梁?jiǎn)卧?縮減積分技術(shù)也是緩解剪切閉鎖的有效手段.數(shù)值算例展示了采用Hu-Washizu三場(chǎng)變分原理與縮減積分技術(shù)緩解梁?jiǎn)卧羟虚]鎖的效果.相比較而言,縮減積分技術(shù)由于其簡(jiǎn)便性在梁?jiǎn)卧袘?yīng)用較為廣泛,但用在板殼單元中會(huì)產(chǎn)生沙漏變形模式.進(jìn)一步,上述Hu-Washizu 三場(chǎng)變分原理的研究將為低階幾何精確板殼單元的閉鎖緩解起到一定的借鑒作用.

3.2 基于應(yīng)變分解法的泊松閉鎖緩解方法

ANCF 全參數(shù)梁?jiǎn)卧奈恢脠?chǎng)和變形梯度可以寫(xiě)為

εc中只有低階應(yīng)變,且εc僅與梁中線(xiàn)參數(shù)X相關(guān);εk中包含高階應(yīng)變,這些高階應(yīng)變會(huì)導(dǎo)致截面變形和彎曲變形.為了緩解閉鎖,泊松耦合只用在正應(yīng)變的低階項(xiàng)中,即僅在εc中考慮泊松效應(yīng),忽略正應(yīng)變的高階項(xiàng)的泊松效應(yīng).因此,梁截面不會(huì)變形為曲面,緩解了泊松閉鎖.第二類(lèi)Piola-Kirchhoff應(yīng)力的Voigt形式可以寫(xiě)為

式中,應(yīng)變也是Voigt 形式.基于平面應(yīng)變假設(shè),彈性系數(shù)矩陣可以被定義為

式中,λ=Eν/[(1+ν)(1 ?2ν)]和μ=E/[2(1+ν)]為拉梅常數(shù),E和ν 分別為楊氏模量和泊松比.應(yīng)變能可以表示為

彈性力可以表示為

材料和幾何剛度可以表示為

4 三類(lèi)常用的有限元插值方法

本節(jié)簡(jiǎn)要介紹三類(lèi)有限元插值方法,包括Lagrange 插值、Hermite 插值以及非均勻有理B-樣條(NURBS).其中,Lagrange 插值在有限元中應(yīng)用最為廣泛;絕對(duì)節(jié)點(diǎn)坐標(biāo)方法則采用Hermite 插值,以保證位移場(chǎng)的C1連續(xù);NURBS 的基函數(shù)可以構(gòu)造任意階連續(xù)的近似函數(shù),建立了計(jì)算機(jī)輔助設(shè)計(jì)(CAD)與計(jì)算機(jī)輔助工程(CAE)間的橋梁,該方法被稱(chēng)為等幾何分析.

為對(duì)比分析不同插值方法的優(yōu)劣勢(shì),本文涉及一階與三階Lagrange 插值、三階Hermite 插值以及三階NURBS 插值.

4.1 Lagrange 插值

給定函數(shù)f(x)在n+1 個(gè)互不相同的點(diǎn)xi(i=0,1,2,···,n)上的函數(shù)值f(xi)=yi,n次Lagrange 插值多項(xiàng)式可以表示為

4.2 Hermite 插值

給定函數(shù)f(x)在n+1 個(gè)互不相同的點(diǎn)xi(i=0,1,2,···,n)上的函數(shù)值yi及一階導(dǎo)數(shù)值,Hermite插值多項(xiàng)式可以表示為

Hermite 插值基函數(shù)Ai(x)和Bi(x)為

4.3 非均勻有理B-樣條

在區(qū)間[a,b]內(nèi),設(shè)U={u0,u1,u2,···,um}是一不減的實(shí)數(shù)序列.以U作為樣條結(jié)點(diǎn),p次的第i個(gè)B-樣條基函數(shù)Ni,p(u)定義為

式中,Pi是控制點(diǎn),ωi是權(quán)系數(shù).

5 直/曲梁結(jié)構(gòu)大變形靜力學(xué)分析

本節(jié)考察XOY平面內(nèi)的直梁與曲梁模型,驗(yàn)證上述閉鎖處理方法的有效性.其中,懸臂直梁長(zhǎng)L=1 m;懸臂曲梁半徑R=1 m,初始構(gòu)型為1/4圓.梁的截面均為矩形,寬和高為0.01 m;材料的楊氏模量為E=2.1×1011Pa,泊松比為ν=0.3.首先,在直梁和曲梁末端分別施加集中力(30,20,400)N (直梁)、(30,20,100)N (曲梁),用ABAQUS 計(jì)算得到梁端點(diǎn)沿z方向的位移分別為0.5143 m (直梁)、0.6138 m(曲梁).計(jì)算的過(guò)程中,上述兩類(lèi)載荷均分為10 個(gè)加載步均勻加載.此外,對(duì)直梁與曲梁進(jìn)行純彎曲測(cè)試.在直梁和曲梁末端均施加z方向的力矩Mz=?πEI/2L,由材料力學(xué)解答知,梁端沿z方向轉(zhuǎn)角位移的解析解為θz=?π/2.

5.1 懸臂直梁和曲梁低階單元閉鎖緩解方法數(shù)值測(cè)試

本節(jié)考查SE(3)幾何精確一階直梁與曲梁?jiǎn)卧谌S復(fù)雜應(yīng)力狀態(tài)下,縮減積分技術(shù)和Hu-Washizu三場(chǎng)變分原理緩解剪切閉鎖的能力.

懸臂直梁和曲梁末端受集中力的作用,變形如圖2 所示.單元位移場(chǎng)采用一階Lagrange 插值.對(duì)比分析完全積分(exact integration,EI)、縮減積分(reduced integration,RI)和Hu-Washizu 三場(chǎng)變分原理(Hu-Washizu variational principle,HWVP)三種情況,在不同自由度(degrees of freedom,Dofs)下懸臂梁末端點(diǎn)z方向的位移uz,計(jì)算相對(duì)誤差,將結(jié)果列于表1和表2 中.同時(shí),采用ABAQUS 軟件中B31 梁?jiǎn)卧?該單元為三維二節(jié)點(diǎn)一階梁?jiǎn)卧?每個(gè)節(jié)點(diǎn)6 個(gè)廣義坐標(biāo),分別是3 個(gè)平移和3 個(gè)轉(zhuǎn)動(dòng)參數(shù))對(duì)本算例進(jìn)行仿真計(jì)算.

圖2 (a)直梁和(b)曲梁變形前后示意圖Fig.2 Initial and deformed configurations of(a)the straight beam and(b)the curved beam

表1 懸臂直梁末端z 方向位移和相對(duì)誤差Table 1 End displacements of z direction and relative error of the straight cantilever beam

表2 懸臂曲梁末端z 方向位移和相對(duì)誤差Table 2 End displacements of z direction and relative error of the curved cantilever beam

表1 和表2 的數(shù)值結(jié)果表明:采用縮減積分技術(shù)的SE(3)幾何精確梁與ABAQUS 中B31 單元收斂結(jié)果基本一致.采用Hu-Washizu 三場(chǎng)變分原理緩解閉鎖的SE(3)幾何精確梁與其他兩類(lèi)方法的收斂結(jié)果不同,其中,直梁模型收斂值偏大,曲梁模型收斂值偏小.在三維復(fù)雜應(yīng)力狀態(tài)下,一階直梁和曲梁均有剪切閉鎖問(wèn)題.緩解閉鎖后,直梁的計(jì)算精度提高了98%以上,曲梁的計(jì)算精度提高了52%以上,效果十分顯著.對(duì)于一階直梁,以上兩種緩解閉鎖的方法計(jì)算精度無(wú)明顯差別.對(duì)于一階曲梁,縮減積分技術(shù)的計(jì)算精度略高于Hu-Washizu 三場(chǎng)變分原理.

5.2 懸臂直梁和曲梁純彎曲測(cè)試

本節(jié)考查SE(3)幾何精確一階直梁與曲梁?jiǎn)卧诩儚澢鷳?yīng)力狀態(tài)下,縮減積分技術(shù)和Hu-Washizu三場(chǎng)變分原理緩解剪切閉鎖的能力,同時(shí)評(píng)價(jià)單元描述純彎曲的能力.

在懸臂直梁和曲梁末端施加集中彎矩,此時(shí),懸臂梁發(fā)生純彎曲,處于簡(jiǎn)單應(yīng)力狀態(tài).單元位移場(chǎng)采用一階Lagrange 插值.對(duì)比分析完全積分、縮減積分和Hu-Washizu 三場(chǎng)變分原理三種情況下,懸臂梁末端點(diǎn)z方向的轉(zhuǎn)角位移θz,計(jì)算相對(duì)誤差,將結(jié)果列于表3 和表4 中.同時(shí),采用ABAQUS 軟件中B31 梁?jiǎn)卧獙?duì)本算例進(jìn)行仿真計(jì)算.

表3 懸臂直梁末端z 方向角位移和相對(duì)誤差Table 3 Angular end displacements of z direction and relative error of the straight cantilever beam

表4 懸臂曲梁末端z 方向角位移和相對(duì)誤差Table 4 Angular displacements of z direction and relative error of the curved cantilever beam

表3 和表4 的數(shù)值結(jié)果表明:在純彎曲應(yīng)力狀態(tài)下,一階直梁和曲梁均有剪切閉鎖問(wèn)題.緩解閉鎖后,直梁和曲梁模型的計(jì)算精度均提高了98%以上,效果十分顯著.而且,以上兩種緩解閉鎖的方法和ABAQUS 中B31 單元的計(jì)算精度完全一致.表3中,一階直梁在任意單元數(shù)目下,轉(zhuǎn)角與解析解的相對(duì)誤差均為零.這是由于SE(3)幾何精確直梁模型能夠精確地構(gòu)造一個(gè)常彎曲單元,因此僅用1 個(gè)單元就能夠精確模擬純彎曲狀態(tài)下的轉(zhuǎn)動(dòng)場(chǎng).表4 中,一階曲梁模型的轉(zhuǎn)角漸近收斂,這是由于一階插值函數(shù)無(wú)法精確地描述一段圓弧.

同時(shí),采用三階Lagrange 插值構(gòu)造了高階幾何精確直梁與曲梁?jiǎn)卧?通過(guò)上述懸臂直梁與曲梁算例,測(cè)試了單元收斂精度,其中同樣采用Hu-Washizu三場(chǎng)變分原理處理剪切閉鎖.數(shù)值結(jié)果表明:處理閉鎖不會(huì)提升單元收斂性能,進(jìn)一步說(shuō)明了高階直梁與曲梁?jiǎn)卧淮嬖诩羟虚]鎖.在5.3 節(jié)中,將進(jìn)一步對(duì)比處理閉鎖后的一階單元與高階單元的收斂性.

5.3 幾類(lèi)局部標(biāo)架梁?jiǎn)卧諗啃詫?duì)比分析

本節(jié)對(duì)比一階和三階SE(3)幾何精確梁?jiǎn)卧?、一階R3×R3幾何精確梁、ANCF 縮減(slope deficiency)梁、ANCF 全參數(shù)梁與ABAQUS 軟件中B31 梁?jiǎn)卧氖諗啃?其中,一階幾何精確梁與ANCF 全參數(shù)梁?jiǎn)卧烟幚黹]鎖.

在懸臂直梁和曲梁末端施加集中力.計(jì)算以上幾類(lèi)建模方法在不同自由度下懸臂梁末端點(diǎn)z方向的位移uz和相對(duì)誤差,將結(jié)果列于表5 和表6 中.ANCF 縮減曲梁?jiǎn)卧獰o(wú)算進(jìn)行本算例的計(jì)算。

表5 懸臂直梁末端z 方向位移和相對(duì)誤差Table 5 End displacements of z direction and relative error of straight cantilever beams

表6 懸臂曲梁末端z 方向位移和相對(duì)誤差Table 6 End displacements of z direction and relative error of curved cantilever beams

表5 和表6 的數(shù)值結(jié)果表明:直梁模型的五種建模方法收斂結(jié)果基本一致.對(duì)于曲梁模型,SE(3)幾何精確梁、R3×R3幾何精確梁和ABAQUS 中B31 單元收斂結(jié)果基本一致,ANCF 全參數(shù)梁?jiǎn)卧諗拷Y(jié)果略大于其他三類(lèi)方法的收斂值.需要指出的是,表5 中三階NURBS 插值的SE(3)幾何精確梁?jiǎn)卧?4 個(gè)自由度時(shí),雖然計(jì)算結(jié)果的相對(duì)誤差為6.6×10?4,但此時(shí)單元還沒(méi)有收斂,60 個(gè)自由度時(shí),單元才收斂.對(duì)于高階單元,三階Lagrange 比三階NURBS 插值的SE(3)幾何精確梁的計(jì)算精度高,幾類(lèi)幾何精確梁均比ANCF 梁?jiǎn)卧挠?jì)算精度高.ANCF 縮減梁比ANCF 全參數(shù)梁?jiǎn)卧挠?jì)算精度高.

對(duì)于一階直梁?jiǎn)卧?兩類(lèi)幾何精確梁具有完全相同的計(jì)算精度.對(duì)于一階曲梁?jiǎn)卧?R3×R3幾何精確梁的計(jì)算精度略高于SE(3)幾何精確梁.對(duì)于一階直梁和曲梁?jiǎn)卧?兩類(lèi)幾何精確方法均比ABAQUS中B31 單元的計(jì)算精度高.顯然高階單元比低階單元精度高,但是高階單元計(jì)算效率低,即自由度相同時(shí),高階單元?jiǎng)偠染仃噧?nèi)非零元素的個(gè)數(shù)多.

6 高轉(zhuǎn)速/大變形柔性多體系統(tǒng)動(dòng)力學(xué)驗(yàn)證算例

本節(jié)通過(guò)高轉(zhuǎn)速曲柄滑塊機(jī)構(gòu)以及空間雙擺的動(dòng)力學(xué)仿真,分析局部標(biāo)架下的梁?jiǎn)卧趧?dòng)力學(xué)仿真中的計(jì)算精度差異,以及消除剛體運(yùn)動(dòng)帶來(lái)的幾何非線(xiàn)性后梁?jiǎn)卧臄?shù)值特性.由于商業(yè)軟件無(wú)法直接進(jìn)行大變形柔性多體系統(tǒng)動(dòng)力學(xué)仿真,因此本節(jié)不再采用商業(yè)軟件進(jìn)行仿真對(duì)比.為了設(shè)計(jì)特定的小變形或大變形的柔性多體系統(tǒng)動(dòng)力學(xué)模型,本節(jié)算例中部分材料參數(shù)未與真實(shí)材料對(duì)應(yīng).

本節(jié)數(shù)值算例共涉及如下3 類(lèi)建模方法:(1)一階和三階Lagrange 插值的SE(3)幾何精確梁?jiǎn)卧?(2)一階Lagrange 插值的R3×R3幾何精確梁?jiǎn)卧?(3)ANCF 全參數(shù)梁?jiǎn)卧?需要指出的是,除了三階SE(3)幾何精確梁?jiǎn)卧?其他單元均已處理閉鎖.

6.1 發(fā)生小變形和中等變形的高速旋轉(zhuǎn)的平面曲柄滑塊機(jī)構(gòu)

本節(jié)對(duì)一個(gè)典型的多體系統(tǒng)——曲柄滑塊機(jī)構(gòu)進(jìn)行動(dòng)力學(xué)仿真.如圖3 所示,平面曲柄滑塊機(jī)構(gòu)由兩根桿與一個(gè)滑塊組成,桿I 一端與地面采用球鉸相連接,另一端與桿II 一端也通過(guò)球鉸連接,桿II 末端與一個(gè)僅能在X軸上運(yùn)動(dòng)的剛體滑塊鉸接.

圖3 曲柄滑塊機(jī)構(gòu)Fig.3 A crank-slider mechanism

本算例中,通過(guò)桿II 中點(diǎn)C到其首尾連線(xiàn)的距離d度量桿II 的彈性變形,變形率為d與桿原長(zhǎng)之間的比值.同時(shí),以1944 個(gè)自由度的R3空間ANCF全參數(shù)梁?jiǎn)卧玫降臄?shù)值結(jié)果作為參考值,建模方法數(shù)值結(jié)果收斂標(biāo)準(zhǔn)設(shè)為變形場(chǎng)最大值與參考解的相對(duì)誤差在1%以?xún)?nèi).

6.1.1 模型I

桿I 與桿II 的長(zhǎng)度分別為0.3 m 和0.5 m,兩桿的截面均為矩形,寬與高為0.05 m,材料參數(shù)設(shè)為:ρ=2000 kg/m3,E=8.2×1012Pa,ν=0.3.桿II 末端滑塊的質(zhì)量為0.1 kg.系統(tǒng)初始靜止放置于XOY水平面內(nèi),桿I 與X軸夾角為0?.在曲柄上施加Z方向力矩900 N·m,持續(xù)時(shí)間為0.2 s,桿I 的最大轉(zhuǎn)速約為1033 rad/s.時(shí)間積分算法采用廣義α 方法,仿真時(shí)間為0.2 s,時(shí)間步長(zhǎng)設(shè)為1.0×10?4s,算法參數(shù)譜半徑設(shè)為0.8.

通過(guò)試算,三階SE(3)幾何精確梁84 個(gè)自由度收斂,采用縮減積分技術(shù)和Hu-Washizu 三場(chǎng)變分原理的一階SE(3)幾何精確梁與一階R3×R3幾何精確梁均在156 個(gè)自由度收斂,ANCF 全參數(shù)梁312 個(gè)自由度收斂.圖4 給出了上述三類(lèi)建模方法收斂時(shí)桿II 變形率隨時(shí)間的變化曲線(xiàn).隨著轉(zhuǎn)速增大,桿II 的變形逐漸增大,最大變形能夠達(dá)到桿長(zhǎng)度的0.04%.

6.1.2 模型II

桿I 與桿II 的長(zhǎng)度分別為0.3 m 和0.5 m,兩桿的截面均為矩形,寬與高為0.03 m,材料參數(shù)設(shè)為:ρ=5600 kg/m3,E=4.1×1011Pa,ν=0.3.桿II 末端滑塊質(zhì)量為0.1 kg.系統(tǒng)初始靜止放置于XOY水平面內(nèi),桿I 與X軸夾角為0?.在曲柄上施加Z方向力矩500 N·m,持續(xù)時(shí)間為0.2 s,桿I 的最大轉(zhuǎn)速約為518 rad/s.時(shí)間積分算法采用廣義α 方法,仿真時(shí)間為0.2 s,時(shí)間步長(zhǎng)設(shè)為3.0×10?5s,算法參數(shù)譜半徑設(shè)為0.8.

通過(guò)試算,三階SE(3)幾何精確梁84 個(gè)自由度收斂,采用縮減積分技術(shù)和Hu-Washizu 三場(chǎng)變分原理的一階SE(3)幾何精確梁與一階R3×R3幾何精確梁均在300 個(gè)自由度收斂,ANCF 全參數(shù)梁504 個(gè)自由度收斂.圖5 給出了上述三類(lèi)建模方法收斂時(shí)桿II 變形率隨時(shí)間的變化曲線(xiàn).隨著轉(zhuǎn)速增大,桿II 的變形逐漸增大,最大變形能夠達(dá)到桿長(zhǎng)度的4%.

圖5 桿II 變形率隨時(shí)間變化曲線(xiàn)Fig.5 Deformation rate of the rod II

由圖4 和圖5 可知,SE(3)幾何精確梁、R3×R3幾何精確梁和ANCF 全參數(shù)梁均能達(dá)到收斂的數(shù)值結(jié)果,驗(yàn)證了這三類(lèi)建模方法在描述高轉(zhuǎn)速、小變形以及中等變形多體系統(tǒng)動(dòng)力學(xué)問(wèn)題時(shí)的正確性.采用縮減積分技術(shù)與Hu-Washizu 三場(chǎng)變分原理可以緩解動(dòng)力學(xué)問(wèn)題中的剪切閉鎖,兩種方法的收斂速度無(wú)明顯差別.一階R3×R3幾何精確梁與一階SE(3)幾何精確梁收斂速度無(wú)明顯差別.

在一個(gè)牛頓迭代步中,最耗時(shí)的兩個(gè)部分包括計(jì)算系統(tǒng)剛度矩陣與迭代矩陣的LU 分解/回代.由于本文中的梁?jiǎn)卧捎镁植繕?biāo)架建模方法,剛度矩陣更新次數(shù)大幅降低.因此,可通過(guò)計(jì)算迭代矩陣中非零元素的個(gè)數(shù)來(lái)評(píng)價(jià)一個(gè)牛頓迭代步的計(jì)算速度.對(duì)于模型I,幾類(lèi)一階幾何精確梁迭代矩陣非零元素的個(gè)數(shù)均為2808,三階SE(3)幾何精確梁和ANCF全參數(shù)梁迭代矩陣非零元素的個(gè)數(shù)分別為3528 和11 232.對(duì)于模型II,幾類(lèi)一階幾何精確梁迭代矩陣非零元素的個(gè)數(shù)均為5400,三階SE(3)幾何精確梁和ANCF 全參數(shù)梁迭代矩陣非零元素的個(gè)數(shù)分別為3528 和18 144.由以上數(shù)據(jù)可知,在達(dá)到相同的計(jì)算精度時(shí),幾類(lèi)一階幾何精確梁計(jì)算效率相同,SE(3)和R3×R3幾何精確梁的計(jì)算效率明顯高于ANCF 全參數(shù)梁.

其次,通過(guò)分析系統(tǒng)質(zhì)量矩陣和剛度矩陣更新情況,驗(yàn)證基于局部標(biāo)架的SE(3)幾何精確方法能夠有效地減輕剛體運(yùn)動(dòng)帶來(lái)的非線(xiàn)性問(wèn)題.若忽略外力對(duì)應(yīng)的Jacobian 矩陣,系統(tǒng)迭代矩陣的表達(dá)式為

式中,h為積分步長(zhǎng),β 為廣義α 方法算法參數(shù),M為廣義質(zhì)量矩陣,K為切線(xiàn)剛度矩陣,Φq為約束方程的Jacobian 矩陣.本算例中,牛頓迭代收斂標(biāo)準(zhǔn)設(shè)為廣義坐標(biāo)修正量||?||2<1.0×10?6.表7 給出了一階SE(3)幾何精確梁迭代矩陣的主要部分Jq=h2β(M+K)的更新次數(shù)(number of update,NOU)以及仿真過(guò)程中牛頓迭代總次數(shù)(number of iterations,NOI).表7 中“NNI ≥i”表示一個(gè)時(shí)間步內(nèi),牛頓迭代次數(shù)大于等于i步時(shí),Jq強(qiáng)制更新.

表7 不同模型Jq 的更新次數(shù)及仿真過(guò)程牛頓迭代總次數(shù)Table 7 Number of updating Jq and number of Newton iterations about different models

模型I 的數(shù)值結(jié)果表明:對(duì)于基于局部標(biāo)架的SE(3)幾何精確梁模型,若每個(gè)牛頓迭代步中Jq均更新,則整個(gè)仿真過(guò)程共需要進(jìn)行4713 次牛頓迭代步,同時(shí)Jq也需要更新4712 次.當(dāng)NNI ≥3 時(shí)(即每個(gè)時(shí)間步中牛頓迭代達(dá)到三步時(shí)更新一次Jq),Jq需要更新357 次,總迭代次數(shù)增加了30%左右.當(dāng)NNI≥4 時(shí)(即每個(gè)時(shí)間步中牛頓迭代達(dá)到四步時(shí)更新一次Jq),Jq需要更新10 次,總迭代次數(shù)增加了37%左右.模型II 與模型I 情況類(lèi)似,這里不再贅述.由此可得,對(duì)于此類(lèi)高速轉(zhuǎn)動(dòng)的多體動(dòng)力學(xué)問(wèn)題,基于局部標(biāo)架的SE(3)幾何精確梁能夠極大地降低剛體運(yùn)動(dòng)帶來(lái)的幾何非線(xiàn)性.

6.2 發(fā)生中等變形和超大變形的空間雙擺

本節(jié)對(duì)空間雙擺模型進(jìn)行動(dòng)力學(xué)仿真.如圖6 所示,空間雙擺模型初始放置于XOY水平面上,并在重力作用下開(kāi)始運(yùn)動(dòng).與算例6.1 一致,本算例通過(guò)桿II 中點(diǎn)C到其首尾連線(xiàn)的距離d度量桿II 的彈性變形,變形率為d與桿原長(zhǎng)之間的比值.同時(shí),以1944個(gè)自由度的R3空間ANCF 全參數(shù)梁?jiǎn)卧玫降臄?shù)值結(jié)果作為參考值,建模方法數(shù)值結(jié)果收斂標(biāo)準(zhǔn)設(shè)為變形場(chǎng)最大值與參考解的相對(duì)誤差在1%以?xún)?nèi).

圖6 空間雙擺Fig.6 A spatial double pendulum

6.2.1 模型III

雙擺兩個(gè)擺臂的長(zhǎng)度分別為0.3 m 和0.5 m,兩個(gè)擺臂的截面均為矩形,寬與高為0.5 mm,材料參數(shù)設(shè)為:ρ=2700 kg/m3,E=2.1×1010Pa,ν=0.3,重力加速度為g=9.81 m/s2.時(shí)間積分算法采用廣義α 方法,仿真時(shí)間1 s,時(shí)間步長(zhǎng)設(shè)為5.0×10?4s,算法參數(shù)譜半徑設(shè)為0.8.

通過(guò)試算,三階SE(3)幾何精確梁84 個(gè)自由度收斂,采用縮減積分技術(shù)和Hu-Washizu 三場(chǎng)變分原理的一階SE(3)幾何精確梁均在300 個(gè)自由度收斂,ANCF 全參數(shù)梁1464 個(gè)自由度收斂.300 個(gè)自由度的一階R3×R3幾何精確梁在t=0.457 5 s 發(fā)散.圖7給出了上述三類(lèi)建模方法收斂時(shí)桿II 變形率隨時(shí)間的變化曲線(xiàn),最大變形能夠達(dá)到桿長(zhǎng)度的5.5%.

圖7 桿II 變形率隨時(shí)間變化曲線(xiàn)Fig.7 Deformation rate of the rod II

6.2.2 模型IV

雙擺兩個(gè)擺臂的長(zhǎng)度分別為0.3 m 和0.5 m,兩個(gè)擺臂的截面均為矩形,寬與高為0.3 mm,材料參數(shù)設(shè)為:ρ=2700 kg/m3,E=1.5×1010Pa,ν=0.3,重力加速度為g=9.81 m/s2.時(shí)間積分算法采用廣義α 方法,仿真時(shí)間0.5 s,時(shí)間步長(zhǎng)設(shè)為5.0×10?4s,算法參數(shù)譜半徑設(shè)為0.8.

通過(guò)試算,三階SE(3)幾何精確梁84 個(gè)自由度收斂,采用縮減積分技術(shù)和Hu-Washizu 三場(chǎng)變分原理的一階SE(3)幾何精確梁156 個(gè)自由度收斂,ANCF 全參數(shù)梁1752 個(gè)自由度收斂.156 個(gè)自由度的一階R3×R3幾何精確梁在t=0.444 5 s 發(fā)散.圖8給出了上述三類(lèi)建模方法收斂時(shí)桿II 變形率隨時(shí)間的變化曲線(xiàn),最大變形能夠達(dá)到桿長(zhǎng)度的15%.

由圖7 和圖8 可知,SE(3)幾何精確梁、R3×R3幾何精確梁和ANCF 全參數(shù)梁均能達(dá)到收斂的數(shù)值結(jié)果,驗(yàn)證了這三類(lèi)建模方法在描述中等變形和超大變形多體系統(tǒng)動(dòng)力學(xué)問(wèn)題時(shí)的正確性.采用縮減積分技術(shù)與Hu-Washizu 三場(chǎng)變分原理均可有效緩解動(dòng)力學(xué)中的剪切閉鎖問(wèn)題,兩種方法的收斂速度無(wú)明顯差別.

圖8 桿II 變形率隨時(shí)間變化曲線(xiàn)Fig.8 Deformation rate of the rod II

下面通過(guò)計(jì)算迭代矩陣中非零元素的個(gè)數(shù)來(lái)評(píng)價(jià)一個(gè)牛頓迭代步的計(jì)算速度.對(duì)于模型III,幾類(lèi)一階幾何精確梁迭代矩陣非零元素的個(gè)數(shù)均為5400,三階SE(3)幾何精確梁和ANCF 全參數(shù)梁迭代矩陣非零元素的個(gè)數(shù)分別為3528 和52 704.對(duì)于模型IV,幾類(lèi)一階幾何精確梁迭代矩陣非零元素的個(gè)數(shù)均為2808,三階SE(3)幾何精確梁和ANCF 全參數(shù)梁迭代矩陣非零元素的個(gè)數(shù)分別為3528 和63 072.由以上數(shù)據(jù)可知,在達(dá)到相同的計(jì)算精度時(shí),幾類(lèi)一階幾何精確梁計(jì)算效率相同,SE(3)幾何精確梁的計(jì)算效率明顯高于ANCF 全參數(shù)梁.

通過(guò)與6.1 節(jié)中等變形的平面多體系統(tǒng)對(duì)比,可以發(fā)現(xiàn),隨著柔性體變形的增大,SE(3)幾何精確梁比ANCF 全參數(shù)梁?jiǎn)卧獌?yōu)勢(shì)更明顯.一階R3×R3幾何精確梁在仿真計(jì)算的過(guò)程中發(fā)散,是由于本文僅采用Rescaling 技術(shù)[32]處理轉(zhuǎn)動(dòng)參數(shù),沒(méi)有進(jìn)一步采用其他方法避免轉(zhuǎn)動(dòng)參數(shù)的奇異性.相比較而言,本文提出的局部標(biāo)架建模方法則可自然地避免轉(zhuǎn)動(dòng)參數(shù)的奇異性問(wèn)題.

與6.1 節(jié)統(tǒng)計(jì)方法類(lèi)似,通過(guò)分析系統(tǒng)質(zhì)量矩陣和剛度矩陣更新情況,以說(shuō)明基于局部標(biāo)架的SE(3)幾何精確建模方法描述大轉(zhuǎn)動(dòng)、大變形非線(xiàn)性運(yùn)動(dòng)的能力.本算例中,牛頓迭代收斂標(biāo)準(zhǔn)設(shè)為廣義坐標(biāo)修正量||?||2<1.0×10?6.表8 中給出了一階SE(3)幾何精確梁Jq的更新次數(shù)以及仿真過(guò)程中的牛頓迭代總次數(shù).

模型IV 的數(shù)值結(jié)果表明:對(duì)于基于局部標(biāo)架的SE(3)幾何精確梁模型,若每個(gè)牛頓迭代步中Jq均更新,則整個(gè)仿真過(guò)程共需要進(jìn)行2479 次牛頓迭代步,同時(shí)Jq也需要更新2478 次.當(dāng)NNI ≥3 時(shí)(即每個(gè)時(shí)間步中牛頓迭代達(dá)到三步時(shí)更新一次Jq),Jq需更新440 次,總迭代次數(shù)增加了40%左右.當(dāng)NNI ≥4 時(shí)(即每個(gè)時(shí)間步中牛頓迭代達(dá)到四步時(shí)更新一次Jq),Jq需要更新283 次,總迭代次數(shù)增加了54%左右.模型III 與模型IV 情況類(lèi)似,這里不再贅述.由此可得,對(duì)于此類(lèi)大變形多體動(dòng)力學(xué)問(wèn)題,基于局部標(biāo)架的SE(3)幾何精確梁仍然能夠明顯地降低剛體運(yùn)動(dòng)帶來(lái)的幾何非線(xiàn)性.

表8 不同模型Jq 的更新次數(shù)及仿真過(guò)程牛頓迭代總次數(shù)Table 8 Number of updating Jq and number of Newton iterations about different models

7 結(jié)論

本文基于李群SE(3)局部標(biāo)架,研究幾類(lèi)梁?jiǎn)卧拈]鎖緩解方法.采用Hu-Washizu 三場(chǎng)變分原理緩解了SE(3)幾何精確梁?jiǎn)卧械募羟虚]鎖,采用應(yīng)變分解法緩解了ANCF 全參數(shù)梁?jiǎn)卧械牟此砷]鎖,并對(duì)SE(3)幾何精確梁、R3×R3幾何精確梁和ANCF全參數(shù)梁的單元性能進(jìn)行了算例對(duì)比.靜力學(xué)與動(dòng)力學(xué)算例結(jié)果均表明,緩解閉鎖后的SE(3)幾何精確梁計(jì)算精度高.在進(jìn)行大轉(zhuǎn)動(dòng)、大變形動(dòng)力學(xué)計(jì)算時(shí),R3×R3幾何精確梁?jiǎn)卧獣?huì)出現(xiàn)轉(zhuǎn)動(dòng)參數(shù)奇異問(wèn)題,若處理不當(dāng)則會(huì)不收斂,而SE(3)幾何精確梁?jiǎn)卧杀苊廪D(zhuǎn)動(dòng)參數(shù)奇異性問(wèn)題.SE(3)幾何精確梁的廣義質(zhì)量矩陣主項(xiàng)為常數(shù),與轉(zhuǎn)動(dòng)參數(shù)相關(guān)的切線(xiàn)剛度矩陣滿(mǎn)足剛體轉(zhuǎn)動(dòng)的不變性,在計(jì)算中無(wú)須更新,可極大地減少迭代矩陣的更新次數(shù),有效地提高計(jì)算效率.基于SE(3)群局部標(biāo)架的幾何精確建模方法更適合用于描述梁結(jié)構(gòu)的大轉(zhuǎn)動(dòng)、大變形運(yùn)動(dòng).