兩類基于局部標架梁單元的閉鎖緩解方法1)

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

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

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

引言

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

已有研究表明,在局部標架思想下,可構造出一類SE(3)幾何精確梁單元[11]以及絕對節點坐標梁單元.局部標架梁單元能夠規避剛體運動帶來的幾何非線性問題,離散數值模型中廣義質量矩陣與切線剛度矩陣滿足剛體變換的不變性,可明顯地提高柔性多體系統動力學問題的計算效率.在有限元方法中,梁、板/殼單元普遍存在剪切、薄膜以及泊松閉鎖問題.多年來,有限元領域學者提出了許多單元技術來有效地緩解閉鎖.本文主要關注如何提高局部標架梁單元的收斂性,關于李群局部標架建模與計算方法涉及的其它細節本文不再贅述,可參見作者的相關研究工作.

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

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

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

本文研究局部標架下幾類梁單元的閉鎖緩解方法.采用Hu-Washizu 三場變分原理緩解局部標架下李群SE(3)幾何精確梁的剪切閉鎖,采用應變分解法緩解局部標架下ANCF 全參數梁的泊松閉鎖.通過對比幾何精確與絕對節點坐標兩類建模方法,分析基于局部標架的李群SE(3)幾何精確梁在計算精度與效率方面的優勢.本文中的李群SE(3)幾何精確梁與ANCF 梁單元均基于局部標架,為了方便起見,下文省去“基于局部標架”幾個字.

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

1.1 梁的運動學描述及本構模型

基于經典的Timoshenko 梁假設,李群幾何精確梁理論融合非線性有限元方法與幾何力學思想,在李群SE(3)內建立系統的平衡方程,可精確地描述柔性梁的大轉動、大變形剛柔耦合動力學特性,同時可規避剛體運動帶來的幾何非線性[31].

梁的初始構型和當前構型如圖1 所示.引入參考坐標系{O;e1,e2,e3}.在當前構型中,在梁中線的每一點定義隨體基矢量(t,n,b),引入梁中線位置矢量x和截面旋轉矩陣R,初始構型對應的量為(·)0.

圖1 幾何精確梁初始構型與當前構型Fig.1 Initial and current configurations of the geometrically exact beam

在梁的當前構型中,截面上任意一點p的位置矢量為

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

1.2 SE3 幾何精確梁控制方程及線性化

系統的動力學平衡方程為

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

為獲得一個對稱的切線剛度矩陣,本文在前一個收斂步所在切平面進行物理量的一次與二次變分.因此,應變的變分可以表示為

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

外部虛功可以表示為

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

式中,v=[UTωT]T,U和ω分別為局部標架下的線速度與角速度,頂標·表示變量對時間t的全導數,ρA為單位長度的材料密度,I3表示3×3 的單位矩陣,J為截面形心慣性矩陣.

對內部虛功做線性化可得

應變算符B 的線性化的轉置為

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

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

最后,內部虛功的線性化可以寫為

需要指出的是,由于δh為局部標架下的平動與轉動虛位移,?h為局部標架下的平動與轉動無限小位移,它們均與該點的任意剛體位移無關.因此,單元的彈性力及其切線剛度矩陣均滿足剛體變換的不變性.對于廣義慣性力及其廣義質量矩陣也同時具有類似的形式.

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

對梁中線位置矢量x和截面相對旋轉偽矢量Θr[11]進行空間離散

式中,NI是節點I的形函數,xI是梁中線上節點I的位置矢量,Rr是參考點所在截面的旋轉矩陣,RI是節點I所在截面的旋轉矩陣,旋轉偽矢量為轉軸單位矢量與轉角的乘積.

離散應變算符B 得到應變?位移矩陣BI,它的表達式為

彈性力可以表示為

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

對于非常細長的梁結構,采用一階插值描述幾何精確梁彎曲時將會產生剪切閉鎖.由于梁非常細長,截面內橫向剪切應變處處接近于0.而由于平動參數x和轉動參數θ采用了相同的插值函數,使式(5)中橫向剪切應變γs中的R(θ)和x′兩項具有不同階次,γs=0 不可能處處滿足,導致能量泛函中剪切變形能項的量級不正確,由此帶來了剪切閉鎖.此時,單元表現得十分剛硬,即計算出來的位移值遠小于參考解,總體剛度矩陣很大.

2 基于局部標架的絕對節點坐標方法的三維全參數梁理論

本節簡要介紹基于局部標架的絕對節點坐標全參數梁單元建模方法.某質點從初始位置X運動到當前位置r,變形梯度可以表示為

Green-Lagrange 應變張量可以表示為

應變能可以表示為

式中,D為彈性矩陣,V為梁的體積.應變能對廣義坐標e求偏導得彈性力為

彈性力對廣義坐標求偏導得剛度陣為

相比于幾何精確梁單元,ANCF 單元沒有轉動參數,不能直接定義局部標架,需要通過平動位移場定義轉動場.此處采用ANCF 斜率矢量定義該點的局部標架基矢量

式中,rX和rY為當前位置矢量r對物質坐標的偏導數.由于平動位移場在局部標架中描述,單元的剛體運動可自然被消除.因此,系統質量矩陣以及切線剛度矩陣滿足剛體運動的不變性,從而可極大地減少系統質量矩陣與剛度矩陣在仿真過程中的更新次數.

ANCF 全參數梁單元會遇到泊松閉鎖問題,閉鎖產生的原因如下.由于ANCF 全參數梁是連續介質單元,梁的中線為三階Hermite 插值,在泊松效應的作用下,截面應該變形為曲面.然而橫截面為一階插值,不能描述截面的曲面變形模式,因此導致泊松閉鎖.

3 兩類局部標架梁單元的閉鎖緩解方法

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

Hu-Washizu 三場變分原理不僅假設位移場,同時假設應力和應變場,且三類場變量的插值相互獨立.假設剪切應變和假設剪切應力的分布為

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

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

式中,

離散方程組可以表示為

式中,Km和Kb分別為薄膜和彎曲部分的剛度矩陣,?hN為節點運動矢量的增量,Res為殘差.將上式靜力縮聚化為單場形式

一階插值的幾何精確梁會產生嚴重的剪切閉鎖,二階及二階以上插值的幾何精確梁剪切閉鎖較輕.因此,本文關于幾何精確梁剪切閉鎖處理方法著重于一階梁單元.Hu-Washizu 三場變分原理中,位移場、應變場和應力場均獨立插值.式(28)假設橫向剪切應變為常數,可避免產生過大的剪切應變能,有效緩解了剪切閉鎖.

此外,對于梁單元,縮減積分技術也是緩解剪切閉鎖的有效手段.數值算例展示了采用Hu-Washizu三場變分原理與縮減積分技術緩解梁單元剪切閉鎖的效果.相比較而言,縮減積分技術由于其簡便性在梁單元中應用較為廣泛,但用在板殼單元中會產生沙漏變形模式.進一步,上述Hu-Washizu 三場變分原理的研究將為低階幾何精確板殼單元的閉鎖緩解起到一定的借鑒作用.

3.2 基于應變分解法的泊松閉鎖緩解方法

ANCF 全參數梁單元的位置場和變形梯度可以寫為

εc中只有低階應變,且εc僅與梁中線參數X相關;εk中包含高階應變,這些高階應變會導致截面變形和彎曲變形.為了緩解閉鎖,泊松耦合只用在正應變的低階項中,即僅在εc中考慮泊松效應,忽略正應變的高階項的泊松效應.因此,梁截面不會變形為曲面,緩解了泊松閉鎖.第二類Piola-Kirchhoff應力的Voigt形式可以寫為

式中,應變也是Voigt 形式.基于平面應變假設,彈性系數矩陣可以被定義為

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

彈性力可以表示為

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

4 三類常用的有限元插值方法

本節簡要介紹三類有限元插值方法,包括Lagrange 插值、Hermite 插值以及非均勻有理B-樣條(NURBS).其中,Lagrange 插值在有限元中應用最為廣泛;絕對節點坐標方法則采用Hermite 插值,以保證位移場的C1連續;NURBS 的基函數可以構造任意階連續的近似函數,建立了計算機輔助設計(CAD)與計算機輔助工程(CAE)間的橋梁,該方法被稱為等幾何分析.

為對比分析不同插值方法的優劣勢,本文涉及一階與三階Lagrange 插值、三階Hermite 插值以及三階NURBS 插值.

4.1 Lagrange 插值

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

4.2 Hermite 插值

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

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

4.3 非均勻有理B-樣條

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

式中,Pi是控制點,ωi是權系數.

5 直/曲梁結構大變形靜力學分析

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

5.1 懸臂直梁和曲梁低階單元閉鎖緩解方法數值測試

本節考查SE(3)幾何精確一階直梁與曲梁單元在三維復雜應力狀態下,縮減積分技術和Hu-Washizu三場變分原理緩解剪切閉鎖的能力.

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

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

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

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

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

5.2 懸臂直梁和曲梁純彎曲測試

本節考查SE(3)幾何精確一階直梁與曲梁單元在純彎曲應力狀態下,縮減積分技術和Hu-Washizu三場變分原理緩解剪切閉鎖的能力,同時評價單元描述純彎曲的能力.

在懸臂直梁和曲梁末端施加集中彎矩,此時,懸臂梁發生純彎曲,處于簡單應力狀態.單元位移場采用一階Lagrange 插值.對比分析完全積分、縮減積分和Hu-Washizu 三場變分原理三種情況下,懸臂梁末端點z方向的轉角位移θz,計算相對誤差,將結果列于表3 和表4 中.同時,采用ABAQUS 軟件中B31 梁單元對本算例進行仿真計算.

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

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

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

同時,采用三階Lagrange 插值構造了高階幾何精確直梁與曲梁單元.通過上述懸臂直梁與曲梁算例,測試了單元收斂精度,其中同樣采用Hu-Washizu三場變分原理處理剪切閉鎖.數值結果表明:處理閉鎖不會提升單元收斂性能,進一步說明了高階直梁與曲梁單元不存在剪切閉鎖.在5.3 節中,將進一步對比處理閉鎖后的一階單元與高階單元的收斂性.

5.3 幾類局部標架梁單元收斂性對比分析

本節對比一階和三階SE(3)幾何精確梁單元、一階R3×R3幾何精確梁、ANCF 縮減(slope deficiency)梁、ANCF 全參數梁與ABAQUS 軟件中B31 梁單元的收斂性.其中,一階幾何精確梁與ANCF 全參數梁單元均已處理閉鎖.

在懸臂直梁和曲梁末端施加集中力.計算以上幾類建模方法在不同自由度下懸臂梁末端點z方向的位移uz和相對誤差,將結果列于表5 和表6 中.ANCF 縮減曲梁單元無算進行本算例的計算。

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

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

表5 和表6 的數值結果表明:直梁模型的五種建模方法收斂結果基本一致.對于曲梁模型,SE(3)幾何精確梁、R3×R3幾何精確梁和ABAQUS 中B31 單元收斂結果基本一致,ANCF 全參數梁單元收斂結果略大于其他三類方法的收斂值.需要指出的是,表5 中三階NURBS 插值的SE(3)幾何精確梁單元在24 個自由度時,雖然計算結果的相對誤差為6.6×10?4,但此時單元還沒有收斂,60 個自由度時,單元才收斂.對于高階單元,三階Lagrange 比三階NURBS 插值的SE(3)幾何精確梁的計算精度高,幾類幾何精確梁均比ANCF 梁單元的計算精度高.ANCF 縮減梁比ANCF 全參數梁單元的計算精度高.

對于一階直梁單元,兩類幾何精確梁具有完全相同的計算精度.對于一階曲梁單元,R3×R3幾何精確梁的計算精度略高于SE(3)幾何精確梁.對于一階直梁和曲梁單元,兩類幾何精確方法均比ABAQUS中B31 單元的計算精度高.顯然高階單元比低階單元精度高,但是高階單元計算效率低,即自由度相同時,高階單元剛度矩陣內非零元素的個數多.

6 高轉速/大變形柔性多體系統動力學驗證算例

本節通過高轉速曲柄滑塊機構以及空間雙擺的動力學仿真,分析局部標架下的梁單元在動力學仿真中的計算精度差異,以及消除剛體運動帶來的幾何非線性后梁單元的數值特性.由于商業軟件無法直接進行大變形柔性多體系統動力學仿真,因此本節不再采用商業軟件進行仿真對比.為了設計特定的小變形或大變形的柔性多體系統動力學模型,本節算例中部分材料參數未與真實材料對應.

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

6.1 發生小變形和中等變形的高速旋轉的平面曲柄滑塊機構

本節對一個典型的多體系統——曲柄滑塊機構進行動力學仿真.如圖3 所示,平面曲柄滑塊機構由兩根桿與一個滑塊組成,桿I 一端與地面采用球鉸相連接,另一端與桿II 一端也通過球鉸連接,桿II 末端與一個僅能在X軸上運動的剛體滑塊鉸接.

圖3 曲柄滑塊機構Fig.3 A crank-slider mechanism

本算例中,通過桿II 中點C到其首尾連線的距離d度量桿II 的彈性變形,變形率為d與桿原長之間的比值.同時,以1944 個自由度的R3空間ANCF全參數梁單元所得到的數值結果作為參考值,建模方法數值結果收斂標準設為變形場最大值與參考解的相對誤差在1%以內.

6.1.1 模型I

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

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

6.1.2 模型II

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

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

圖5 桿II 變形率隨時間變化曲線Fig.5 Deformation rate of the rod II

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

在一個牛頓迭代步中,最耗時的兩個部分包括計算系統剛度矩陣與迭代矩陣的LU 分解/回代.由于本文中的梁單元均采用局部標架建模方法,剛度矩陣更新次數大幅降低.因此,可通過計算迭代矩陣中非零元素的個數來評價一個牛頓迭代步的計算速度.對于模型I,幾類一階幾何精確梁迭代矩陣非零元素的個數均為2808,三階SE(3)幾何精確梁和ANCF全參數梁迭代矩陣非零元素的個數分別為3528 和11 232.對于模型II,幾類一階幾何精確梁迭代矩陣非零元素的個數均為5400,三階SE(3)幾何精確梁和ANCF 全參數梁迭代矩陣非零元素的個數分別為3528 和18 144.由以上數據可知,在達到相同的計算精度時,幾類一階幾何精確梁計算效率相同,SE(3)和R3×R3幾何精確梁的計算效率明顯高于ANCF 全參數梁.

其次,通過分析系統質量矩陣和剛度矩陣更新情況,驗證基于局部標架的SE(3)幾何精確方法能夠有效地減輕剛體運動帶來的非線性問題.若忽略外力對應的Jacobian 矩陣,系統迭代矩陣的表達式為

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

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

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

6.2 發生中等變形和超大變形的空間雙擺

本節對空間雙擺模型進行動力學仿真.如圖6 所示,空間雙擺模型初始放置于XOY水平面上,并在重力作用下開始運動.與算例6.1 一致,本算例通過桿II 中點C到其首尾連線的距離d度量桿II 的彈性變形,變形率為d與桿原長之間的比值.同時,以1944個自由度的R3空間ANCF 全參數梁單元所得到的數值結果作為參考值,建模方法數值結果收斂標準設為變形場最大值與參考解的相對誤差在1%以內.

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

6.2.1 模型III

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

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

圖7 桿II 變形率隨時間變化曲線Fig.7 Deformation rate of the rod II

6.2.2 模型IV

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

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

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

圖8 桿II 變形率隨時間變化曲線Fig.8 Deformation rate of the rod II

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

通過與6.1 節中等變形的平面多體系統對比,可以發現,隨著柔性體變形的增大,SE(3)幾何精確梁比ANCF 全參數梁單元優勢更明顯.一階R3×R3幾何精確梁在仿真計算的過程中發散,是由于本文僅采用Rescaling 技術[32]處理轉動參數,沒有進一步采用其他方法避免轉動參數的奇異性.相比較而言,本文提出的局部標架建模方法則可自然地避免轉動參數的奇異性問題.

與6.1 節統計方法類似,通過分析系統質量矩陣和剛度矩陣更新情況,以說明基于局部標架的SE(3)幾何精確建模方法描述大轉動、大變形非線性運動的能力.本算例中,牛頓迭代收斂標準設為廣義坐標修正量||?||2<1.0×10?6.表8 中給出了一階SE(3)幾何精確梁Jq的更新次數以及仿真過程中的牛頓迭代總次數.

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

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

7 結論

本文基于李群SE(3)局部標架,研究幾類梁單元的閉鎖緩解方法.采用Hu-Washizu 三場變分原理緩解了SE(3)幾何精確梁單元中的剪切閉鎖,采用應變分解法緩解了ANCF 全參數梁單元中的泊松閉鎖,并對SE(3)幾何精確梁、R3×R3幾何精確梁和ANCF全參數梁的單元性能進行了算例對比.靜力學與動力學算例結果均表明,緩解閉鎖后的SE(3)幾何精確梁計算精度高.在進行大轉動、大變形動力學計算時,R3×R3幾何精確梁單元會出現轉動參數奇異問題,若處理不當則會不收斂,而SE(3)幾何精確梁單元可避免轉動參數奇異性問題.SE(3)幾何精確梁的廣義質量矩陣主項為常數,與轉動參數相關的切線剛度矩陣滿足剛體轉動的不變性,在計算中無須更新,可極大地減少迭代矩陣的更新次數,有效地提高計算效率.基于SE(3)群局部標架的幾何精確建模方法更適合用于描述梁結構的大轉動、大變形運動.